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 */
 95:   PetscInitialize(&argc, &argv, (char *)0, help);
 96:   analyze = PETSC_FALSE;
 97:   PetscOptionsGetBool(NULL, NULL, "-analyze", &analyze, NULL);
 98:   extractTransferOperators = PETSC_TRUE;
 99:   PetscOptionsGetBool(NULL, NULL, "-extractTransferOperators", &extractTransferOperators, NULL);
100:   extractSystem = PETSC_TRUE;
101:   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:     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:     DMSetFromOptions(dmSol);
109:     DMSetUp(dmSol);
110:   }

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

116:   /* Create compatible DMStags with only velocity dof (which we may or may not
117:      directly use) */
118:   DMStagCreateCompatibleDMStag(dmSol, 0, 1, 0, 0, &dmuu); /* vel-only */
119:   DMSetUp(dmuu);
120:   DMCoarsen(dmuu, MPI_COMM_NULL, &dmuuc);
121:   DMSetUp(dmuuc);

123:   /* Define uniform coordinates as a product of 1D arrays */
124:   DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
125:   DMStagSetUniformCoordinatesProduct(dmSolc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
126:   DMStagSetUniformCoordinatesProduct(dmuu, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
127:   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:   DMStagCreateISFromStencils(dmSol, 2, stencil_set, &isuf);
142:   DMStagCreateISFromStencils(dmSol, 1, &stencil_set[2], &ispf);
143:   DMStagCreateISFromStencils(dmSolc, 2, stencil_set, &isuc);
144:   DMStagCreateISFromStencils(dmSolc, 1, &stencil_set[2], &ispc);

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

151:     CreateSystem(dmSol, &A, &rhs);
152:     CreateSystem(dmSolc, &Ac, NULL);
153:     MatCreateSubMatrix(A, isuf, isuf, MAT_INITIAL_MATRIX, &Auu);
154:     MatCreateSubMatrix(Ac, isuc, isuc, MAT_INITIAL_MATRIX, &Auuc);
155:     MatCreateVecs(Auu, &xu, &fu);
156:     VecGetSubVector(rhs, isuf, &tmp);
157:     VecCopy(tmp, fu);
158:     VecRestoreSubVector(rhs, isuf, &tmp);
159:     MatDestroy(&Ac);
160:     VecDestroy(&rhs);
161:     MatDestroy(&A);
162:   } else {
163:     CreateSystem(dmuu, &Auu, &fu);
164:     CreateSystem(dmuuc, &Auuc, NULL);
165:     MatCreateVecs(Auu, &xu, NULL);
166:   }
167:   PetscObjectSetName((PetscObject)Auu, "Auu");
168:   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:     DMCreateInterpolation(dmSolc, dmSol, &II, NULL);
176:     DMCreateRestriction(dmSolc, dmSol, &R);
177:     MatCreateSubMatrix(II, isuf, isuc, MAT_INITIAL_MATRIX, &IIu);
178:     MatCreateSubMatrix(R, isuc, isuf, MAT_INITIAL_MATRIX, &Ru);
179:     DMCreateInterpolationScale(dmSolc, dmSol, II, &s);
180:     MatCreateVecs(IIu, &su, NULL);
181:     VecGetSubVector(s, isuc, &tmp);
182:     VecCopy(tmp, su);
183:     VecRestoreSubVector(s, isuc, &tmp);
184:     MatDestroy(&R);
185:     MatDestroy(&II);
186:     VecDestroy(&s);
187:   } else {
188:     DMCreateInterpolation(dmuuc, dmuu, &IIu, NULL);
189:     DMCreateInterpolationScale(dmuuc, dmuu, IIu, &su);
190:     DMCreateRestriction(dmuuc, dmuu, &Ru);
191:   }

193:   /* Create and configure solver */
194:   KSPCreate(PETSC_COMM_WORLD, &ksp);
195:   KSPSetOperators(ksp, Auu, Auu);
196:   KSPGetPC(ksp, &pc);
197:   PCSetType(pc, PCMG);
198:   PCMGSetLevels(pc, 2, NULL);
199:   PCMGSetInterpolation(pc, 1, IIu);
200:   /* PCMGSetRestriction(pc,1,Ru); */
201:   PCMGSetRScale(pc, 1, su);

203:   PCMGGetCoarseSolve(pc, &kspc);
204:   KSPSetOperators(kspc, Auuc, Auuc);
205:   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:     CreateNumericalReferenceSolution(Auu, fu, &mctx.solRefNum);
212:     KSPMonitorSet(ksp, DMStagAnalysisKSPMonitor, &mctx, NULL);
213:   }

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

218:   /* Clean up and finalize PETSc */
219:   if (analyze) {
220:     VecDestroy(&mctx.solPrev);
221:     VecDestroy(&mctx.solRef);
222:     VecDestroy(&mctx.solRefNum);
223:   }
224:   DMDestroy(&dmuu);
225:   DMDestroy(&dmuuc);
226:   VecDestroy(&su);
227:   ISDestroy(&ispc);
228:   ISDestroy(&ispf);
229:   ISDestroy(&isuc);
230:   ISDestroy(&isuf);
231:   MatDestroy(&Ru);
232:   MatDestroy(&IIu);
233:   MatDestroy(&Auuc);
234:   MatDestroy(&Auu);
235:   VecDestroy(&xu);
236:   VecDestroy(&fu);
237:   KSPDestroy(&ksp);
238:   DMDestroy(&dmSolc);
239:   DMDestroy(&dmSol);
240:   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;


262:   /* Determine whether or not to create system including pressure dof (on elements) */
263:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);
265:   hasPressure = (PetscBool)(dof[2] == 1);

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

270:   DMCreateMatrix(dm, pA);
271:   A = *pA;
272:   if (pRhs) {
273:     DMCreateGlobalVector(dm, pRhs);
274:     rhs = *pRhs;
275:   } else {
276:     rhs = NULL;
277:   }
278:   DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL);
279:   DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
280:   DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL);
281:   DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL);
282:   hx = 1.0 / N[0];
283:   hy = 1.0 / N[1];
284:   dv = hx * hy;
285:   DMStagGetProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL);
286:   DMStagGetProductCoordinateLocationSlot(dm, LEFT, &iprev);
287:   DMStagGetProductCoordinateLocationSlot(dm, RIGHT, &inext);
288:   DMStagGetProductCoordinateLocationSlot(dm, ELEMENT, &icenter);

290:   /* Loop over all local elements. Note that it may be more efficient in real
291:      applications to loop over each boundary separately */
292:   for (ey = starty; ey < starty + ny; ++ey) {
293:     for (ex = startx; ex < startx + nx; ++ex) {
294:       if (ex == N[0] - 1) {
295:         /* Right Boundary velocity Dirichlet */
296:         DMStagStencil row;
297:         PetscScalar   valRhs;

299:         const PetscScalar valA = bogusScale * 1.0;
300:         row.i                  = ex;
301:         row.j                  = ey;
302:         row.loc                = RIGHT;
303:         row.c                  = 0;
304:         DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
305:         if (rhs) {
306:           valRhs = bogusScale * uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
307:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
308:         }
309:       }
310:       if (ey == N[1] - 1) {
311:         /* Top boundary velocity Dirichlet */
312:         DMStagStencil     row;
313:         PetscScalar       valRhs;
314:         const PetscScalar valA = 1.0;
315:         row.i                  = ex;
316:         row.j                  = ey;
317:         row.loc                = UP;
318:         row.c                  = 0;
319:         DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
320:         if (rhs) {
321:           valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
322:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
323:         }
324:       }

326:       if (ey == 0) {
327:         /* Bottom boundary velocity Dirichlet */
328:         DMStagStencil     row;
329:         PetscScalar       valRhs;
330:         const PetscScalar valA = 1.0;
331:         row.i                  = ex;
332:         row.j                  = ey;
333:         row.loc                = DOWN;
334:         row.c                  = 0;
335:         DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
336:         if (rhs) {
337:           valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
338:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
339:         }
340:       } else {
341:         /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
342:         DMStagStencil row, col[7];
343:         PetscScalar   valA[7], valRhs;
344:         PetscInt      nEntries;

346:         row.i   = ex;
347:         row.j   = ey;
348:         row.loc = DOWN;
349:         row.c   = 0;
350:         if (ex == 0) {
351:           nEntries   = 4;
352:           col[0].i   = ex;
353:           col[0].j   = ey;
354:           col[0].loc = DOWN;
355:           col[0].c   = 0;
356:           valA[0]    = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
357:           col[1].i   = ex;
358:           col[1].j   = ey - 1;
359:           col[1].loc = DOWN;
360:           col[1].c   = 0;
361:           valA[1]    = dv * 1.0 / (hy * hy);
362:           col[2].i   = ex;
363:           col[2].j   = ey + 1;
364:           col[2].loc = DOWN;
365:           col[2].c   = 0;
366:           valA[2]    = dv * 1.0 / (hy * hy);
367:           /* Missing left element */
368:           col[3].i   = ex + 1;
369:           col[3].j   = ey;
370:           col[3].loc = DOWN;
371:           col[3].c   = 0;
372:           valA[3]    = dv * 1.0 / (hx * hx);
373:           if (hasPressure) {
374:             nEntries += 2;
375:             col[4].i   = ex;
376:             col[4].j   = ey - 1;
377:             col[4].loc = ELEMENT;
378:             col[4].c   = 0;
379:             valA[4]    = dv * 1.0 / hy;
380:             col[5].i   = ex;
381:             col[5].j   = ey;
382:             col[5].loc = ELEMENT;
383:             col[5].c   = 0;
384:             valA[5]    = -dv * 1.0 / hy;
385:           }
386:         } else if (ex == N[0] - 1) {
387:           /* Right boundary y velocity stencil */
388:           nEntries   = 4;
389:           col[0].i   = ex;
390:           col[0].j   = ey;
391:           col[0].loc = DOWN;
392:           col[0].c   = 0;
393:           valA[0]    = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
394:           col[1].i   = ex;
395:           col[1].j   = ey - 1;
396:           col[1].loc = DOWN;
397:           col[1].c   = 0;
398:           valA[1]    = dv * 1.0 / (hy * hy);
399:           col[2].i   = ex;
400:           col[2].j   = ey + 1;
401:           col[2].loc = DOWN;
402:           col[2].c   = 0;
403:           valA[2]    = dv * 1.0 / (hy * hy);
404:           col[3].i   = ex - 1;
405:           col[3].j   = ey;
406:           col[3].loc = DOWN;
407:           col[3].c   = 0;
408:           valA[3]    = dv * 1.0 / (hx * hx);
409:           /* Missing right element */
410:           if (hasPressure) {
411:             nEntries += 2;
412:             col[4].i   = ex;
413:             col[4].j   = ey - 1;
414:             col[4].loc = ELEMENT;
415:             col[4].c   = 0;
416:             valA[4]    = dv * 1.0 / hy;
417:             col[5].i   = ex;
418:             col[5].j   = ey;
419:             col[5].loc = ELEMENT;
420:             col[5].c   = 0;
421:             valA[5]    = -dv * 1.0 / hy;
422:           }
423:         } else {
424:           nEntries   = 5;
425:           col[0].i   = ex;
426:           col[0].j   = ey;
427:           col[0].loc = DOWN;
428:           col[0].c   = 0;
429:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
430:           col[1].i   = ex;
431:           col[1].j   = ey - 1;
432:           col[1].loc = DOWN;
433:           col[1].c   = 0;
434:           valA[1]    = dv * 1.0 / (hy * hy);
435:           col[2].i   = ex;
436:           col[2].j   = ey + 1;
437:           col[2].loc = DOWN;
438:           col[2].c   = 0;
439:           valA[2]    = dv * 1.0 / (hy * hy);
440:           col[3].i   = ex - 1;
441:           col[3].j   = ey;
442:           col[3].loc = DOWN;
443:           col[3].c   = 0;
444:           valA[3]    = dv * 1.0 / (hx * hx);
445:           col[4].i   = ex + 1;
446:           col[4].j   = ey;
447:           col[4].loc = DOWN;
448:           col[4].c   = 0;
449:           valA[4]    = dv * 1.0 / (hx * hx);
450:           if (hasPressure) {
451:             nEntries += 2;
452:             col[5].i   = ex;
453:             col[5].j   = ey - 1;
454:             col[5].loc = ELEMENT;
455:             col[5].c   = 0;
456:             valA[5]    = dv * 1.0 / hy;
457:             col[6].i   = ex;
458:             col[6].j   = ey;
459:             col[6].loc = ELEMENT;
460:             col[6].c   = 0;
461:             valA[6]    = -dv * 1.0 / hy;
462:           }
463:         }
464:         DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
465:         if (rhs) {
466:           valRhs = dv * fy(cArrX[ex][icenter], cArrY[ey][iprev]);
467:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
468:         }
469:       }

471:       if (ex == 0) {
472:         /* Left velocity Dirichlet */
473:         DMStagStencil     row;
474:         PetscScalar       valRhs;
475:         const PetscScalar valA = 1.0;
476:         row.i                  = ex;
477:         row.j                  = ey;
478:         row.loc                = LEFT;
479:         row.c                  = 0;
480:         DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
481:         if (rhs) {
482:           valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
483:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
484:         }
485:       } else {
486:         /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
487:         DMStagStencil row, col[7];
488:         PetscScalar   valA[7], valRhs;
489:         PetscInt      nEntries;
490:         row.i   = ex;
491:         row.j   = ey;
492:         row.loc = LEFT;
493:         row.c   = 0;

495:         if (ey == 0) {
496:           nEntries   = 4;
497:           col[0].i   = ex;
498:           col[0].j   = ey;
499:           col[0].loc = LEFT;
500:           col[0].c   = 0;
501:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
502:           /* missing term from element below */
503:           col[1].i   = ex;
504:           col[1].j   = ey + 1;
505:           col[1].loc = LEFT;
506:           col[1].c   = 0;
507:           valA[1]    = dv * 1.0 / (hy * hy);
508:           col[2].i   = ex - 1;
509:           col[2].j   = ey;
510:           col[2].loc = LEFT;
511:           col[2].c   = 0;
512:           valA[2]    = dv * 1.0 / (hx * hx);
513:           col[3].i   = ex + 1;
514:           col[3].j   = ey;
515:           col[3].loc = LEFT;
516:           col[3].c   = 0;
517:           valA[3]    = dv * 1.0 / (hx * hx);
518:           if (hasPressure) {
519:             nEntries += 2;
520:             col[4].i   = ex - 1;
521:             col[4].j   = ey;
522:             col[4].loc = ELEMENT;
523:             col[4].c   = 0;
524:             valA[4]    = dv * 1.0 / hx;
525:             col[5].i   = ex;
526:             col[5].j   = ey;
527:             col[5].loc = ELEMENT;
528:             col[5].c   = 0;
529:             valA[5]    = -dv * 1.0 / hx;
530:           }
531:         } else if (ey == N[1] - 1) {
532:           /* Top boundary x velocity stencil */
533:           nEntries   = 4;
534:           row.i      = ex;
535:           row.j      = ey;
536:           row.loc    = LEFT;
537:           row.c      = 0;
538:           col[0].i   = ex;
539:           col[0].j   = ey;
540:           col[0].loc = LEFT;
541:           col[0].c   = 0;
542:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
543:           col[1].i   = ex;
544:           col[1].j   = ey - 1;
545:           col[1].loc = LEFT;
546:           col[1].c   = 0;
547:           valA[1]    = dv * 1.0 / (hy * hy);
548:           /* Missing element above term */
549:           col[2].i   = ex - 1;
550:           col[2].j   = ey;
551:           col[2].loc = LEFT;
552:           col[2].c   = 0;
553:           valA[2]    = dv * 1.0 / (hx * hx);
554:           col[3].i   = ex + 1;
555:           col[3].j   = ey;
556:           col[3].loc = LEFT;
557:           col[3].c   = 0;
558:           valA[3]    = dv * 1.0 / (hx * hx);
559:           if (hasPressure) {
560:             nEntries += 2;
561:             col[4].i   = ex - 1;
562:             col[4].j   = ey;
563:             col[4].loc = ELEMENT;
564:             col[4].c   = 0;
565:             valA[4]    = dv * 1.0 / hx;
566:             col[5].i   = ex;
567:             col[5].j   = ey;
568:             col[5].loc = ELEMENT;
569:             col[5].c   = 0;
570:             valA[5]    = -dv * 1.0 / hx;
571:           }
572:         } else {
573:           /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
574:           nEntries   = 5;
575:           col[0].i   = ex;
576:           col[0].j   = ey;
577:           col[0].loc = LEFT;
578:           col[0].c   = 0;
579:           valA[0]    = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
580:           col[1].i   = ex;
581:           col[1].j   = ey - 1;
582:           col[1].loc = LEFT;
583:           col[1].c   = 0;
584:           valA[1]    = dv * 1.0 / (hy * hy);
585:           col[2].i   = ex;
586:           col[2].j   = ey + 1;
587:           col[2].loc = LEFT;
588:           col[2].c   = 0;
589:           valA[2]    = dv * 1.0 / (hy * hy);
590:           col[3].i   = ex - 1;
591:           col[3].j   = ey;
592:           col[3].loc = LEFT;
593:           col[3].c   = 0;
594:           valA[3]    = dv * 1.0 / (hx * hx);
595:           col[4].i   = ex + 1;
596:           col[4].j   = ey;
597:           col[4].loc = LEFT;
598:           col[4].c   = 0;
599:           valA[4]    = dv * 1.0 / (hx * hx);
600:           if (hasPressure) {
601:             nEntries += 2;
602:             col[5].i   = ex - 1;
603:             col[5].j   = ey;
604:             col[5].loc = ELEMENT;
605:             col[5].c   = 0;
606:             valA[5]    = dv * 1.0 / hx;
607:             col[6].i   = ex;
608:             col[6].j   = ey;
609:             col[6].loc = ELEMENT;
610:             col[6].c   = 0;
611:             valA[6]    = -dv * 1.0 / hx;
612:           }
613:         }
614:         DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
615:         if (rhs) {
616:           valRhs = dv * fx(cArrX[ex][iprev], cArrY[ey][icenter]);
617:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
618:         }
619:       }

621:       /* P equation : u_x + v_y = g
622:          Note that this includes an explicit zero on the diagonal. This is only needed for
623:          direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
624:       if (hasPressure) {
625:         DMStagStencil row, col[5];
626:         PetscScalar   valA[5], valRhs;

628:         row.i   = ex;
629:         row.j   = ey;
630:         row.loc = ELEMENT;
631:         row.c   = 0;
632:         /* Note: the scaling by dv here may not be optimal (but this test isn't concerned with these equations) */
633:         col[0].i   = ex;
634:         col[0].j   = ey;
635:         col[0].loc = LEFT;
636:         col[0].c   = 0;
637:         valA[0]    = -dv * 1.0 / hx;
638:         col[1].i   = ex;
639:         col[1].j   = ey;
640:         col[1].loc = RIGHT;
641:         col[1].c   = 0;
642:         valA[1]    = dv * 1.0 / hx;
643:         col[2].i   = ex;
644:         col[2].j   = ey;
645:         col[2].loc = DOWN;
646:         col[2].c   = 0;
647:         valA[2]    = -dv * 1.0 / hy;
648:         col[3].i   = ex;
649:         col[3].j   = ey;
650:         col[3].loc = UP;
651:         col[3].c   = 0;
652:         valA[3]    = dv * 1.0 / hy;
653:         col[4]     = row;
654:         valA[4]    = dv * 0.0;
655:         DMStagMatSetValuesStencil(dm, A, 1, &row, 5, col, valA, INSERT_VALUES);
656:         if (rhs) {
657:           valRhs = dv * g(cArrX[ex][icenter], cArrY[ey][icenter]);
658:           DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
659:         }
660:       }
661:     }
662:   }
663:   DMStagRestoreProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL);
664:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
665:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
666:   if (rhs) {
667:     VecAssemblyBegin(rhs);
668:     VecAssemblyEnd(rhs);
669:   }

671:   return 0;
672: }

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

683:   KSPBuildSolution(ksp, NULL, &sol); /* don't destroy sol */
684:   KSPBuildResidual(ksp, NULL, NULL, &r);
685:   dm = ctx->dm; /* Would typically get this with VecGetDM(), KSPGetDM() */
687:   PetscPrintf(PetscObjectComm((PetscObject)dm), " *** DMStag Analysis KSP Monitor (it. %" PetscInt_FMT ") ***\n", it);
688:   PetscPrintf(PetscObjectComm((PetscObject)dm), " Residual Norm: %g\n", (double)rnorm);
689:   PetscPrintf(PetscObjectComm((PetscObject)dm), " Dumping files...\n");

691:   /* Note: these blocks are almost entirely duplicated */
692:   {
693:     const DMStagStencilLocation loc = DMSTAG_LEFT;
694:     const PetscInt              c   = 0;
695:     Vec                         vec;
696:     DM                          da;
697:     PetscViewer                 viewer;
698:     char                        filename[PETSC_MAX_PATH_LEN];

700:     DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec);
701:     PetscSNPrintf(filename, sizeof(filename), "res_vx_%" PetscInt_FMT ".vtr", it);
702:     PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
703:     VecView(vec, viewer);
704:     PetscViewerDestroy(&viewer);
705:     VecDestroy(&vec);
706:     DMDestroy(&da);
707:   }

709:   {
710:     const DMStagStencilLocation loc = DMSTAG_DOWN;
711:     const PetscInt              c   = 0;
712:     Vec                         vec;
713:     DM                          da;
714:     PetscViewer                 viewer;
715:     char                        filename[PETSC_MAX_PATH_LEN];

717:     DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec);
718:     PetscSNPrintf(filename, sizeof(filename), "res_vy_%" PetscInt_FMT ".vtr", it);
719:     PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
720:     VecView(vec, viewer);
721:     PetscViewerDestroy(&viewer);
722:     VecDestroy(&vec);
723:     DMDestroy(&da);
724:   }

726:   if (ctx->solRef) {
727:     Vec e;
728:     VecDuplicate(ctx->solRef, &e);
729:     VecCopy(ctx->solRef, e);
730:     VecAXPY(e, -1.0, sol);

732:     {
733:       const DMStagStencilLocation loc = DMSTAG_LEFT;
734:       const PetscInt              c   = 0;
735:       Vec                         vec;
736:       DM                          da;
737:       PetscViewer                 viewer;
738:       char                        filename[PETSC_MAX_PATH_LEN];

740:       DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
741:       PetscSNPrintf(filename, sizeof(filename), "err_vx_%" PetscInt_FMT ".vtr", it);
742:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
743:       VecView(vec, viewer);
744:       PetscViewerDestroy(&viewer);
745:       VecDestroy(&vec);
746:       DMDestroy(&da);
747:     }

749:     {
750:       const DMStagStencilLocation loc = DMSTAG_DOWN;
751:       const PetscInt              c   = 0;
752:       Vec                         vec;
753:       DM                          da;
754:       PetscViewer                 viewer;
755:       char                        filename[PETSC_MAX_PATH_LEN];

757:       DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
758:       PetscSNPrintf(filename, sizeof(filename), "err_vy_%" PetscInt_FMT ".vtr", it);
759:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
760:       VecView(vec, viewer);
761:       PetscViewerDestroy(&viewer);
762:       VecDestroy(&vec);
763:       DMDestroy(&da);
764:     }

766:     VecDestroy(&e);
767:   }

769:   /* Repeat error computations wrt an "exact" solution to the discrete equations */
770:   if (ctx->solRefNum) {
771:     Vec e;
772:     VecDuplicate(ctx->solRefNum, &e);
773:     VecCopy(ctx->solRefNum, e);
774:     VecAXPY(e, -1.0, sol);

776:     {
777:       const DMStagStencilLocation loc = DMSTAG_LEFT;
778:       const PetscInt              c   = 0;
779:       Vec                         vec;
780:       DM                          da;
781:       PetscViewer                 viewer;
782:       char                        filename[PETSC_MAX_PATH_LEN];

784:       DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
785:       PetscSNPrintf(filename, sizeof(filename), "err_num_vx_%" PetscInt_FMT ".vtr", it);
786:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
787:       VecView(vec, viewer);
788:       PetscViewerDestroy(&viewer);
789:       VecDestroy(&vec);
790:       DMDestroy(&da);
791:     }

793:     {
794:       const DMStagStencilLocation loc = DMSTAG_DOWN;
795:       const PetscInt              c   = 0;
796:       Vec                         vec;
797:       DM                          da;
798:       PetscViewer                 viewer;
799:       char                        filename[PETSC_MAX_PATH_LEN];

801:       DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
802:       PetscSNPrintf(filename, sizeof(filename), "err_num_vy_%" PetscInt_FMT ".vtr", it);
803:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
804:       VecView(vec, viewer);
805:       PetscViewerDestroy(&viewer);
806:       VecDestroy(&vec);
807:       DMDestroy(&da);
808:     }

810:     VecDestroy(&e);
811:   }

813:   /* Step */
814:   if (!ctx->solPrev) {
815:     VecDuplicate(sol, &ctx->solPrev);
816:   } else {
817:     Vec diff;
818:     VecDuplicate(sol, &diff);
819:     VecCopy(sol, diff);
820:     VecAXPY(diff, -1.0, ctx->solPrev);
821:     {
822:       const DMStagStencilLocation loc = DMSTAG_LEFT;
823:       const PetscInt              c   = 0;
824:       Vec                         vec;
825:       DM                          da;
826:       PetscViewer                 viewer;
827:       char                        filename[PETSC_MAX_PATH_LEN];

829:       DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec);
830:       PetscSNPrintf(filename, sizeof(filename), "diff_vx_%" PetscInt_FMT ".vtr", it);
831:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
832:       VecView(vec, viewer);
833:       PetscViewerDestroy(&viewer);
834:       VecDestroy(&vec);
835:       DMDestroy(&da);
836:     }

838:     {
839:       const DMStagStencilLocation loc = DMSTAG_DOWN;
840:       const PetscInt              c   = 0;
841:       Vec                         vec;
842:       DM                          da;
843:       PetscViewer                 viewer;
844:       char                        filename[PETSC_MAX_PATH_LEN];

846:       DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec);
847:       PetscSNPrintf(filename, sizeof(filename), "diff_vy_%" PetscInt_FMT ".vtr", it);
848:       PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
849:       VecView(vec, viewer);
850:       PetscViewerDestroy(&viewer);
851:       VecDestroy(&vec);
852:       DMDestroy(&da);
853:     }
854:     VecDestroy(&diff);
855:   }
856:   VecCopy(sol, ctx->solPrev);

858:   VecDestroy(&r);
859:   return 0;
860: }

862: /* Use a direct solver to create an "exact" solution to the discrete system
863:    useful for testing solvers (in that it doesn't include discretization error) */
864: static PetscErrorCode CreateNumericalReferenceSolution(Mat A, Vec rhs, Vec *px)
865: {
866:   KSP ksp;
867:   PC  pc;
868:   Vec x;

871:   KSPCreate(PetscObjectComm((PetscObject)A), &ksp);
872:   KSPSetType(ksp, KSPPREONLY);
873:   KSPGetPC(ksp, &pc);
874:   PCSetType(pc, PCLU);
875:   PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK);
876:   KSPSetOptionsPrefix(ksp, "numref_");
877:   KSPSetFromOptions(ksp);
878:   KSPSetOperators(ksp, A, A);
879:   VecDuplicate(rhs, px);
880:   x = *px;
881:   KSPSolve(ksp, rhs, x);
882:   KSPDestroy(&ksp);
883:   return 0;
884: }

886: /*TEST

888:    test:
889:       suffix: gmg_1
890:       nsize: 1
891:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason

893:    test:
894:       suffix: gmg_1_b
895:       nsize: 1
896:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false

898:    test:
899:       suffix: gmg_1_c
900:       nsize: 1
901:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false

903:    test:
904:       suffix: gmg_1_d
905:       nsize: 1
906:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false

908:    test:
909:       suffix: gmg_1_bigger
910:       requires: !complex
911:       nsize: 1
912:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32

914:    test:
915:       suffix: gmg_8
916:       requires: !single !complex
917:       nsize: 8
918:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason

920:    test:
921:       suffix: gmg_8_b
922:       nsize: 8
923:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false

925:    test:
926:       suffix: gmg_8_c
927:       nsize: 8
928:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false

930:    test:
931:       suffix: gmg_8_d
932:       nsize: 8
933:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false

935:    test:
936:       suffix: gmg_8_galerkin
937:       nsize: 8
938:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false -pc_mg_galerkin

940:    test:
941:       suffix: gmg_8_bigger
942:       requires: !complex
943:       nsize: 8
944:       args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32

946: TEST*/