Actual source code: ex11.c

  1: static char help[] = "Test DMStag ghosted boundaries in 2d\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, nx, ny, ex, ey, nExtrax, nExtray;
 21:   PetscInt       iux, iuy, ip;
 22:   PetscInt       dof0, dof1, dof2;
 23:   PetscScalar ***arrSol, ***arrRHS;

 26:   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 = 1;
 31:   dof2 = 1;
 32:   PetscOptionsGetInt(NULL, NULL, "-dof2", &dof2, NULL);
 33:   DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, 3, 3, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dmSol);
 34:   DMSetFromOptions(dmSol);
 35:   DMSetUp(dmSol);

 37:   /* Compute reference solution on the grid, using direct array access */
 38:   DMCreateGlobalVector(dmSol, &rhs);
 39:   DMCreateGlobalVector(dmSol, &solRef);
 40:   DMGetLocalVector(dmSol, &solRefLocal);
 41:   DMGetLocalVector(dmSol, &rhsLocal);
 42:   DMStagVecGetArray(dmSol, solRefLocal, &arrSol);

 44:   DMStagGetCorners(dmSol, &startx, &starty, NULL, &nx, &ny, NULL, &nExtrax, &nExtray, NULL);
 45:   DMStagVecGetArray(dmSol, rhsLocal, &arrRHS);

 47:   /* Get the correct entries for each of our variables in local element-wise storage */
 48:   DMStagGetLocationSlot(dmSol, DMSTAG_LEFT, 0, &iux);
 49:   DMStagGetLocationSlot(dmSol, DMSTAG_DOWN, 0, &iuy);
 50:   DMStagGetLocationSlot(dmSol, DMSTAG_ELEMENT, 0, &ip);
 51:   for (ey = starty; ey < starty + ny + nExtray; ++ey) {
 52:     for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
 53:       arrSol[ey][ex][iux] = 2 * PRESSURE_CONST;
 54:       arrRHS[ey][ex][iux] = 0.0;
 55:       arrSol[ey][ex][iuy] = 2 * PRESSURE_CONST;
 56:       arrRHS[ey][ex][iuy] = 0.0;
 57:       if (ex < startx + nx && ey < starty + ny) {
 58:         arrSol[ey][ex][ip] = PRESSURE_CONST;
 59:         arrRHS[ey][ex][ip] = PRESSURE_CONST;
 60:       }
 61:     }
 62:   }
 63:   DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS);
 64:   DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs);
 65:   DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs);
 66:   DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol);
 67:   DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef);
 68:   DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef);
 69:   DMRestoreLocalVector(dmSol, &solRefLocal);
 70:   DMRestoreLocalVector(dmSol, &rhsLocal);

 72:   /* Matrix-free Operator */
 73:   DMSetMatType(dmSol, MATSHELL);
 74:   DMCreateMatrix(dmSol, &A);
 75:   MatShellSetOperation(A, MATOP_MULT, (void (*)(void))ApplyOperator);

 77:   /* Solve */
 78:   DMCreateGlobalVector(dmSol, &sol);
 79:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 80:   KSPSetOperators(ksp, A, A);
 81:   KSPGetPC(ksp, &pc);
 82:   PCSetType(pc, PCNONE);
 83:   KSPSetFromOptions(ksp);
 84:   KSPSolve(ksp, rhs, sol);

 86:   /* Check Solution */
 87:   {
 88:     Vec       diff;
 89:     PetscReal normsolRef, errAbs, errRel;

 91:     VecDuplicate(sol, &diff);
 92:     VecCopy(sol, diff);
 93:     VecAXPY(diff, -1.0, solRef);
 94:     VecNorm(diff, NORM_2, &errAbs);
 95:     VecNorm(solRef, NORM_2, &normsolRef);
 96:     errRel = errAbs / normsolRef;
 97:     if (errAbs > 1e14 || errRel > 1e14) {
 98:       PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel);
 99:       PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Non-zero error. Probable failure.\n");
100:     }
101:     VecDestroy(&diff);
102:   }

104:   KSPDestroy(&ksp);
105:   VecDestroy(&sol);
106:   VecDestroy(&solRef);
107:   VecDestroy(&rhs);
108:   MatDestroy(&A);
109:   DMDestroy(&dmSol);
110:   PetscFinalize();
111:   return 0;
112: }

114: PetscErrorCode ApplyOperator(Mat A, Vec in, Vec out)
115: {
116:   DM             dm;
117:   Vec            inLocal, outLocal;
118:   PetscScalar ***arrIn;
119:   PetscScalar ***arrOut;
120:   PetscInt       startx, starty, nx, ny, nExtrax, nExtray, ex, ey, idxP, idxUx, idxUy, startGhostx, startGhosty, nGhostx, nGhosty;
121:   PetscBool      isFirstx, isFirsty, isFirstz, isLastx, isLasty, isLastz;

124:   MatGetDM(A, &dm);
125:   DMGetLocalVector(dm, &inLocal);
126:   DMGetLocalVector(dm, &outLocal);
127:   DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal);
128:   DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal);
129:   DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, &nExtrax, &nExtray, NULL);
130:   DMStagGetGhostCorners(dm, &startGhostx, &startGhosty, NULL, &nGhostx, &nGhosty, NULL);
131:   DMStagVecGetArrayRead(dm, inLocal, &arrIn);
132:   DMStagVecGetArray(dm, outLocal, &arrOut);
133:   DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &idxUx);
134:   DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &idxUy);
135:   DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &idxP);
136:   DMStagGetIsFirstRank(dm, &isFirstx, &isFirsty, &isFirstz);
137:   DMStagGetIsLastRank(dm, &isLastx, &isLasty, &isLastz);

139:   /* Set "pressures" on ghost boundaries by copying neighboring values*/
140:   if (isFirstx) {
141:     for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ey][-1][idxP] = arrIn[ey][0][idxP];
142:   }
143:   if (isLastx) {
144:     for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ey][startx + nx][idxP] = arrIn[ey][startx + nx - 1][idxP];
145:   }
146:   if (isFirsty) {
147:     for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[-1][ex][idxP] = arrIn[0][ex][idxP];
148:   }
149:   if (isLasty) {
150:     for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[starty + ny][ex][idxP] = arrIn[starty + ny - 1][ex][idxP];
151:   }

153:   /* Apply operator on physical points */
154:   for (ey = starty; ey < starty + ny + nExtray; ++ey) {
155:     for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
156:       if (ex < startx + nx && ey < starty + ny) { /* Don't compute pressure outside domain */
157:         arrOut[ey][ex][idxP] = arrIn[ey][ex][idxP];
158:       }
159:       arrOut[ey][ex][idxUx] = arrIn[ey][ex][idxP] + arrIn[ey][ex - 1][idxP] - arrIn[ey][ex][idxUx];
160:       arrOut[ey][ex][idxUy] = arrIn[ey][ex][idxP] + arrIn[ey - 1][ex][idxP] - arrIn[ey][ex][idxUy];
161:     }
162:   }
163:   DMStagVecRestoreArrayRead(dm, inLocal, &arrIn);
164:   DMStagVecRestoreArray(dm, outLocal, &arrOut);
165:   DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out);
166:   DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out);
167:   DMRestoreLocalVector(dm, &inLocal);
168:   DMRestoreLocalVector(dm, &outLocal);
169:   return 0;
170: }

172: /*TEST

174:    test:
175:       suffix: 1
176:       nsize: 1

178:    test:
179:       suffix: 2
180:       nsize: 4

182:    test:
183:       suffix: 3
184:       nsize: 1
185:       args: -stag_stencil_width 2

187:    test:
188:       suffix: 4
189:       nsize: 4
190:       args: -stag_grid_x 4 -stag_grid_y 5 -stag_stencil_width 2

192:    test:
193:       suffix: 5
194:       nsize: 4
195:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6

197: TEST*/