1 static char help[] = "Test DMStag ghosted boundaries in 3d\n\n";
2 /* This solves a very contrived problem - the "pressure" terms are set to a constant function
3 and the "velocity" terms are just the sum of neighboring values of these, hence twice the
4 constant */
5 #include <petscdm.h>
6 #include <petscksp.h>
7 #include <petscdmstag.h>
8
9 #define PRESSURE_CONST 2.0
10
11 PetscErrorCode ApplyOperator(Mat, Vec, Vec);
12
main(int argc,char ** argv)13 int main(int argc, char **argv)
14 {
15 DM dmSol;
16 Vec sol, solRef, solRefLocal, rhs, rhsLocal;
17 Mat A;
18 KSP ksp;
19 PC pc;
20 PetscInt startx, starty, startz, nx, ny, nz, ex, ey, ez, nExtrax, nExtray, nExtraz;
21 PetscInt iux, iuy, iuz, ip;
22 PetscInt dof0, dof1, dof2, dof3;
23 PetscScalar ****arrSol, ****arrRHS;
24
25 PetscFunctionBeginUser;
26 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
27 /* Note: these defaults are chosen to suit the problem. We allow adjusting
28 them to check that things still work when you add unused extra dof */
29 dof0 = 0;
30 dof1 = 0;
31 dof2 = 1;
32 dof3 = 1;
33 PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, 3, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, 1, NULL, NULL, NULL, &dmSol));
34 PetscCall(DMSetFromOptions(dmSol));
35 PetscCall(DMSetUp(dmSol));
36
37 /* Compute reference solution on the grid, using direct array access */
38 PetscCall(DMCreateGlobalVector(dmSol, &rhs));
39 PetscCall(DMCreateGlobalVector(dmSol, &solRef));
40 PetscCall(DMGetLocalVector(dmSol, &solRefLocal));
41 PetscCall(DMGetLocalVector(dmSol, &rhsLocal));
42 PetscCall(DMStagVecGetArray(dmSol, solRefLocal, &arrSol));
43
44 PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, &nExtrax, &nExtray, &nExtraz));
45 PetscCall(DMStagVecGetArray(dmSol, rhsLocal, &arrRHS));
46
47 /* Get the correct entries for each of our variables in local element-wise storage */
48 PetscCall(DMStagGetLocationSlot(dmSol, DMSTAG_LEFT, 0, &iux));
49 PetscCall(DMStagGetLocationSlot(dmSol, DMSTAG_DOWN, 0, &iuy));
50 PetscCall(DMStagGetLocationSlot(dmSol, DMSTAG_BACK, 0, &iuz));
51 PetscCall(DMStagGetLocationSlot(dmSol, DMSTAG_ELEMENT, 0, &ip));
52 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
53 for (ey = starty; ey < starty + ny + nExtray; ++ey) {
54 for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
55 arrSol[ez][ey][ex][iux] = 2 * PRESSURE_CONST;
56 arrRHS[ez][ey][ex][iux] = 0.0;
57 arrSol[ez][ey][ex][iuy] = 2 * PRESSURE_CONST;
58 arrRHS[ez][ey][ex][iuy] = 0.0;
59 arrSol[ez][ey][ex][iuz] = 2 * PRESSURE_CONST;
60 arrRHS[ez][ey][ex][iuz] = 0.0;
61 if (ex < startx + nx && ey < starty + ny && ez < startz + nz) {
62 arrSol[ez][ey][ex][ip] = PRESSURE_CONST;
63 arrRHS[ez][ey][ex][ip] = PRESSURE_CONST;
64 }
65 }
66 }
67 }
68 PetscCall(DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS));
69 PetscCall(DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs));
70 PetscCall(DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs));
71 PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
72 PetscCall(DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef));
73 PetscCall(DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef));
74 PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
75 PetscCall(DMRestoreLocalVector(dmSol, &rhsLocal));
76
77 /* Matrix-free Operator */
78 PetscCall(DMSetMatType(dmSol, MATSHELL));
79 PetscCall(DMCreateMatrix(dmSol, &A));
80 PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)ApplyOperator));
81
82 /* Solve */
83 PetscCall(DMCreateGlobalVector(dmSol, &sol));
84 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
85 PetscCall(KSPSetOperators(ksp, A, A));
86 PetscCall(KSPGetPC(ksp, &pc));
87 PetscCall(PCSetType(pc, PCNONE));
88 PetscCall(KSPSetFromOptions(ksp));
89 PetscCall(KSPSolve(ksp, rhs, sol));
90
91 /* Check Solution */
92 {
93 Vec diff;
94 PetscReal normsolRef, errAbs, errRel;
95
96 PetscCall(VecDuplicate(sol, &diff));
97 PetscCall(VecCopy(sol, diff));
98 PetscCall(VecAXPY(diff, -1.0, solRef));
99 PetscCall(VecNorm(diff, NORM_2, &errAbs));
100 PetscCall(VecNorm(solRef, NORM_2, &normsolRef));
101 errRel = errAbs / normsolRef;
102 if (errAbs > 1e14 || errRel > 1e14) {
103 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel));
104 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Non-zero error. Probable failure.\n"));
105 }
106 PetscCall(VecDestroy(&diff));
107 }
108
109 PetscCall(KSPDestroy(&ksp));
110 PetscCall(VecDestroy(&sol));
111 PetscCall(VecDestroy(&solRef));
112 PetscCall(VecDestroy(&rhs));
113 PetscCall(MatDestroy(&A));
114 PetscCall(DMDestroy(&dmSol));
115 PetscCall(PetscFinalize());
116 return 0;
117 }
118
ApplyOperator(Mat A,Vec in,Vec out)119 PetscErrorCode ApplyOperator(Mat A, Vec in, Vec out)
120 {
121 DM dm;
122 Vec inLocal, outLocal;
123 PetscScalar ****arrIn;
124 PetscScalar ****arrOut;
125 PetscInt startx, starty, startz, nx, ny, nz, nExtrax, nExtray, nExtraz, ex, ey, ez, idxP, idxUx, idxUy, idxUz, startGhostx, startGhosty, startGhostz, nGhostx, nGhosty, nGhostz;
126 PetscBool isFirstx, isFirsty, isFirstz, isLastx, isLasty, isLastz;
127
128 PetscFunctionBeginUser;
129 PetscCall(MatGetDM(A, &dm));
130 PetscCall(DMGetLocalVector(dm, &inLocal));
131 PetscCall(DMGetLocalVector(dm, &outLocal));
132 PetscCall(DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal));
133 PetscCall(DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal));
134 PetscCall(DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, &nExtrax, &nExtray, &nExtraz));
135 PetscCall(DMStagGetGhostCorners(dm, &startGhostx, &startGhosty, &startGhostz, &nGhostx, &nGhosty, &nGhostz));
136 PetscCall(DMStagVecGetArrayRead(dm, inLocal, &arrIn));
137 PetscCall(DMStagVecGetArray(dm, outLocal, &arrOut));
138 PetscCall(DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &idxUx));
139 PetscCall(DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &idxUy));
140 PetscCall(DMStagGetLocationSlot(dm, DMSTAG_BACK, 0, &idxUz));
141 PetscCall(DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &idxP));
142 PetscCall(DMStagGetIsFirstRank(dm, &isFirstx, &isFirsty, &isFirstz));
143 PetscCall(DMStagGetIsLastRank(dm, &isLastx, &isLasty, &isLastz));
144
145 /* Set "pressures" on ghost boundaries by copying neighboring values*/
146 if (isFirstx) {
147 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
148 for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ez][ey][-1][idxP] = arrIn[ez][ey][0][idxP];
149 }
150 }
151 if (isLastx) {
152 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
153 for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ez][ey][startx + nx][idxP] = arrIn[ez][ey][startx + nx - 1][idxP];
154 }
155 }
156 if (isFirsty) {
157 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
158 for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[ez][-1][ex][idxP] = arrIn[ez][0][ex][idxP];
159 }
160 }
161 if (isLasty) {
162 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
163 for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[ez][starty + ny][ex][idxP] = arrIn[ez][starty + ny - 1][ex][idxP];
164 }
165 }
166
167 if (isFirstz) {
168 for (ey = starty; ey < starty + ny + nExtray; ++ey) {
169 for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[-1][ey][ex][idxP] = arrIn[0][ey][ex][idxP];
170 }
171 }
172 if (isLastz) {
173 for (ey = starty; ey < starty + ny + nExtray; ++ey) {
174 for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[startz + nz][ey][ex][idxP] = arrIn[startz + nz - 1][ey][ex][idxP];
175 }
176 }
177
178 /* Apply operator on physical points */
179 for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
180 for (ey = starty; ey < starty + ny + nExtray; ++ey) {
181 for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
182 if (ex < startx + nx && ey < starty + ny && ez < startz + nz) { /* Don't compute pressure outside domain */
183 arrOut[ez][ey][ex][idxP] = arrIn[ez][ey][ex][idxP];
184 }
185 arrOut[ez][ey][ex][idxUx] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey][ex - 1][idxP] - arrIn[ez][ey][ex][idxUx];
186 arrOut[ez][ey][ex][idxUy] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey - 1][ex][idxP] - arrIn[ez][ey][ex][idxUy];
187 arrOut[ez][ey][ex][idxUz] = arrIn[ez][ey][ex][idxP] + arrIn[ez - 1][ey][ex][idxP] - arrIn[ez][ey][ex][idxUz];
188 }
189 }
190 }
191 PetscCall(DMStagVecRestoreArrayRead(dm, inLocal, &arrIn));
192 PetscCall(DMStagVecRestoreArray(dm, outLocal, &arrOut));
193 PetscCall(DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out));
194 PetscCall(DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out));
195 PetscCall(DMRestoreLocalVector(dm, &inLocal));
196 PetscCall(DMRestoreLocalVector(dm, &outLocal));
197 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
200 /*TEST
201
202 test:
203 suffix: 1
204 nsize: 1
205 output_file: output/empty.out
206
207 test:
208 suffix: 2
209 nsize: 8
210 output_file: output/empty.out
211
212 test:
213 suffix: 3
214 nsize: 1
215 args: -stag_stencil_width 2
216 output_file: output/empty.out
217
218 test:
219 suffix: 4
220 nsize: 8
221 args: -stag_grid_x 4 -stag_grid_y 5 -stag_grid_z 4 -stag_stencil_width 2
222 output_file: output/empty.out
223
224 test:
225 suffix: 5
226 nsize: 8
227 args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6 -stag_grid_z 6
228 output_file: output/empty.out
229
230 TEST*/
231