Actual source code: ex8.c
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>
9: #define PRESSURE_CONST 2.0
11: PetscErrorCode ApplyOperator(Mat, Vec, Vec);
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;
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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));
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));
44: PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, &nExtrax, &nExtray, &nExtraz));
45: PetscCall(DMStagVecGetArray(dmSol, rhsLocal, &arrRHS));
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));
77: /* Matrix-free Operator */
78: PetscCall(DMSetMatType(dmSol, MATSHELL));
79: PetscCall(DMCreateMatrix(dmSol, &A));
80: PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))ApplyOperator));
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));
91: /* Check Solution */
92: {
93: Vec diff;
94: PetscReal normsolRef, errAbs, errRel;
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: }
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: }
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;
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));
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: }
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: }
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: }
200: /*TEST
202: test:
203: suffix: 1
204: nsize: 1
206: test:
207: suffix: 2
208: nsize: 8
210: test:
211: suffix: 3
212: nsize: 1
213: args: -stag_stencil_width 2
215: test:
216: suffix: 4
217: nsize: 8
218: args: -stag_grid_x 4 -stag_grid_y 5 -stag_grid_z 4 -stag_stencil_width 2
220: test:
221: suffix: 5
222: nsize: 8
223: 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
225: TEST*/