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, 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));

 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*/