Actual source code: ex8.c
1: static char help[] = "Test DMStag ghosted boundaries in 3d\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, startz, nx, ny, nz, ex, ey, ez, nExtrax, nExtray, nExtraz;
21: PetscInt iux, iuy, iuz, ip;
22: PetscInt dof0, dof1, dof2, dof3;
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 = 0;
31: dof2 = 1;
32: dof3 = 1;
33: DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, 3, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, 1, NULL, 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, &startz, &nx, &ny, &nz, &nExtrax, &nExtray, &nExtraz);
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_BACK, 0, &iuz);
51: DMStagGetLocationSlot(dmSol, DMSTAG_ELEMENT, 0, &ip);
52: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
53: for (ey = starty; ey < starty + ny + nExtray; ++ey) {
54: for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
55: arrSol[ez][ey][ex][iux] = 2 * PRESSURE_CONST;
56: arrRHS[ez][ey][ex][iux] = 0.0;
57: arrSol[ez][ey][ex][iuy] = 2 * PRESSURE_CONST;
58: arrRHS[ez][ey][ex][iuy] = 0.0;
59: arrSol[ez][ey][ex][iuz] = 2 * PRESSURE_CONST;
60: arrRHS[ez][ey][ex][iuz] = 0.0;
61: if (ex < startx + nx && ey < starty + ny && ez < startz + nz) {
62: arrSol[ez][ey][ex][ip] = PRESSURE_CONST;
63: arrRHS[ez][ey][ex][ip] = PRESSURE_CONST;
64: }
65: }
66: }
67: }
68: DMStagVecRestoreArray(dmSol, rhsLocal, &arrRHS);
69: DMLocalToGlobalBegin(dmSol, rhsLocal, INSERT_VALUES, rhs);
70: DMLocalToGlobalEnd(dmSol, rhsLocal, INSERT_VALUES, rhs);
71: DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol);
72: DMLocalToGlobalBegin(dmSol, solRefLocal, INSERT_VALUES, solRef);
73: DMLocalToGlobalEnd(dmSol, solRefLocal, INSERT_VALUES, solRef);
74: DMRestoreLocalVector(dmSol, &solRefLocal);
75: DMRestoreLocalVector(dmSol, &rhsLocal);
77: /* Matrix-free Operator */
78: DMSetMatType(dmSol, MATSHELL);
79: DMCreateMatrix(dmSol, &A);
80: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))ApplyOperator);
82: /* Solve */
83: DMCreateGlobalVector(dmSol, &sol);
84: KSPCreate(PETSC_COMM_WORLD, &ksp);
85: KSPSetOperators(ksp, A, A);
86: KSPGetPC(ksp, &pc);
87: PCSetType(pc, PCNONE);
88: KSPSetFromOptions(ksp);
89: KSPSolve(ksp, rhs, sol);
91: /* Check Solution */
92: {
93: Vec diff;
94: PetscReal normsolRef, errAbs, errRel;
96: VecDuplicate(sol, &diff);
97: VecCopy(sol, diff);
98: VecAXPY(diff, -1.0, solRef);
99: VecNorm(diff, NORM_2, &errAbs);
100: VecNorm(solRef, NORM_2, &normsolRef);
101: errRel = errAbs / normsolRef;
102: if (errAbs > 1e14 || errRel > 1e14) {
103: PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel);
104: PetscPrintf(PetscObjectComm((PetscObject)dmSol), "Non-zero error. Probable failure.\n");
105: }
106: VecDestroy(&diff);
107: }
109: KSPDestroy(&ksp);
110: VecDestroy(&sol);
111: VecDestroy(&solRef);
112: VecDestroy(&rhs);
113: MatDestroy(&A);
114: DMDestroy(&dmSol);
115: PetscFinalize();
116: return 0;
117: }
119: PetscErrorCode ApplyOperator(Mat A, Vec in, Vec out)
120: {
121: DM dm;
122: Vec inLocal, outLocal;
123: PetscScalar ****arrIn;
124: PetscScalar ****arrOut;
125: PetscInt startx, starty, startz, nx, ny, nz, nExtrax, nExtray, nExtraz, ex, ey, ez, idxP, idxUx, idxUy, idxUz, startGhostx, startGhosty, startGhostz, nGhostx, nGhosty, nGhostz;
126: PetscBool isFirstx, isFirsty, isFirstz, isLastx, isLasty, isLastz;
129: MatGetDM(A, &dm);
130: DMGetLocalVector(dm, &inLocal);
131: DMGetLocalVector(dm, &outLocal);
132: DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal);
133: DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal);
134: DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, &nExtrax, &nExtray, &nExtraz);
135: DMStagGetGhostCorners(dm, &startGhostx, &startGhosty, &startGhostz, &nGhostx, &nGhosty, &nGhostz);
136: DMStagVecGetArrayRead(dm, inLocal, &arrIn);
137: DMStagVecGetArray(dm, outLocal, &arrOut);
138: DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &idxUx);
139: DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &idxUy);
140: DMStagGetLocationSlot(dm, DMSTAG_BACK, 0, &idxUz);
141: DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &idxP);
142: DMStagGetIsFirstRank(dm, &isFirstx, &isFirsty, &isFirstz);
143: DMStagGetIsLastRank(dm, &isLastx, &isLasty, &isLastz);
145: /* Set "pressures" on ghost boundaries by copying neighboring values*/
146: if (isFirstx) {
147: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
148: for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ez][ey][-1][idxP] = arrIn[ez][ey][0][idxP];
149: }
150: }
151: if (isLastx) {
152: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
153: for (ey = starty; ey < starty + ny + nExtray; ++ey) arrIn[ez][ey][startx + nx][idxP] = arrIn[ez][ey][startx + nx - 1][idxP];
154: }
155: }
156: if (isFirsty) {
157: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
158: for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[ez][-1][ex][idxP] = arrIn[ez][0][ex][idxP];
159: }
160: }
161: if (isLasty) {
162: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
163: for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[ez][starty + ny][ex][idxP] = arrIn[ez][starty + ny - 1][ex][idxP];
164: }
165: }
167: if (isFirstz) {
168: for (ey = starty; ey < starty + ny + nExtray; ++ey) {
169: for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[-1][ey][ex][idxP] = arrIn[0][ey][ex][idxP];
170: }
171: }
172: if (isLastz) {
173: for (ey = starty; ey < starty + ny + nExtray; ++ey) {
174: for (ex = startx; ex < startx + nx + nExtrax; ++ex) arrIn[startz + nz][ey][ex][idxP] = arrIn[startz + nz - 1][ey][ex][idxP];
175: }
176: }
178: /* Apply operator on physical points */
179: for (ez = startz; ez < startz + nz + nExtraz; ++ez) {
180: for (ey = starty; ey < starty + ny + nExtray; ++ey) {
181: for (ex = startx; ex < startx + nx + nExtrax; ++ex) {
182: if (ex < startx + nx && ey < starty + ny && ez < startz + nz) { /* Don't compute pressure outside domain */
183: arrOut[ez][ey][ex][idxP] = arrIn[ez][ey][ex][idxP];
184: }
185: arrOut[ez][ey][ex][idxUx] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey][ex - 1][idxP] - arrIn[ez][ey][ex][idxUx];
186: arrOut[ez][ey][ex][idxUy] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey - 1][ex][idxP] - arrIn[ez][ey][ex][idxUy];
187: arrOut[ez][ey][ex][idxUz] = arrIn[ez][ey][ex][idxP] + arrIn[ez - 1][ey][ex][idxP] - arrIn[ez][ey][ex][idxUz];
188: }
189: }
190: }
191: DMStagVecRestoreArrayRead(dm, inLocal, &arrIn);
192: DMStagVecRestoreArray(dm, outLocal, &arrOut);
193: DMLocalToGlobalBegin(dm, outLocal, INSERT_VALUES, out);
194: DMLocalToGlobalEnd(dm, outLocal, INSERT_VALUES, out);
195: DMRestoreLocalVector(dm, &inLocal);
196: DMRestoreLocalVector(dm, &outLocal);
197: return 0;
198: }
200: /*TEST
202: test:
203: suffix: 1
204: nsize: 1
206: test:
207: suffix: 2
208: nsize: 8
210: test:
211: suffix: 3
212: nsize: 1
213: args: -stag_stencil_width 2
215: test:
216: suffix: 4
217: nsize: 8
218: args: -stag_grid_x 4 -stag_grid_y 5 -stag_grid_z 4 -stag_stencil_width 2
220: test:
221: suffix: 5
222: nsize: 8
223: args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6 -stag_grid_z 6
225: TEST*/