Actual source code: ex18.c

  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: /*

  5:   Solves only the velocity block of a isoviscous incompressible Stokes problem on a rectangular 2D
  6:   domain.

  8:   u_xx + u_yy - p_x = f^x
  9:   v_xx + v_yy - p_y = f^y
 10:   u_x + v_y         = g

 12:   g is 0 in the physical case.

 14:   Boundary conditions give prescribed flow perpendicular to the boundaries,
 15:   and zero derivative perpendicular to them (free slip).

 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).

 22:   This is for testing purposes, and uses some routines to make
 23:   sure that transfer operators are consistent with extracting submatrices.

 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.

 29: */
 30: #include <petscdm.h>
 31: #include <petscksp.h>
 32: #include <petscdmstag.h>

 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

 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;

 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. */
 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  */
 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  */
 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  */
 69: static PetscScalar fy(PetscScalar x, PetscScalar y)
 70: {
 71:   return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y;
 72: }
 73: static PetscScalar g(PetscScalar x, PetscScalar y)
 74: {
 75:   return 0.0 * x * y;
 76: } /* identically zero */

 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;

 90:   DMStagAnalysisKSPMonitorContext mctx;
 91:   PetscBool                       analyze;

 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));

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:   }

112:   /* Create coarse DMStag (which we may or may not directly use) */
113:   PetscCall(DMCoarsen(dmSol, MPI_COMM_NULL, &dmSolc));
114:   PetscCall(DMSetUp(dmSolc));

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));

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));

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 */

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));

146:   /* Assemble velocity-velocity system */
147:   if (extractSystem) {
148:     Mat A, Ac;
149:     Vec tmp, rhs;

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"));

170:   /* Create Transfer Operators and scaling for the velocity-velocity block */
171:   if (extractTransferOperators) {
172:     Mat II, R;
173:     Vec s, tmp;

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:   }

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));

203:   PetscCall(PCMGGetCoarseSolve(pc, &kspc));
204:   PetscCall(KSPSetOperators(kspc, Auuc, Auuc));
205:   PetscCall(KSPSetFromOptions(ksp));

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:   }

215:   /* Solve */
216:   PetscCall(KSPSolve(ksp, fu, xu));

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: }

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: */
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;

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);

266:   bogusScale = 1.0;
267:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-bogus", &bogusScale, NULL)); /* Use this to break MG (try small values) */

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));

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;

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:       }

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;

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:       }

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;

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:       }

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;

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: }

672: /* A custom monitor function for analysis purposes. Computes and dumps
673:    residuals and errors for each KSP iteration */
674: PetscErrorCode DMStagAnalysisKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
675: {
676:   DM                               dm;
677:   Vec                              r, sol;
678:   DMStagAnalysisKSPMonitorContext *ctx = (DMStagAnalysisKSPMonitorContext *)mctx;

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, "Analaysis 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"));

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];

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:   }

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];

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:   }

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));

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];

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:     }

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];

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:     }

764:     PetscCall(VecDestroy(&e));
765:   }

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));

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];

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:     }

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];

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:     }

808:     PetscCall(VecDestroy(&e));
809:   }

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];

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:     }

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];

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));

856:   PetscCall(VecDestroy(&r));
857:   PetscFunctionReturn(PETSC_SUCCESS);
858: }

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) */
862: static PetscErrorCode CreateNumericalReferenceSolution(Mat A, Vec rhs, Vec *px)
863: {
864:   KSP ksp;
865:   PC  pc;
866:   Vec x;

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: }

884: /*TEST

886:    test:
887:       suffix: gmg_1
888:       nsize: 1
889:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason

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

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

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

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"

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

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

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

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

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

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

945: TEST*/