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;

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

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

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

 44:   /* Get the correct entries for each of our variables in local element-wise storage */
 45:   DMStagGetLocationSlot(dmSol, LEFT, 0, &iu);
 46:   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:   DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS);
 58:   DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs);
 59:   DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs);
 60:   DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol);
 61:   DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef);
 62:   DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef);
 63:   DMRestoreLocalVector(dmSol, &solRefLocal);
 64:   DMRestoreLocalVector(dmSol, &rhsLocal);

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

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

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

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

 98:   KSPDestroy(&ksp);
 99:   VecDestroy(&sol);
100:   VecDestroy(&solRef);
101:   VecDestroy(&rhs);
102:   MatDestroy(&A);
103:   DMDestroy(&dmSol);
104:   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;

119:   MatGetDM(A, &dm);
120:   DMStagGetBoundaryTypes(dm, &boundaryType, NULL, NULL);
122:   DMGetLocalVector(dm, &inLocal);
123:   DMGetLocalVector(dm, &outLocal);
124:   DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal);
125:   DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal);
126:   DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &nExtra, NULL, NULL);
127:   DMStagGetGhostCorners(dm, &startGhost, NULL, NULL, &nGhost, NULL, NULL);
128:   DMStagVecGetArrayRead(dm, inLocal, &arrIn);
129:   DMStagVecGetArray(dm, outLocal, &arrOut);
130:   DMStagGetLocationSlot(dm, LEFT, 0, &idxU);
131:   DMStagGetLocationSlot(dm, ELEMENT, 0, &idxP);
132:   DMStagGetIsFirstRank(dm, &isFirst, NULL, NULL);
133:   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:   DMStagVecRestoreArrayRead(dm, inLocal, &arrIn);
147:   DMStagVecRestoreArray(dm, outLocal, &arrOut);
148:   DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out);
149:   DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out);
150:   DMRestoreLocalVector(dm, &inLocal);
151:   DMRestoreLocalVector(dm, &outLocal);
152:   return 0;
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*/