Actual source code: ex5.c

  1: static char help[] = "Test DMStag ghosted boundaries in 1d\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 LEFT    DMSTAG_LEFT
 10: #define RIGHT   DMSTAG_RIGHT
 11: #define ELEMENT DMSTAG_ELEMENT

 13: #define PRESSURE_CONST 2.0

 15: PetscErrorCode ApplyOperator(Mat, Vec, Vec);

 17: int main(int argc, char **argv)
 18: {
 19:   DM            dmSol;
 20:   Vec           sol, solRef, solRefLocal, rhs, rhsLocal;
 21:   Mat           A;
 22:   KSP           ksp;
 23:   PC            pc;
 24:   PetscInt      start, n, e, nExtra;
 25:   PetscInt      iu, ip;
 26:   PetscScalar **arrSol, **arrRHS;

 28:   PetscFunctionBeginUser;
 29:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 30:   PetscCall(DMStagCreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, 3, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, &dmSol));
 31:   PetscCall(DMSetFromOptions(dmSol));
 32:   PetscCall(DMSetUp(dmSol));

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

 41:   PetscCall(DMStagGetCorners(dmSol, &start, NULL, NULL, &n, NULL, NULL, &nExtra, NULL, NULL));
 42:   PetscCall(DMStagVecGetArray(dmSol, rhsLocal, &arrRHS));

 44:   /* Get the correct entries for each of our variables in local element-wise storage */
 45:   PetscCall(DMStagGetLocationSlot(dmSol, LEFT, 0, &iu));
 46:   PetscCall(DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip));
 47:   for (e = start; e < start + n + nExtra; ++e) {
 48:     {
 49:       arrSol[e][iu] = 2 * PRESSURE_CONST;
 50:       arrRHS[e][iu] = 0.0;
 51:     }
 52:     if (e < start + n) {
 53:       arrSol[e][ip] = PRESSURE_CONST;
 54:       arrRHS[e][ip] = PRESSURE_CONST;
 55:     }
 56:   }
 57:   PetscCall(DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS));
 58:   PetscCall(DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs));
 59:   PetscCall(DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs));
 60:   PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
 61:   PetscCall(DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef));
 62:   PetscCall(DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef));
 63:   PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
 64:   PetscCall(DMRestoreLocalVector(dmSol, &rhsLocal));

 66:   /* Matrix-free Operator */
 67:   PetscCall(DMSetMatType(dmSol, MATSHELL));
 68:   PetscCall(DMCreateMatrix(dmSol, &A));
 69:   PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))ApplyOperator));

 71:   /* Solve */
 72:   PetscCall(DMCreateGlobalVector(dmSol, &sol));
 73:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 74:   PetscCall(KSPSetOperators(ksp, A, A));
 75:   PetscCall(KSPGetPC(ksp, &pc));
 76:   PetscCall(PCSetType(pc, PCNONE));
 77:   PetscCall(KSPSetFromOptions(ksp));
 78:   PetscCall(KSPSolve(ksp, rhs, sol));

 80:   /* Check Solution */
 81:   {
 82:     Vec       diff;
 83:     PetscReal normsolRef, errAbs, errRel;

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

 98:   PetscCall(KSPDestroy(&ksp));
 99:   PetscCall(VecDestroy(&sol));
100:   PetscCall(VecDestroy(&solRef));
101:   PetscCall(VecDestroy(&rhs));
102:   PetscCall(MatDestroy(&A));
103:   PetscCall(DMDestroy(&dmSol));
104:   PetscCall(PetscFinalize());
105:   return 0;
106: }

108: PetscErrorCode ApplyOperator(Mat A, Vec in, Vec out)
109: {
110:   DM             dm;
111:   Vec            inLocal, outLocal;
112:   PetscScalar  **arrIn;
113:   PetscScalar  **arrOut;
114:   PetscInt       start, n, nExtra, ex, idxP, idxU, startGhost, nGhost;
115:   DMBoundaryType boundaryType;
116:   PetscBool      isFirst, isLast;

118:   PetscFunctionBeginUser;
119:   PetscCall(MatGetDM(A, &dm));
120:   PetscCall(DMStagGetBoundaryTypes(dm, &boundaryType, NULL, NULL));
121:   PetscCheck(boundaryType == DM_BOUNDARY_GHOSTED, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "Ghosted boundaries required");
122:   PetscCall(DMGetLocalVector(dm, &inLocal));
123:   PetscCall(DMGetLocalVector(dm, &outLocal));
124:   PetscCall(DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal));
125:   PetscCall(DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal));
126:   PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &nExtra, NULL, NULL));
127:   PetscCall(DMStagGetGhostCorners(dm, &startGhost, NULL, NULL, &nGhost, NULL, NULL));
128:   PetscCall(DMStagVecGetArrayRead(dm, inLocal, &arrIn));
129:   PetscCall(DMStagVecGetArray(dm, outLocal, &arrOut));
130:   PetscCall(DMStagGetLocationSlot(dm, LEFT, 0, &idxU));
131:   PetscCall(DMStagGetLocationSlot(dm, ELEMENT, 0, &idxP));
132:   PetscCall(DMStagGetIsFirstRank(dm, &isFirst, NULL, NULL));
133:   PetscCall(DMStagGetIsLastRank(dm, &isLast, NULL, NULL));

135:   /* Set "pressures" on ghost boundaries by copying neighboring values*/
136:   if (isFirst) arrIn[-1][idxP] = arrIn[0][idxP];
137:   if (isLast) arrIn[start + n][idxP] = arrIn[start + n - 1][idxP];

139:   /* Apply operator on physical points */
140:   for (ex = start; ex < start + n + nExtra; ++ex) {
141:     if (ex < start + n) { /* Don't compute pressure outside domain */
142:       arrOut[ex][idxP] = arrIn[ex][idxP];
143:     }
144:     arrOut[ex][idxU] = arrIn[ex][idxP] + arrIn[ex - 1][idxP] - arrIn[ex][idxU];
145:   }
146:   PetscCall(DMStagVecRestoreArrayRead(dm, inLocal, &arrIn));
147:   PetscCall(DMStagVecRestoreArray(dm, outLocal, &arrOut));
148:   PetscCall(DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out));
149:   PetscCall(DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out));
150:   PetscCall(DMRestoreLocalVector(dm, &inLocal));
151:   PetscCall(DMRestoreLocalVector(dm, &outLocal));
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: /*TEST

157:    test:
158:       suffix: 1
159:       nsize: 1

161:    test:
162:       suffix: 2
163:       nsize: 2

165:    test:
166:       suffix: 3
167:       nsize: 3
168:       args: -stag_grid_x 19

170:    test:
171:       suffix: 4
172:       nsize: 5
173:       args: -stag_grid_x 21 -stag_stencil_width 2

175: TEST*/