Actual source code: ex2.c

  1: static char help[] = "Solve a toy 2D problem on a staggered grid\n\n";
  2: /*

  4:   To demonstrate the basic functionality of DMStag, solves an isoviscous
  5:   incompressible Stokes problem on a rectangular 2D domain, using a manufactured
  6:   solution.

  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:   Use the -pinpressure option to fix a pressure node, instead of providing
 18:   a constant-pressure nullspace. This allows use of direct solvers, e.g. to
 19:   use UMFPACK,

 21:      ./ex2 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack

 23:   This example demonstrates the use of DMProduct to efficiently store coordinates
 24:   on an orthogonal grid.

 26: */
 27: #include <petscdm.h>
 28: #include <petscksp.h>
 29: #include <petscdmstag.h>

 31: /* Shorter, more convenient names for DMStagStencilLocation entries */
 32: #define DOWN_LEFT  DMSTAG_DOWN_LEFT
 33: #define DOWN       DMSTAG_DOWN
 34: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
 35: #define LEFT       DMSTAG_LEFT
 36: #define ELEMENT    DMSTAG_ELEMENT
 37: #define RIGHT      DMSTAG_RIGHT
 38: #define UP_LEFT    DMSTAG_UP_LEFT
 39: #define UP         DMSTAG_UP
 40: #define UP_RIGHT   DMSTAG_UP_RIGHT

 42: static PetscErrorCode CreateSystem(DM, Mat *, Vec *, PetscBool);
 43: static PetscErrorCode CreateReferenceSolution(DM, Vec *);
 44: static PetscErrorCode AttachNullspace(DM, Mat);
 45: static PetscErrorCode CheckSolution(Vec, Vec);

 47: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
 48: and to have a zero derivative for flow parallel to the boundaries. That is,
 49: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
 50: and left boundaries. */
 51: static PetscScalar uxRef(PetscScalar x, PetscScalar y)
 52: {
 53:   return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y;
 54: } /* no x-dependence  */
 55: static PetscScalar uyRef(PetscScalar x, PetscScalar y)
 56: {
 57:   return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y;
 58: } /* no y-dependence  */
 59: static PetscScalar pRef(PetscScalar x, PetscScalar y)
 60: {
 61:   return -1.0 * (x - 0.5) + -3.0 / 2.0 * y * y + 0.5;
 62: } /* zero integral    */
 63: static PetscScalar fx(PetscScalar x, PetscScalar y)
 64: {
 65:   return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 1.0;
 66: } /* no x-dependence  */
 67: static PetscScalar fy(PetscScalar x, PetscScalar y)
 68: {
 69:   return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y;
 70: }
 71: static PetscScalar g(PetscScalar x, PetscScalar y)
 72: {
 73:   return 0.0 * x * y;
 74: } /* identically zero */

 76: int main(int argc, char **argv)
 77: {
 78:   DM        dmSol;
 79:   Vec       sol, solRef, rhs;
 80:   Mat       A;
 81:   KSP       ksp;
 82:   PC        pc;
 83:   PetscBool pinPressure;

 85:   /* Initialize PETSc and process command line arguments */
 86:   PetscFunctionBeginUser;
 87:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 88:   pinPressure = PETSC_TRUE;
 89:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pinpressure", &pinPressure, NULL));

 91:   /* Create 2D DMStag for the solution, and set up. */
 92:   {
 93:     const PetscInt dof0 = 0, dof1 = 1, dof2 = 1; /* 1 dof on each edge and element center */
 94:     const PetscInt stencilWidth = 1;
 95:     PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 7, 9, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, &dmSol));
 96:     PetscCall(DMSetFromOptions(dmSol));
 97:     PetscCall(DMSetUp(dmSol));
 98:   }

100:   /* Define uniform coordinates as a product of 1D arrays */
101:   PetscCall(DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));

103:   /* Compute (manufactured) reference solution */
104:   PetscCall(CreateReferenceSolution(dmSol, &solRef));

106:   /* Assemble system */
107:   PetscCall(CreateSystem(dmSol, &A, &rhs, pinPressure));

109:   /* Attach a constant-pressure nullspace to the operator
110:   (logically, this should be in CreateSystem, but we separate it here to highlight
111:    the features of DMStag exposed, in particular DMStagMigrateVec()) */
112:   if (!pinPressure) PetscCall(AttachNullspace(dmSol, A));

114:   /* Solve, using the default FieldSplit (Approximate Block Factorization) Preconditioner
115:      This is not intended to be an example of a good solver!  */
116:   PetscCall(DMCreateGlobalVector(dmSol, &sol));
117:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
118:   PetscCall(KSPSetType(ksp, KSPFGMRES));
119:   PetscCall(KSPSetOperators(ksp, A, A));
120:   PetscCall(KSPGetPC(ksp, &pc));
121:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
122:   PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE));
123:   PetscCall(KSPSetFromOptions(ksp));
124:   PetscCall(KSPSolve(ksp, rhs, sol));

126:   /* Check Solution */
127:   PetscCall(CheckSolution(sol, solRef));

129:   /* Clean up and finalize PETSc */
130:   PetscCall(KSPDestroy(&ksp));
131:   PetscCall(VecDestroy(&sol));
132:   PetscCall(VecDestroy(&solRef));
133:   PetscCall(VecDestroy(&rhs));
134:   PetscCall(MatDestroy(&A));
135:   PetscCall(DMDestroy(&dmSol));
136:   PetscCall(PetscFinalize());
137:   return 0;
138: }

140: /*
141: Note: this system is not well-scaled! Generally one would adjust the equations
142:  to try to get matrix entries to be of comparable order, regardless of grid spacing
143:  or choice of coefficients.
144: */
145: static PetscErrorCode CreateSystem(DM dmSol, Mat *pA, Vec *pRhs, PetscBool pinPressure)
146: {
147:   PetscInt      N[2];
148:   PetscInt      ex, ey, startx, starty, nx, ny;
149:   PetscInt      iprev, icenter, inext;
150:   Mat           A;
151:   Vec           rhs;
152:   PetscReal     hx, hy;
153:   PetscScalar **cArrX, **cArrY;

155:   /* Here, we showcase two different methods for manipulating local vector entries.
156:      One can use DMStagStencil objects with DMStagVecSetValuesStencil(),
157:      making sure to call VecAssemble[Begin/End]() after all values are set.
158:      Alternately, one can use DMStagVecGetArray[Read]() and DMStagVecRestoreArray[Read]().
159:      The first approach is used to build the rhs, and the second is used to
160:      obtain coordinate values. Working with the array is almost certainly more efficient,
161:      but only allows setting local entries, requires understanding which "slot" to use,
162:      and doesn't correspond as precisely to the matrix assembly process using DMStagStencil objects */

164:   PetscFunctionBeginUser;
165:   PetscCall(DMCreateMatrix(dmSol, pA));
166:   A = *pA;
167:   PetscCall(DMCreateGlobalVector(dmSol, pRhs));
168:   rhs = *pRhs;
169:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
170:   PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], NULL));
171:   hx = 1.0 / N[0];
172:   hy = 1.0 / N[1];
173:   PetscCall(DMStagGetProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
174:   PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, ELEMENT, &icenter));
175:   PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, LEFT, &iprev));
176:   PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, RIGHT, &inext));

178:   /* Loop over all local elements. Note that it may be more efficient in real
179:      applications to loop over each boundary separately */
180:   for (ey = starty; ey < starty + ny; ++ey) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
181:     for (ex = startx; ex < startx + nx; ++ex) {
182:       if (ex == N[0] - 1) {
183:         /* Right Boundary velocity Dirichlet */
184:         DMStagStencil     row;
185:         PetscScalar       valRhs;
186:         const PetscScalar valA = 1.0;
187:         row.i                  = ex;
188:         row.j                  = ey;
189:         row.loc                = RIGHT;
190:         row.c                  = 0;
191:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
192:         valRhs = uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
193:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
194:       }
195:       if (ey == N[1] - 1) {
196:         /* Top boundary velocity Dirichlet */
197:         DMStagStencil     row;
198:         PetscScalar       valRhs;
199:         const PetscScalar valA = 1.0;
200:         row.i                  = ex;
201:         row.j                  = ey;
202:         row.loc                = UP;
203:         row.c                  = 0;
204:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
205:         valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
206:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
207:       }

209:       if (ey == 0) {
210:         /* Bottom boundary velocity Dirichlet */
211:         DMStagStencil     row;
212:         PetscScalar       valRhs;
213:         const PetscScalar valA = 1.0;
214:         row.i                  = ex;
215:         row.j                  = ey;
216:         row.loc                = DOWN;
217:         row.c                  = 0;
218:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
219:         valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
220:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
221:       } else {
222:         /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
223:         DMStagStencil row, col[7];
224:         PetscScalar   valA[7], valRhs;
225:         PetscInt      nEntries;

227:         row.i   = ex;
228:         row.j   = ey;
229:         row.loc = DOWN;
230:         row.c   = 0;
231:         if (ex == 0) {
232:           nEntries   = 6;
233:           col[0].i   = ex;
234:           col[0].j   = ey;
235:           col[0].loc = DOWN;
236:           col[0].c   = 0;
237:           valA[0]    = -1.0 / (hx * hx) - 2.0 / (hy * hy);
238:           col[1].i   = ex;
239:           col[1].j   = ey - 1;
240:           col[1].loc = DOWN;
241:           col[1].c   = 0;
242:           valA[1]    = 1.0 / (hy * hy);
243:           col[2].i   = ex;
244:           col[2].j   = ey + 1;
245:           col[2].loc = DOWN;
246:           col[2].c   = 0;
247:           valA[2]    = 1.0 / (hy * hy);
248:           /* Missing left element */
249:           col[3].i   = ex + 1;
250:           col[3].j   = ey;
251:           col[3].loc = DOWN;
252:           col[3].c   = 0;
253:           valA[3]    = 1.0 / (hx * hx);
254:           col[4].i   = ex;
255:           col[4].j   = ey - 1;
256:           col[4].loc = ELEMENT;
257:           col[4].c   = 0;
258:           valA[4]    = 1.0 / hy;
259:           col[5].i   = ex;
260:           col[5].j   = ey;
261:           col[5].loc = ELEMENT;
262:           col[5].c   = 0;
263:           valA[5]    = -1.0 / hy;
264:         } else if (ex == N[0] - 1) {
265:           /* Right boundary y velocity stencil */
266:           nEntries   = 6;
267:           col[0].i   = ex;
268:           col[0].j   = ey;
269:           col[0].loc = DOWN;
270:           col[0].c   = 0;
271:           valA[0]    = -1.0 / (hx * hx) - 2.0 / (hy * hy);
272:           col[1].i   = ex;
273:           col[1].j   = ey - 1;
274:           col[1].loc = DOWN;
275:           col[1].c   = 0;
276:           valA[1]    = 1.0 / (hy * hy);
277:           col[2].i   = ex;
278:           col[2].j   = ey + 1;
279:           col[2].loc = DOWN;
280:           col[2].c   = 0;
281:           valA[2]    = 1.0 / (hy * hy);
282:           col[3].i   = ex - 1;
283:           col[3].j   = ey;
284:           col[3].loc = DOWN;
285:           col[3].c   = 0;
286:           valA[3]    = 1.0 / (hx * hx);
287:           /* Missing right element */
288:           col[4].i   = ex;
289:           col[4].j   = ey - 1;
290:           col[4].loc = ELEMENT;
291:           col[4].c   = 0;
292:           valA[4]    = 1.0 / hy;
293:           col[5].i   = ex;
294:           col[5].j   = ey;
295:           col[5].loc = ELEMENT;
296:           col[5].c   = 0;
297:           valA[5]    = -1.0 / hy;
298:         } else {
299:           nEntries   = 7;
300:           col[0].i   = ex;
301:           col[0].j   = ey;
302:           col[0].loc = DOWN;
303:           col[0].c   = 0;
304:           valA[0]    = -2.0 / (hx * hx) - 2.0 / (hy * hy);
305:           col[1].i   = ex;
306:           col[1].j   = ey - 1;
307:           col[1].loc = DOWN;
308:           col[1].c   = 0;
309:           valA[1]    = 1.0 / (hy * hy);
310:           col[2].i   = ex;
311:           col[2].j   = ey + 1;
312:           col[2].loc = DOWN;
313:           col[2].c   = 0;
314:           valA[2]    = 1.0 / (hy * hy);
315:           col[3].i   = ex - 1;
316:           col[3].j   = ey;
317:           col[3].loc = DOWN;
318:           col[3].c   = 0;
319:           valA[3]    = 1.0 / (hx * hx);
320:           col[4].i   = ex + 1;
321:           col[4].j   = ey;
322:           col[4].loc = DOWN;
323:           col[4].c   = 0;
324:           valA[4]    = 1.0 / (hx * hx);
325:           col[5].i   = ex;
326:           col[5].j   = ey - 1;
327:           col[5].loc = ELEMENT;
328:           col[5].c   = 0;
329:           valA[5]    = 1.0 / hy;
330:           col[6].i   = ex;
331:           col[6].j   = ey;
332:           col[6].loc = ELEMENT;
333:           col[6].c   = 0;
334:           valA[6]    = -1.0 / hy;
335:         }
336:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
337:         valRhs = fy(cArrX[ex][icenter], cArrY[ey][iprev]);
338:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
339:       }

341:       if (ex == 0) {
342:         /* Left velocity Dirichlet */
343:         DMStagStencil     row;
344:         PetscScalar       valRhs;
345:         const PetscScalar valA = 1.0;
346:         row.i                  = ex;
347:         row.j                  = ey;
348:         row.loc                = LEFT;
349:         row.c                  = 0;
350:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
351:         valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
352:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
353:       } else {
354:         /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
355:         DMStagStencil row, col[7];
356:         PetscScalar   valA[7], valRhs;
357:         PetscInt      nEntries;
358:         row.i   = ex;
359:         row.j   = ey;
360:         row.loc = LEFT;
361:         row.c   = 0;

363:         if (ey == 0) {
364:           nEntries   = 6;
365:           col[0].i   = ex;
366:           col[0].j   = ey;
367:           col[0].loc = LEFT;
368:           col[0].c   = 0;
369:           valA[0]    = -2.0 / (hx * hx) - 1.0 / (hy * hy);
370:           /* missing term from element below */
371:           col[1].i   = ex;
372:           col[1].j   = ey + 1;
373:           col[1].loc = LEFT;
374:           col[1].c   = 0;
375:           valA[1]    = 1.0 / (hy * hy);
376:           col[2].i   = ex - 1;
377:           col[2].j   = ey;
378:           col[2].loc = LEFT;
379:           col[2].c   = 0;
380:           valA[2]    = 1.0 / (hx * hx);
381:           col[3].i   = ex + 1;
382:           col[3].j   = ey;
383:           col[3].loc = LEFT;
384:           col[3].c   = 0;
385:           valA[3]    = 1.0 / (hx * hx);
386:           col[4].i   = ex - 1;
387:           col[4].j   = ey;
388:           col[4].loc = ELEMENT;
389:           col[4].c   = 0;
390:           valA[4]    = 1.0 / hx;
391:           col[5].i   = ex;
392:           col[5].j   = ey;
393:           col[5].loc = ELEMENT;
394:           col[5].c   = 0;
395:           valA[5]    = -1.0 / hx;
396:         } else if (ey == N[1] - 1) {
397:           /* Top boundary x velocity stencil */
398:           nEntries   = 6;
399:           row.i      = ex;
400:           row.j      = ey;
401:           row.loc    = LEFT;
402:           row.c      = 0;
403:           col[0].i   = ex;
404:           col[0].j   = ey;
405:           col[0].loc = LEFT;
406:           col[0].c   = 0;
407:           valA[0]    = -2.0 / (hx * hx) - 1.0 / (hy * hy);
408:           col[1].i   = ex;
409:           col[1].j   = ey - 1;
410:           col[1].loc = LEFT;
411:           col[1].c   = 0;
412:           valA[1]    = 1.0 / (hy * hy);
413:           /* Missing element above term */
414:           col[2].i   = ex - 1;
415:           col[2].j   = ey;
416:           col[2].loc = LEFT;
417:           col[2].c   = 0;
418:           valA[2]    = 1.0 / (hx * hx);
419:           col[3].i   = ex + 1;
420:           col[3].j   = ey;
421:           col[3].loc = LEFT;
422:           col[3].c   = 0;
423:           valA[3]    = 1.0 / (hx * hx);
424:           col[4].i   = ex - 1;
425:           col[4].j   = ey;
426:           col[4].loc = ELEMENT;
427:           col[4].c   = 0;
428:           valA[4]    = 1.0 / hx;
429:           col[5].i   = ex;
430:           col[5].j   = ey;
431:           col[5].loc = ELEMENT;
432:           col[5].c   = 0;
433:           valA[5]    = -1.0 / hx;
434:         } else {
435:           /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
436:           nEntries   = 7;
437:           col[0].i   = ex;
438:           col[0].j   = ey;
439:           col[0].loc = LEFT;
440:           col[0].c   = 0;
441:           valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy);
442:           col[1].i   = ex;
443:           col[1].j   = ey - 1;
444:           col[1].loc = LEFT;
445:           col[1].c   = 0;
446:           valA[1]    = 1.0 / (hy * hy);
447:           col[2].i   = ex;
448:           col[2].j   = ey + 1;
449:           col[2].loc = LEFT;
450:           col[2].c   = 0;
451:           valA[2]    = 1.0 / (hy * hy);
452:           col[3].i   = ex - 1;
453:           col[3].j   = ey;
454:           col[3].loc = LEFT;
455:           col[3].c   = 0;
456:           valA[3]    = 1.0 / (hx * hx);
457:           col[4].i   = ex + 1;
458:           col[4].j   = ey;
459:           col[4].loc = LEFT;
460:           col[4].c   = 0;
461:           valA[4]    = 1.0 / (hx * hx);
462:           col[5].i   = ex - 1;
463:           col[5].j   = ey;
464:           col[5].loc = ELEMENT;
465:           col[5].c   = 0;
466:           valA[5]    = 1.0 / hx;
467:           col[6].i   = ex;
468:           col[6].j   = ey;
469:           col[6].loc = ELEMENT;
470:           col[6].c   = 0;
471:           valA[6]    = -1.0 / hx;
472:         }
473:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
474:         valRhs = fx(cArrX[ex][iprev], cArrY[ey][icenter]);
475:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
476:       }

478:       /* P equation : u_x + v_y = g
479:          Note that this includes an explicit zero on the diagonal. This is only needed for
480:          direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
481:       if (pinPressure && ex == 0 && ey == 0) { /* Pin the first pressure node, if requested */
482:         DMStagStencil row;
483:         PetscScalar   valA, valRhs;
484:         row.i   = ex;
485:         row.j   = ey;
486:         row.loc = ELEMENT;
487:         row.c   = 0;
488:         valA    = 1.0;
489:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
490:         valRhs = pRef(cArrX[ex][icenter], cArrY[ey][icenter]);
491:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
492:       } else {
493:         DMStagStencil row, col[5];
494:         PetscScalar   valA[5], valRhs;

496:         row.i      = ex;
497:         row.j      = ey;
498:         row.loc    = ELEMENT;
499:         row.c      = 0;
500:         col[0].i   = ex;
501:         col[0].j   = ey;
502:         col[0].loc = LEFT;
503:         col[0].c   = 0;
504:         valA[0]    = -1.0 / hx;
505:         col[1].i   = ex;
506:         col[1].j   = ey;
507:         col[1].loc = RIGHT;
508:         col[1].c   = 0;
509:         valA[1]    = 1.0 / hx;
510:         col[2].i   = ex;
511:         col[2].j   = ey;
512:         col[2].loc = DOWN;
513:         col[2].c   = 0;
514:         valA[2]    = -1.0 / hy;
515:         col[3].i   = ex;
516:         col[3].j   = ey;
517:         col[3].loc = UP;
518:         col[3].c   = 0;
519:         valA[3]    = 1.0 / hy;
520:         col[4]     = row;
521:         valA[4]    = 0.0;
522:         PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 5, col, valA, INSERT_VALUES));
523:         valRhs = g(cArrX[ex][icenter], cArrY[ey][icenter]);
524:         PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
525:       }
526:     }
527:   }
528:   PetscCall(DMStagRestoreProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
529:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
530:   PetscCall(VecAssemblyBegin(rhs));
531:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
532:   PetscCall(VecAssemblyEnd(rhs));
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: /* Create a pressure-only DMStag and use it to generate a nullspace vector
537:    - Create a compatible DMStag with one dof per element (and nothing else).
538:    - Create a constant vector and normalize it
539:    - Migrate it to a vector on the original dmSol (making use of the fact
540:    that this will fill in zeros for "extra" dof)
541:    - Set the nullspace for the operator
542:    - Destroy everything (the operator keeps the references it needs) */
543: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
544: {
545:   DM           dmPressure;
546:   Vec          constantPressure, basis;
547:   PetscReal    nrm;
548:   MatNullSpace matNullSpace;

550:   PetscFunctionBeginUser;
551:   PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 0, 1, 0, &dmPressure));
552:   PetscCall(DMGetGlobalVector(dmPressure, &constantPressure));
553:   PetscCall(VecSet(constantPressure, 1.0));
554:   PetscCall(VecNorm(constantPressure, NORM_2, &nrm));
555:   PetscCall(VecScale(constantPressure, 1.0 / nrm));
556:   PetscCall(DMCreateGlobalVector(dmSol, &basis));
557:   PetscCall(DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis));
558:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace));
559:   PetscCall(VecDestroy(&basis));
560:   PetscCall(VecDestroy(&constantPressure));
561:   PetscCall(MatSetNullSpace(A, matNullSpace));
562:   PetscCall(MatNullSpaceDestroy(&matNullSpace));
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: /* Create a reference solution.
567:    Here, we use the more direct method of iterating over arrays.  */
568: static PetscErrorCode CreateReferenceSolution(DM dmSol, Vec *pSolRef)
569: {
570:   PetscInt       startx, starty, nx, ny, nExtra[2], ex, ey;
571:   PetscInt       iuy, iux, ip, iprev, icenter;
572:   PetscScalar ***arrSol, **cArrX, **cArrY;
573:   Vec            solRefLocal;

575:   PetscFunctionBeginUser;
576:   PetscCall(DMCreateGlobalVector(dmSol, pSolRef));
577:   PetscCall(DMGetLocalVector(dmSol, &solRefLocal));

579:   /* Obtain indices to use in the raw arrays */
580:   PetscCall(DMStagGetLocationSlot(dmSol, DOWN, 0, &iuy));
581:   PetscCall(DMStagGetLocationSlot(dmSol, LEFT, 0, &iux));
582:   PetscCall(DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip));

584:   /* Use high-level convenience functions to get raw arrays and indices for 1d coordinates */
585:   PetscCall(DMStagGetProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
586:   PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, ELEMENT, &icenter));
587:   PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, LEFT, &iprev));

589:   PetscCall(DMStagVecGetArray(dmSol, solRefLocal, &arrSol));
590:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, NULL, &nx, &ny, NULL, &nExtra[0], &nExtra[1], NULL));
591:   for (ey = starty; ey < starty + ny + nExtra[1]; ++ey) {
592:     for (ex = startx; ex < startx + nx + nExtra[0]; ++ex) {
593:       arrSol[ey][ex][iuy] = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
594:       arrSol[ey][ex][iux] = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
595:       if (ey < starty + ny && ex < startx + nx) { /* Don't fill on the dummy elements (though you could, and these values would just be ignored) */
596:         arrSol[ey][ex][ip] = pRef(cArrX[ex][icenter], cArrY[ey][icenter]);
597:       }
598:     }
599:   }
600:   PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
601:   PetscCall(DMStagRestoreProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
602:   PetscCall(DMLocalToGlobal(dmSol, solRefLocal, INSERT_VALUES, *pSolRef));
603:   PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
604:   PetscFunctionReturn(PETSC_SUCCESS);
605: }

607: static PetscErrorCode CheckSolution(Vec sol, Vec solRef)
608: {
609:   Vec       diff;
610:   PetscReal normsolRef, errAbs, errRel;

612:   PetscFunctionBeginUser;
613:   PetscCall(VecDuplicate(sol, &diff));
614:   PetscCall(VecCopy(sol, diff));
615:   PetscCall(VecAXPY(diff, -1.0, solRef));
616:   PetscCall(VecNorm(diff, NORM_2, &errAbs));
617:   PetscCall(VecNorm(solRef, NORM_2, &normsolRef));
618:   errRel = errAbs / normsolRef;
619:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel));
620:   PetscCall(VecDestroy(&diff));
621:   PetscFunctionReturn(PETSC_SUCCESS);
622: }

624: /*TEST

626:    test:
627:       suffix: 1
628:       nsize: 4
629:       args: -ksp_monitor_short -ksp_converged_reason

631:    test:
632:       suffix: direct_umfpack
633:       requires: suitesparse
634:       nsize: 1
635:       args: -pinpressure 1 -stag_grid_x 8 -stag_grid_y 6 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack

637: TEST*/