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, (char *)0, 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*/