1 static char help[] = "Test: Solve a toy 2D problem on a staggered grid using 2-level multigrid. \n"
2 "The solve is only applied to the u-u block.\n\n";
3 /*
4
5 Solves only the velocity block of a isoviscous incompressible Stokes problem on a rectangular 2D
6 domain.
7
8 u_xx + u_yy - p_x = f^x
9 v_xx + v_yy - p_y = f^y
10 u_x + v_y = g
11
12 g is 0 in the physical case.
13
14 Boundary conditions give prescribed flow perpendicular to the boundaries,
15 and zero derivative perpendicular to them (free slip).
16
17 Supply the -analyze flag to activate a custom KSP monitor. Note that
18 this does an additional direct solve of the velocity block of the system to have an "exact"
19 solution to the discrete system available (see KSPSetOptionsPrefix
20 call below for the prefix).
21
22 This is for testing purposes, and uses some routines to make
23 sure that transfer operators are consistent with extracting submatrices.
24
25 -extractTransferOperators (true by default) defines transfer operators for
26 the velocity-velocity system by extracting submatrices from the operators for
27 the full system.
28
29 */
30 #include <petscdm.h>
31 #include <petscksp.h>
32 #include <petscdmstag.h> /* Includes petscdmproduct.h */
33
34 /* Shorter, more convenient names for DMStagStencilLocation entries */
35 #define DOWN_LEFT DMSTAG_DOWN_LEFT
36 #define DOWN DMSTAG_DOWN
37 #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
38 #define LEFT DMSTAG_LEFT
39 #define ELEMENT DMSTAG_ELEMENT
40 #define RIGHT DMSTAG_RIGHT
41 #define UP_LEFT DMSTAG_UP_LEFT
42 #define UP DMSTAG_UP
43 #define UP_RIGHT DMSTAG_UP_RIGHT
44
45 static PetscErrorCode CreateSystem(DM, Mat *, Vec *);
46 static PetscErrorCode CreateNumericalReferenceSolution(Mat, Vec, Vec *);
47 static PetscErrorCode DMStagAnalysisKSPMonitor(KSP, PetscInt, PetscReal, void *);
48 typedef struct {
49 DM dm;
50 Vec solRef, solRefNum, solPrev;
51 } DMStagAnalysisKSPMonitorContext;
52
53 /* Manufactured solution. Chosen to be higher order than can be solved exactly,
54 and to have a zero derivative for flow parallel to the boundaries. That is,
55 d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
56 and left boundaries. */
uxRef(PetscScalar x,PetscScalar y)57 static PetscScalar uxRef(PetscScalar x, PetscScalar y)
58 {
59 return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y;
60 } /* no x-dependence */
uyRef(PetscScalar x,PetscScalar y)61 static PetscScalar uyRef(PetscScalar x, PetscScalar y)
62 {
63 return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y;
64 } /* no y-dependence */
fx(PetscScalar x,PetscScalar y)65 static PetscScalar fx(PetscScalar x, PetscScalar y)
66 {
67 return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 1.0;
68 } /* no x-dependence */
fy(PetscScalar x,PetscScalar y)69 static PetscScalar fy(PetscScalar x, PetscScalar y)
70 {
71 return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y;
72 }
g(PetscScalar x,PetscScalar y)73 static PetscScalar g(PetscScalar x, PetscScalar y)
74 {
75 return 0.0 * x * y;
76 } /* identically zero */
77
main(int argc,char ** argv)78 int main(int argc, char **argv)
79 {
80 DM dmSol, dmSolc, dmuu, dmuuc;
81 KSP ksp, kspc;
82 PC pc;
83 Mat IIu, Ru, Auu, Auuc;
84 Vec su, xu, fu;
85 IS isuf, ispf, isuc, ispc;
86 PetscInt cnt = 0;
87 DMStagStencil stencil_set[1 + 1 + 1];
88 PetscBool extractTransferOperators, extractSystem;
89
90 DMStagAnalysisKSPMonitorContext mctx;
91 PetscBool analyze;
92
93 /* Initialize PETSc and process command line arguments */
94 PetscFunctionBeginUser;
95 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96 analyze = PETSC_FALSE;
97 PetscCall(PetscOptionsGetBool(NULL, NULL, "-analyze", &analyze, NULL));
98 extractTransferOperators = PETSC_TRUE;
99 PetscCall(PetscOptionsGetBool(NULL, NULL, "-extractTransferOperators", &extractTransferOperators, NULL));
100 extractSystem = PETSC_TRUE;
101 PetscCall(PetscOptionsGetBool(NULL, NULL, "-extractSystem", &extractSystem, NULL));
102
103 /* Create 2D DMStag for the solution, and set up. */
104 {
105 const PetscInt dof0 = 0, dof1 = 1, dof2 = 1; /* 1 dof on each edge and element center */
106 const PetscInt stencilWidth = 1;
107 PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 8, 8, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, &dmSol));
108 PetscCall(DMSetFromOptions(dmSol));
109 PetscCall(DMSetUp(dmSol));
110 }
111
112 /* Create coarse DMStag (which we may or may not directly use) */
113 PetscCall(DMCoarsen(dmSol, MPI_COMM_NULL, &dmSolc));
114 PetscCall(DMSetUp(dmSolc));
115
116 /* Create compatible DMStags with only velocity dof (which we may or may not
117 directly use) */
118 PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 1, 0, 0, &dmuu)); /* vel-only */
119 PetscCall(DMSetUp(dmuu));
120 PetscCall(DMCoarsen(dmuu, MPI_COMM_NULL, &dmuuc));
121 PetscCall(DMSetUp(dmuuc));
122
123 /* Define uniform coordinates as a product of 1D arrays */
124 PetscCall(DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
125 PetscCall(DMStagSetUniformCoordinatesProduct(dmSolc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
126 PetscCall(DMStagSetUniformCoordinatesProduct(dmuu, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
127 PetscCall(DMStagSetUniformCoordinatesProduct(dmuuc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
128
129 /* Create ISs for the velocity and pressure blocks */
130 /* i,j will be ignored */
131 stencil_set[cnt].loc = DMSTAG_DOWN;
132 stencil_set[cnt].c = 0;
133 cnt++; /* u */
134 stencil_set[cnt].loc = DMSTAG_LEFT;
135 stencil_set[cnt].c = 0;
136 cnt++; /* v */
137 stencil_set[cnt].loc = DMSTAG_ELEMENT;
138 stencil_set[cnt].c = 0;
139 cnt++; /* p */
140
141 PetscCall(DMStagCreateISFromStencils(dmSol, 2, stencil_set, &isuf));
142 PetscCall(DMStagCreateISFromStencils(dmSol, 1, &stencil_set[2], &ispf));
143 PetscCall(DMStagCreateISFromStencils(dmSolc, 2, stencil_set, &isuc));
144 PetscCall(DMStagCreateISFromStencils(dmSolc, 1, &stencil_set[2], &ispc));
145
146 /* Assemble velocity-velocity system */
147 if (extractSystem) {
148 Mat A, Ac;
149 Vec tmp, rhs;
150
151 PetscCall(CreateSystem(dmSol, &A, &rhs));
152 PetscCall(CreateSystem(dmSolc, &Ac, NULL));
153 PetscCall(MatCreateSubMatrix(A, isuf, isuf, MAT_INITIAL_MATRIX, &Auu));
154 PetscCall(MatCreateSubMatrix(Ac, isuc, isuc, MAT_INITIAL_MATRIX, &Auuc));
155 PetscCall(MatCreateVecs(Auu, &xu, &fu));
156 PetscCall(VecGetSubVector(rhs, isuf, &tmp));
157 PetscCall(VecCopy(tmp, fu));
158 PetscCall(VecRestoreSubVector(rhs, isuf, &tmp));
159 PetscCall(MatDestroy(&Ac));
160 PetscCall(VecDestroy(&rhs));
161 PetscCall(MatDestroy(&A));
162 } else {
163 PetscCall(CreateSystem(dmuu, &Auu, &fu));
164 PetscCall(CreateSystem(dmuuc, &Auuc, NULL));
165 PetscCall(MatCreateVecs(Auu, &xu, NULL));
166 }
167 PetscCall(PetscObjectSetName((PetscObject)Auu, "Auu"));
168 PetscCall(PetscObjectSetName((PetscObject)Auuc, "Auuc"));
169
170 /* Create Transfer Operators and scaling for the velocity-velocity block */
171 if (extractTransferOperators) {
172 Mat II, R;
173 Vec s, tmp;
174
175 PetscCall(DMCreateInterpolation(dmSolc, dmSol, &II, NULL));
176 PetscCall(DMCreateRestriction(dmSolc, dmSol, &R));
177 PetscCall(MatCreateSubMatrix(II, isuf, isuc, MAT_INITIAL_MATRIX, &IIu));
178 PetscCall(MatCreateSubMatrix(R, isuc, isuf, MAT_INITIAL_MATRIX, &Ru));
179 PetscCall(DMCreateInterpolationScale(dmSolc, dmSol, II, &s));
180 PetscCall(MatCreateVecs(IIu, &su, NULL));
181 PetscCall(VecGetSubVector(s, isuc, &tmp));
182 PetscCall(VecCopy(tmp, su));
183 PetscCall(VecRestoreSubVector(s, isuc, &tmp));
184 PetscCall(MatDestroy(&R));
185 PetscCall(MatDestroy(&II));
186 PetscCall(VecDestroy(&s));
187 } else {
188 PetscCall(DMCreateInterpolation(dmuuc, dmuu, &IIu, NULL));
189 PetscCall(DMCreateInterpolationScale(dmuuc, dmuu, IIu, &su));
190 PetscCall(DMCreateRestriction(dmuuc, dmuu, &Ru));
191 }
192
193 /* Create and configure solver */
194 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
195 PetscCall(KSPSetOperators(ksp, Auu, Auu));
196 PetscCall(KSPGetPC(ksp, &pc));
197 PetscCall(PCSetType(pc, PCMG));
198 PetscCall(PCMGSetLevels(pc, 2, NULL));
199 PetscCall(PCMGSetInterpolation(pc, 1, IIu));
200 /* PetscCall(PCMGSetRestriction(pc,1,Ru)); */
201 PetscCall(PCMGSetRScale(pc, 1, su));
202
203 PetscCall(PCMGGetCoarseSolve(pc, &kspc));
204 PetscCall(KSPSetOperators(kspc, Auuc, Auuc));
205 PetscCall(KSPSetFromOptions(ksp));
206
207 if (analyze) {
208 mctx.dm = dmuu;
209 mctx.solRef = NULL; /* Reference solution not computed for u-u only */
210 mctx.solPrev = NULL; /* Populated automatically */
211 PetscCall(CreateNumericalReferenceSolution(Auu, fu, &mctx.solRefNum));
212 PetscCall(KSPMonitorSet(ksp, DMStagAnalysisKSPMonitor, &mctx, NULL));
213 }
214
215 /* Solve */
216 PetscCall(KSPSolve(ksp, fu, xu));
217
218 /* Clean up and finalize PETSc */
219 if (analyze) {
220 PetscCall(VecDestroy(&mctx.solPrev));
221 PetscCall(VecDestroy(&mctx.solRef));
222 PetscCall(VecDestroy(&mctx.solRefNum));
223 }
224 PetscCall(DMDestroy(&dmuu));
225 PetscCall(DMDestroy(&dmuuc));
226 PetscCall(VecDestroy(&su));
227 PetscCall(ISDestroy(&ispc));
228 PetscCall(ISDestroy(&ispf));
229 PetscCall(ISDestroy(&isuc));
230 PetscCall(ISDestroy(&isuf));
231 PetscCall(MatDestroy(&Ru));
232 PetscCall(MatDestroy(&IIu));
233 PetscCall(MatDestroy(&Auuc));
234 PetscCall(MatDestroy(&Auu));
235 PetscCall(VecDestroy(&xu));
236 PetscCall(VecDestroy(&fu));
237 PetscCall(KSPDestroy(&ksp));
238 PetscCall(DMDestroy(&dmSolc));
239 PetscCall(DMDestroy(&dmSol));
240 PetscCall(PetscFinalize());
241 return 0;
242 }
243
244 /*
245 Note: in this system all stencil coefficients which are not related to the Dirichlet boundary are scaled by dv = dx*dy.
246 This scaling is necessary for multigrid to converge.
247 */
CreateSystem(DM dm,Mat * pA,Vec * pRhs)248 static PetscErrorCode CreateSystem(DM dm, Mat *pA, Vec *pRhs)
249 {
250 PetscInt N[2], dof[3];
251 PetscBool isLastRankx, isLastRanky, isFirstRankx, isFirstRanky;
252 PetscInt ex, ey, startx, starty, nx, ny;
253 PetscInt iprev, icenter, inext;
254 Mat A;
255 Vec rhs;
256 PetscReal hx, hy, dv, bogusScale;
257 PetscScalar **cArrX, **cArrY;
258 PetscBool hasPressure;
259
260 PetscFunctionBeginUser;
261 /* Determine whether or not to create system including pressure dof (on elements) */
262 PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL));
263 PetscCheck(dof[0] == 0 && dof[1] == 1 && (dof[2] == 1 || dof[2] == 0), PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CreateSystem only implemented for velocity-only or velocity+pressure grids");
264 hasPressure = (PetscBool)(dof[2] == 1);
265
266 bogusScale = 1.0;
267 PetscCall(PetscOptionsGetReal(NULL, NULL, "-bogus", &bogusScale, NULL)); /* Use this to break MG (try small values) */
268
269 PetscCall(DMCreateMatrix(dm, pA));
270 A = *pA;
271 if (pRhs) {
272 PetscCall(DMCreateGlobalVector(dm, pRhs));
273 rhs = *pRhs;
274 } else {
275 rhs = NULL;
276 }
277 PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
278 PetscCall(DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL));
279 PetscCall(DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL));
280 PetscCall(DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL));
281 hx = 1.0 / N[0];
282 hy = 1.0 / N[1];
283 dv = hx * hy;
284 PetscCall(DMStagGetProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL));
285 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, LEFT, &iprev));
286 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, RIGHT, &inext));
287 PetscCall(DMStagGetProductCoordinateLocationSlot(dm, ELEMENT, &icenter));
288
289 /* Loop over all local elements. Note that it may be more efficient in real
290 applications to loop over each boundary separately */
291 for (ey = starty; ey < starty + ny; ++ey) {
292 for (ex = startx; ex < startx + nx; ++ex) {
293 if (ex == N[0] - 1) {
294 /* Right Boundary velocity Dirichlet */
295 DMStagStencil row;
296 PetscScalar valRhs;
297
298 const PetscScalar valA = bogusScale * 1.0;
299 row.i = ex;
300 row.j = ey;
301 row.loc = RIGHT;
302 row.c = 0;
303 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
304 if (rhs) {
305 valRhs = bogusScale * uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
306 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
307 }
308 }
309 if (ey == N[1] - 1) {
310 /* Top boundary velocity Dirichlet */
311 DMStagStencil row;
312 PetscScalar valRhs;
313 const PetscScalar valA = 1.0;
314 row.i = ex;
315 row.j = ey;
316 row.loc = UP;
317 row.c = 0;
318 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
319 if (rhs) {
320 valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
321 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
322 }
323 }
324
325 if (ey == 0) {
326 /* Bottom boundary velocity Dirichlet */
327 DMStagStencil row;
328 PetscScalar valRhs;
329 const PetscScalar valA = 1.0;
330 row.i = ex;
331 row.j = ey;
332 row.loc = DOWN;
333 row.c = 0;
334 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
335 if (rhs) {
336 valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
337 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
338 }
339 } else {
340 /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
341 DMStagStencil row, col[7];
342 PetscScalar valA[7], valRhs;
343 PetscInt nEntries;
344
345 row.i = ex;
346 row.j = ey;
347 row.loc = DOWN;
348 row.c = 0;
349 if (ex == 0) {
350 nEntries = 4;
351 col[0].i = ex;
352 col[0].j = ey;
353 col[0].loc = DOWN;
354 col[0].c = 0;
355 valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
356 col[1].i = ex;
357 col[1].j = ey - 1;
358 col[1].loc = DOWN;
359 col[1].c = 0;
360 valA[1] = dv * 1.0 / (hy * hy);
361 col[2].i = ex;
362 col[2].j = ey + 1;
363 col[2].loc = DOWN;
364 col[2].c = 0;
365 valA[2] = dv * 1.0 / (hy * hy);
366 /* Missing left element */
367 col[3].i = ex + 1;
368 col[3].j = ey;
369 col[3].loc = DOWN;
370 col[3].c = 0;
371 valA[3] = dv * 1.0 / (hx * hx);
372 if (hasPressure) {
373 nEntries += 2;
374 col[4].i = ex;
375 col[4].j = ey - 1;
376 col[4].loc = ELEMENT;
377 col[4].c = 0;
378 valA[4] = dv * 1.0 / hy;
379 col[5].i = ex;
380 col[5].j = ey;
381 col[5].loc = ELEMENT;
382 col[5].c = 0;
383 valA[5] = -dv * 1.0 / hy;
384 }
385 } else if (ex == N[0] - 1) {
386 /* Right boundary y velocity stencil */
387 nEntries = 4;
388 col[0].i = ex;
389 col[0].j = ey;
390 col[0].loc = DOWN;
391 col[0].c = 0;
392 valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
393 col[1].i = ex;
394 col[1].j = ey - 1;
395 col[1].loc = DOWN;
396 col[1].c = 0;
397 valA[1] = dv * 1.0 / (hy * hy);
398 col[2].i = ex;
399 col[2].j = ey + 1;
400 col[2].loc = DOWN;
401 col[2].c = 0;
402 valA[2] = dv * 1.0 / (hy * hy);
403 col[3].i = ex - 1;
404 col[3].j = ey;
405 col[3].loc = DOWN;
406 col[3].c = 0;
407 valA[3] = dv * 1.0 / (hx * hx);
408 /* Missing right element */
409 if (hasPressure) {
410 nEntries += 2;
411 col[4].i = ex;
412 col[4].j = ey - 1;
413 col[4].loc = ELEMENT;
414 col[4].c = 0;
415 valA[4] = dv * 1.0 / hy;
416 col[5].i = ex;
417 col[5].j = ey;
418 col[5].loc = ELEMENT;
419 col[5].c = 0;
420 valA[5] = -dv * 1.0 / hy;
421 }
422 } else {
423 nEntries = 5;
424 col[0].i = ex;
425 col[0].j = ey;
426 col[0].loc = DOWN;
427 col[0].c = 0;
428 valA[0] = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
429 col[1].i = ex;
430 col[1].j = ey - 1;
431 col[1].loc = DOWN;
432 col[1].c = 0;
433 valA[1] = dv * 1.0 / (hy * hy);
434 col[2].i = ex;
435 col[2].j = ey + 1;
436 col[2].loc = DOWN;
437 col[2].c = 0;
438 valA[2] = dv * 1.0 / (hy * hy);
439 col[3].i = ex - 1;
440 col[3].j = ey;
441 col[3].loc = DOWN;
442 col[3].c = 0;
443 valA[3] = dv * 1.0 / (hx * hx);
444 col[4].i = ex + 1;
445 col[4].j = ey;
446 col[4].loc = DOWN;
447 col[4].c = 0;
448 valA[4] = dv * 1.0 / (hx * hx);
449 if (hasPressure) {
450 nEntries += 2;
451 col[5].i = ex;
452 col[5].j = ey - 1;
453 col[5].loc = ELEMENT;
454 col[5].c = 0;
455 valA[5] = dv * 1.0 / hy;
456 col[6].i = ex;
457 col[6].j = ey;
458 col[6].loc = ELEMENT;
459 col[6].c = 0;
460 valA[6] = -dv * 1.0 / hy;
461 }
462 }
463 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
464 if (rhs) {
465 valRhs = dv * fy(cArrX[ex][icenter], cArrY[ey][iprev]);
466 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
467 }
468 }
469
470 if (ex == 0) {
471 /* Left velocity Dirichlet */
472 DMStagStencil row;
473 PetscScalar valRhs;
474 const PetscScalar valA = 1.0;
475 row.i = ex;
476 row.j = ey;
477 row.loc = LEFT;
478 row.c = 0;
479 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
480 if (rhs) {
481 valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
482 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
483 }
484 } else {
485 /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
486 DMStagStencil row, col[7];
487 PetscScalar valA[7], valRhs;
488 PetscInt nEntries;
489 row.i = ex;
490 row.j = ey;
491 row.loc = LEFT;
492 row.c = 0;
493
494 if (ey == 0) {
495 nEntries = 4;
496 col[0].i = ex;
497 col[0].j = ey;
498 col[0].loc = LEFT;
499 col[0].c = 0;
500 valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
501 /* missing term from element below */
502 col[1].i = ex;
503 col[1].j = ey + 1;
504 col[1].loc = LEFT;
505 col[1].c = 0;
506 valA[1] = dv * 1.0 / (hy * hy);
507 col[2].i = ex - 1;
508 col[2].j = ey;
509 col[2].loc = LEFT;
510 col[2].c = 0;
511 valA[2] = dv * 1.0 / (hx * hx);
512 col[3].i = ex + 1;
513 col[3].j = ey;
514 col[3].loc = LEFT;
515 col[3].c = 0;
516 valA[3] = dv * 1.0 / (hx * hx);
517 if (hasPressure) {
518 nEntries += 2;
519 col[4].i = ex - 1;
520 col[4].j = ey;
521 col[4].loc = ELEMENT;
522 col[4].c = 0;
523 valA[4] = dv * 1.0 / hx;
524 col[5].i = ex;
525 col[5].j = ey;
526 col[5].loc = ELEMENT;
527 col[5].c = 0;
528 valA[5] = -dv * 1.0 / hx;
529 }
530 } else if (ey == N[1] - 1) {
531 /* Top boundary x velocity stencil */
532 nEntries = 4;
533 row.i = ex;
534 row.j = ey;
535 row.loc = LEFT;
536 row.c = 0;
537 col[0].i = ex;
538 col[0].j = ey;
539 col[0].loc = LEFT;
540 col[0].c = 0;
541 valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
542 col[1].i = ex;
543 col[1].j = ey - 1;
544 col[1].loc = LEFT;
545 col[1].c = 0;
546 valA[1] = dv * 1.0 / (hy * hy);
547 /* Missing element above term */
548 col[2].i = ex - 1;
549 col[2].j = ey;
550 col[2].loc = LEFT;
551 col[2].c = 0;
552 valA[2] = dv * 1.0 / (hx * hx);
553 col[3].i = ex + 1;
554 col[3].j = ey;
555 col[3].loc = LEFT;
556 col[3].c = 0;
557 valA[3] = dv * 1.0 / (hx * hx);
558 if (hasPressure) {
559 nEntries += 2;
560 col[4].i = ex - 1;
561 col[4].j = ey;
562 col[4].loc = ELEMENT;
563 col[4].c = 0;
564 valA[4] = dv * 1.0 / hx;
565 col[5].i = ex;
566 col[5].j = ey;
567 col[5].loc = ELEMENT;
568 col[5].c = 0;
569 valA[5] = -dv * 1.0 / hx;
570 }
571 } else {
572 /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
573 nEntries = 5;
574 col[0].i = ex;
575 col[0].j = ey;
576 col[0].loc = LEFT;
577 col[0].c = 0;
578 valA[0] = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
579 col[1].i = ex;
580 col[1].j = ey - 1;
581 col[1].loc = LEFT;
582 col[1].c = 0;
583 valA[1] = dv * 1.0 / (hy * hy);
584 col[2].i = ex;
585 col[2].j = ey + 1;
586 col[2].loc = LEFT;
587 col[2].c = 0;
588 valA[2] = dv * 1.0 / (hy * hy);
589 col[3].i = ex - 1;
590 col[3].j = ey;
591 col[3].loc = LEFT;
592 col[3].c = 0;
593 valA[3] = dv * 1.0 / (hx * hx);
594 col[4].i = ex + 1;
595 col[4].j = ey;
596 col[4].loc = LEFT;
597 col[4].c = 0;
598 valA[4] = dv * 1.0 / (hx * hx);
599 if (hasPressure) {
600 nEntries += 2;
601 col[5].i = ex - 1;
602 col[5].j = ey;
603 col[5].loc = ELEMENT;
604 col[5].c = 0;
605 valA[5] = dv * 1.0 / hx;
606 col[6].i = ex;
607 col[6].j = ey;
608 col[6].loc = ELEMENT;
609 col[6].c = 0;
610 valA[6] = -dv * 1.0 / hx;
611 }
612 }
613 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
614 if (rhs) {
615 valRhs = dv * fx(cArrX[ex][iprev], cArrY[ey][icenter]);
616 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
617 }
618 }
619
620 /* P equation : u_x + v_y = g
621 Note that this includes an explicit zero on the diagonal. This is only needed for
622 direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
623 if (hasPressure) {
624 DMStagStencil row, col[5];
625 PetscScalar valA[5], valRhs;
626
627 row.i = ex;
628 row.j = ey;
629 row.loc = ELEMENT;
630 row.c = 0;
631 /* Note: the scaling by dv here may not be optimal (but this test isn't concerned with these equations) */
632 col[0].i = ex;
633 col[0].j = ey;
634 col[0].loc = LEFT;
635 col[0].c = 0;
636 valA[0] = -dv * 1.0 / hx;
637 col[1].i = ex;
638 col[1].j = ey;
639 col[1].loc = RIGHT;
640 col[1].c = 0;
641 valA[1] = dv * 1.0 / hx;
642 col[2].i = ex;
643 col[2].j = ey;
644 col[2].loc = DOWN;
645 col[2].c = 0;
646 valA[2] = -dv * 1.0 / hy;
647 col[3].i = ex;
648 col[3].j = ey;
649 col[3].loc = UP;
650 col[3].c = 0;
651 valA[3] = dv * 1.0 / hy;
652 col[4] = row;
653 valA[4] = dv * 0.0;
654 PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 5, col, valA, INSERT_VALUES));
655 if (rhs) {
656 valRhs = dv * g(cArrX[ex][icenter], cArrY[ey][icenter]);
657 PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
658 }
659 }
660 }
661 }
662 PetscCall(DMStagRestoreProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL));
663 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
664 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
665 if (rhs) {
666 PetscCall(VecAssemblyBegin(rhs));
667 PetscCall(VecAssemblyEnd(rhs));
668 }
669 PetscFunctionReturn(PETSC_SUCCESS);
670 }
671
672 /* A custom monitor function for analysis purposes. Computes and dumps
673 residuals and errors for each KSP iteration */
DMStagAnalysisKSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm,void * mctx)674 PetscErrorCode DMStagAnalysisKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
675 {
676 DM dm;
677 Vec r, sol;
678 DMStagAnalysisKSPMonitorContext *ctx = (DMStagAnalysisKSPMonitorContext *)mctx;
679
680 PetscFunctionBeginUser;
681 PetscCall(KSPBuildSolution(ksp, NULL, &sol)); /* don't destroy sol */
682 PetscCall(KSPBuildResidual(ksp, NULL, NULL, &r));
683 dm = ctx->dm; /* Would typically get this with VecGetDM(), KSPGetDM() */
684 PetscCheck(dm, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Analysis monitor requires a DM which is properly associated with the solution Vec");
685 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " *** DMStag Analysis KSP Monitor (it. %" PetscInt_FMT ") ***\n", it));
686 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " Residual Norm: %g\n", (double)rnorm));
687 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " Dumping files...\n"));
688
689 /* Note: these blocks are almost entirely duplicated */
690 {
691 const DMStagStencilLocation loc = DMSTAG_LEFT;
692 const PetscInt c = 0;
693 Vec vec;
694 DM da;
695 PetscViewer viewer;
696 char filename[PETSC_MAX_PATH_LEN];
697
698 PetscCall(DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec));
699 PetscCall(PetscSNPrintf(filename, sizeof(filename), "res_vx_%" PetscInt_FMT ".vtr", it));
700 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
701 PetscCall(VecView(vec, viewer));
702 PetscCall(PetscViewerDestroy(&viewer));
703 PetscCall(VecDestroy(&vec));
704 PetscCall(DMDestroy(&da));
705 }
706
707 {
708 const DMStagStencilLocation loc = DMSTAG_DOWN;
709 const PetscInt c = 0;
710 Vec vec;
711 DM da;
712 PetscViewer viewer;
713 char filename[PETSC_MAX_PATH_LEN];
714
715 PetscCall(DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec));
716 PetscCall(PetscSNPrintf(filename, sizeof(filename), "res_vy_%" PetscInt_FMT ".vtr", it));
717 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
718 PetscCall(VecView(vec, viewer));
719 PetscCall(PetscViewerDestroy(&viewer));
720 PetscCall(VecDestroy(&vec));
721 PetscCall(DMDestroy(&da));
722 }
723
724 if (ctx->solRef) {
725 Vec e;
726 PetscCall(VecDuplicate(ctx->solRef, &e));
727 PetscCall(VecCopy(ctx->solRef, e));
728 PetscCall(VecAXPY(e, -1.0, sol));
729
730 {
731 const DMStagStencilLocation loc = DMSTAG_LEFT;
732 const PetscInt c = 0;
733 Vec vec;
734 DM da;
735 PetscViewer viewer;
736 char filename[PETSC_MAX_PATH_LEN];
737
738 PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
739 PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_vx_%" PetscInt_FMT ".vtr", it));
740 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
741 PetscCall(VecView(vec, viewer));
742 PetscCall(PetscViewerDestroy(&viewer));
743 PetscCall(VecDestroy(&vec));
744 PetscCall(DMDestroy(&da));
745 }
746
747 {
748 const DMStagStencilLocation loc = DMSTAG_DOWN;
749 const PetscInt c = 0;
750 Vec vec;
751 DM da;
752 PetscViewer viewer;
753 char filename[PETSC_MAX_PATH_LEN];
754
755 PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
756 PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_vy_%" PetscInt_FMT ".vtr", it));
757 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
758 PetscCall(VecView(vec, viewer));
759 PetscCall(PetscViewerDestroy(&viewer));
760 PetscCall(VecDestroy(&vec));
761 PetscCall(DMDestroy(&da));
762 }
763
764 PetscCall(VecDestroy(&e));
765 }
766
767 /* Repeat error computations wrt an "exact" solution to the discrete equations */
768 if (ctx->solRefNum) {
769 Vec e;
770 PetscCall(VecDuplicate(ctx->solRefNum, &e));
771 PetscCall(VecCopy(ctx->solRefNum, e));
772 PetscCall(VecAXPY(e, -1.0, sol));
773
774 {
775 const DMStagStencilLocation loc = DMSTAG_LEFT;
776 const PetscInt c = 0;
777 Vec vec;
778 DM da;
779 PetscViewer viewer;
780 char filename[PETSC_MAX_PATH_LEN];
781
782 PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
783 PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_num_vx_%" PetscInt_FMT ".vtr", it));
784 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
785 PetscCall(VecView(vec, viewer));
786 PetscCall(PetscViewerDestroy(&viewer));
787 PetscCall(VecDestroy(&vec));
788 PetscCall(DMDestroy(&da));
789 }
790
791 {
792 const DMStagStencilLocation loc = DMSTAG_DOWN;
793 const PetscInt c = 0;
794 Vec vec;
795 DM da;
796 PetscViewer viewer;
797 char filename[PETSC_MAX_PATH_LEN];
798
799 PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
800 PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_num_vy_%" PetscInt_FMT ".vtr", it));
801 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
802 PetscCall(VecView(vec, viewer));
803 PetscCall(PetscViewerDestroy(&viewer));
804 PetscCall(VecDestroy(&vec));
805 PetscCall(DMDestroy(&da));
806 }
807
808 PetscCall(VecDestroy(&e));
809 }
810
811 /* Step */
812 if (!ctx->solPrev) {
813 PetscCall(VecDuplicate(sol, &ctx->solPrev));
814 } else {
815 Vec diff;
816 PetscCall(VecDuplicate(sol, &diff));
817 PetscCall(VecCopy(sol, diff));
818 PetscCall(VecAXPY(diff, -1.0, ctx->solPrev));
819 {
820 const DMStagStencilLocation loc = DMSTAG_LEFT;
821 const PetscInt c = 0;
822 Vec vec;
823 DM da;
824 PetscViewer viewer;
825 char filename[PETSC_MAX_PATH_LEN];
826
827 PetscCall(DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec));
828 PetscCall(PetscSNPrintf(filename, sizeof(filename), "diff_vx_%" PetscInt_FMT ".vtr", it));
829 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
830 PetscCall(VecView(vec, viewer));
831 PetscCall(PetscViewerDestroy(&viewer));
832 PetscCall(VecDestroy(&vec));
833 PetscCall(DMDestroy(&da));
834 }
835
836 {
837 const DMStagStencilLocation loc = DMSTAG_DOWN;
838 const PetscInt c = 0;
839 Vec vec;
840 DM da;
841 PetscViewer viewer;
842 char filename[PETSC_MAX_PATH_LEN];
843
844 PetscCall(DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec));
845 PetscCall(PetscSNPrintf(filename, sizeof(filename), "diff_vy_%" PetscInt_FMT ".vtr", it));
846 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
847 PetscCall(VecView(vec, viewer));
848 PetscCall(PetscViewerDestroy(&viewer));
849 PetscCall(VecDestroy(&vec));
850 PetscCall(DMDestroy(&da));
851 }
852 PetscCall(VecDestroy(&diff));
853 }
854 PetscCall(VecCopy(sol, ctx->solPrev));
855
856 PetscCall(VecDestroy(&r));
857 PetscFunctionReturn(PETSC_SUCCESS);
858 }
859
860 /* Use a direct solver to create an "exact" solution to the discrete system
861 useful for testing solvers (in that it doesn't include discretization error) */
CreateNumericalReferenceSolution(Mat A,Vec rhs,Vec * px)862 static PetscErrorCode CreateNumericalReferenceSolution(Mat A, Vec rhs, Vec *px)
863 {
864 KSP ksp;
865 PC pc;
866 Vec x;
867
868 PetscFunctionBeginUser;
869 PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
870 PetscCall(KSPSetType(ksp, KSPPREONLY));
871 PetscCall(KSPGetPC(ksp, &pc));
872 PetscCall(PCSetType(pc, PCLU));
873 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK));
874 PetscCall(KSPSetOptionsPrefix(ksp, "numref_"));
875 PetscCall(KSPSetFromOptions(ksp));
876 PetscCall(KSPSetOperators(ksp, A, A));
877 PetscCall(VecDuplicate(rhs, px));
878 x = *px;
879 PetscCall(KSPSolve(ksp, rhs, x));
880 PetscCall(KSPDestroy(&ksp));
881 PetscFunctionReturn(PETSC_SUCCESS);
882 }
883
884 /*TEST
885
886 test:
887 suffix: gmg_1
888 nsize: 1
889 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
890
891 test:
892 suffix: gmg_1_b
893 nsize: 1
894 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
895
896 test:
897 suffix: gmg_1_c
898 nsize: 1
899 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
900
901 test:
902 suffix: gmg_1_d
903 nsize: 1
904 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
905
906 test:
907 suffix: gmg_1_bigger
908 requires: !complex
909 nsize: 1
910 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
911 filter: sed -e "s/ iterations 8/ iterations 7/g"
912
913 test:
914 suffix: gmg_8
915 requires: !single !complex
916 nsize: 8
917 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
918
919 test:
920 suffix: gmg_8_b
921 nsize: 8
922 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
923
924 test:
925 suffix: gmg_8_c
926 nsize: 8
927 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
928
929 test:
930 suffix: gmg_8_d
931 nsize: 8
932 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
933
934 test:
935 suffix: gmg_8_galerkin
936 nsize: 8
937 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false -pc_mg_galerkin
938
939 test:
940 suffix: gmg_8_bigger
941 requires: !complex
942 nsize: 8
943 args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
944
945 TEST*/
946