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;
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 = 1;
31: dof2 = 1;
32: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof2", &dof2, NULL));
33: PetscCall(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: 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, NULL, &nx, &ny, NULL, &nExtrax, &nExtray, NULL));
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_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: PetscCall(DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS));
64: PetscCall(DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs));
65: PetscCall(DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs));
66: PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
67: PetscCall(DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef));
68: PetscCall(DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef));
69: PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
70: PetscCall(DMRestoreLocalVector(dmSol, &rhsLocal));
72: /* Matrix-free Operator */
73: PetscCall(DMSetMatType(dmSol, MATSHELL));
74: PetscCall(DMCreateMatrix(dmSol, &A));
75: PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))ApplyOperator));
77: /* Solve */
78: PetscCall(DMCreateGlobalVector(dmSol, &sol));
79: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
80: PetscCall(KSPSetOperators(ksp, A, A));
81: PetscCall(KSPGetPC(ksp, &pc));
82: PetscCall(PCSetType(pc, PCNONE));
83: PetscCall(KSPSetFromOptions(ksp));
84: PetscCall(KSPSolve(ksp, rhs, sol));
86: /* Check Solution */
87: {
88: Vec diff;
89: PetscReal normsolRef, errAbs, errRel;
91: PetscCall(VecDuplicate(sol, &diff));
92: PetscCall(VecCopy(sol, diff));
93: PetscCall(VecAXPY(diff, -1.0, solRef));
94: PetscCall(VecNorm(diff, NORM_2, &errAbs));
95: PetscCall(VecNorm(solRef, NORM_2, &normsolRef));
96: errRel = errAbs / normsolRef;
97: if (errAbs > 1e14 || errRel > 1e14) {
98: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel));
99: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Non-zero error. Probable failure.\n"));
100: }
101: PetscCall(VecDestroy(&diff));
102: }
104: PetscCall(KSPDestroy(&ksp));
105: PetscCall(VecDestroy(&sol));
106: PetscCall(VecDestroy(&solRef));
107: PetscCall(VecDestroy(&rhs));
108: PetscCall(MatDestroy(&A));
109: PetscCall(DMDestroy(&dmSol));
110: PetscCall(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;
123: PetscFunctionBeginUser;
124: PetscCall(MatGetDM(A, &dm));
125: PetscCall(DMGetLocalVector(dm, &inLocal));
126: PetscCall(DMGetLocalVector(dm, &outLocal));
127: PetscCall(DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal));
128: PetscCall(DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal));
129: PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, &nExtrax, &nExtray, NULL));
130: PetscCall(DMStagGetGhostCorners(dm, &startGhostx, &startGhosty, NULL, &nGhostx, &nGhosty, NULL));
131: PetscCall(DMStagVecGetArrayRead(dm, inLocal, &arrIn));
132: PetscCall(DMStagVecGetArray(dm, outLocal, &arrOut));
133: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &idxUx));
134: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &idxUy));
135: PetscCall(DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &idxP));
136: PetscCall(DMStagGetIsFirstRank(dm, &isFirstx, &isFirsty, &isFirstz));
137: PetscCall(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: PetscCall(DMStagVecRestoreArrayRead(dm, inLocal, &arrIn));
164: PetscCall(DMStagVecRestoreArray(dm, outLocal, &arrOut));
165: PetscCall(DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out));
166: PetscCall(DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out));
167: PetscCall(DMRestoreLocalVector(dm, &inLocal));
168: PetscCall(DMRestoreLocalVector(dm, &outLocal));
169: PetscFunctionReturn(PETSC_SUCCESS);
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*/