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>
8
9 #define LEFT DMSTAG_LEFT
10 #define RIGHT DMSTAG_RIGHT
11 #define ELEMENT DMSTAG_ELEMENT
12
13 #define PRESSURE_CONST 2.0
14
15 PetscErrorCode ApplyOperator(Mat, Vec, Vec);
16
main(int argc,char ** argv)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;
27
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));
33
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));
40
41 PetscCall(DMStagGetCorners(dmSol, &start, NULL, NULL, &n, NULL, NULL, &nExtra, NULL, NULL));
42 PetscCall(DMStagVecGetArray(dmSol, rhsLocal, &arrRHS));
43
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));
65
66 /* Matrix-free Operator */
67 PetscCall(DMSetMatType(dmSol, MATSHELL));
68 PetscCall(DMCreateMatrix(dmSol, &A));
69 PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)ApplyOperator));
70
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));
79
80 /* Check Solution */
81 {
82 Vec diff;
83 PetscReal normsolRef, errAbs, errRel;
84
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 }
97
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 }
107
ApplyOperator(Mat A,Vec in,Vec out)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;
117
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));
134
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];
138
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 }
154
155 /*TEST
156
157 test:
158 suffix: 1
159 nsize: 1
160 output_file: output/empty.out
161
162 test:
163 suffix: 2
164 nsize: 2
165 output_file: output/empty.out
166
167 test:
168 suffix: 3
169 nsize: 3
170 args: -stag_grid_x 19
171 output_file: output/empty.out
172
173 test:
174 suffix: 4
175 nsize: 5
176 args: -stag_grid_x 21 -stag_stencil_width 2
177 output_file: output/empty.out
178
179 TEST*/
180