Actual source code: ex15.c

  1: static char help[] = "Test DMStag default MG components, on a Stokes-like system.\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmstag.h>
  5: #include <petscksp.h>

  7: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);

  9: int main(int argc, char **argv)
 10: {
 11:   DM        dm;
 12:   PetscInt  dim;
 13:   PetscBool flg;
 14:   KSP       ksp;
 15:   Mat       A;
 16:   Vec       x, b;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 20:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, &flg));
 21:   if (!flg) {
 22:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Supply -dim option\n"));
 23:     return 1;
 24:   }
 25:   if (dim == 1) {
 26:     PetscCall(DMStagCreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, &dm));
 27:   } else if (dim == 2) {
 28:     PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dm));
 29:   } else if (dim == 3) {
 30:     PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 0, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, NULL, &dm));
 31:   } else {
 32:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Supply -dim option with value 1, 2, or 3\n"));
 33:     return 1;
 34:   }
 35:   PetscCall(DMSetFromOptions(dm));
 36:   PetscCall(DMSetUp(dm));

 38:   /* Directly create a coarsened DM and transfer operators */
 39:   {
 40:     DM dmCoarse;
 41:     PetscCall(DMCoarsen(dm, MPI_COMM_NULL, &dmCoarse));
 42:     {
 43:       Mat Ai;
 44:       PetscCall(DMCreateInterpolation(dmCoarse, dm, &Ai, NULL));
 45:       PetscCall(MatDestroy(&Ai));
 46:     }
 47:     {
 48:       Mat Ar;
 49:       PetscCall(DMCreateRestriction(dmCoarse, dm, &Ar));
 50:       PetscCall(MatDestroy(&Ar));
 51:     }
 52:     PetscCall(DMDestroy(&dmCoarse));
 53:   }

 55:   /* Create and solve a system */
 56:   PetscCall(CreateSystem(dm, &A, &b));
 57:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
 58:   PetscCall(KSPSetOperators(ksp, A, A));
 59:   PetscCall(KSPSetDM(ksp, dm));
 60:   PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
 61:   PetscCall(KSPSetFromOptions(ksp));

 63:   PetscCall(VecDuplicate(b, &x));
 64:   PetscCall(KSPSolve(ksp, b, x));

 66:   PetscCall(VecDestroy(&x));
 67:   PetscCall(VecDestroy(&b));
 68:   PetscCall(KSPDestroy(&ksp));
 69:   PetscCall(MatDestroy(&A));
 70:   PetscCall(DMDestroy(&dm));
 71:   PetscCall(PetscFinalize());
 72:   return 0;
 73: }

 75: /* Note: unlike in the 2D case, this does not include reasonable scaling and so will not work well */
 76: PetscErrorCode CreateSystem1d(DM dm, Mat *A, Vec *b)
 77: {
 78:   PetscInt dim;

 80:   PetscFunctionBeginUser;
 81:   PetscCall(DMGetDimension(dm, &dim));
 82:   PetscCall(DMCreateMatrix(dm, A));
 83:   PetscCall(DMCreateGlobalVector(dm, b));
 84:   if (dim == 1) {
 85:     PetscInt    e, start, n, N;
 86:     PetscBool   isFirstRank, isLastRank;
 87:     PetscScalar h;
 88:     PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, NULL, NULL, NULL));
 89:     PetscCall(DMStagGetGlobalSizes(dm, &N, NULL, NULL));
 90:     h = 1.0 / N;
 91:     PetscCall(DMStagGetIsFirstRank(dm, &isFirstRank, NULL, NULL));
 92:     PetscCall(DMStagGetIsLastRank(dm, &isLastRank, NULL, NULL));
 93:     for (e = start; e < start + n; ++e) {
 94:       DMStagStencil pos[3];
 95:       PetscScalar   val[3];
 96:       PetscInt      idxLoc;

 98:       if (isFirstRank && e == start) {
 99:         /* Fix first pressure node to eliminate nullspace */
100:         idxLoc          = 0;
101:         pos[idxLoc].i   = e;
102:         pos[idxLoc].loc = DMSTAG_ELEMENT;
103:         pos[idxLoc].c   = 0;
104:         val[idxLoc]     = 1.0; /* 0 pressure forcing term (physical) */
105:         ++idxLoc;
106:       } else {
107:         idxLoc          = 0;
108:         pos[idxLoc].i   = e;
109:         pos[idxLoc].loc = DMSTAG_ELEMENT;
110:         pos[idxLoc].c   = 0;
111:         val[idxLoc]     = 0.0; /* 0 pressure forcing term (physical) */
112:         ++idxLoc;
113:       }

115:       if (isFirstRank && e == start) {
116:         pos[idxLoc].i   = e;
117:         pos[idxLoc].loc = DMSTAG_LEFT;
118:         pos[idxLoc].c   = 0;
119:         val[idxLoc]     = 3.0; /* fixed left BC */
120:         ++idxLoc;
121:       } else {
122:         pos[idxLoc].i   = e;
123:         pos[idxLoc].loc = DMSTAG_LEFT;
124:         pos[idxLoc].c   = 0;
125:         val[idxLoc]     = 1.0; /* constant forcing */
126:         ++idxLoc;
127:       }
128:       if (isLastRank && e == start + n - 1) {
129:         /* Special case on right boundary (in addition to usual case) */
130:         pos[idxLoc].i   = e; /* This element in the 1d ordering */
131:         pos[idxLoc].loc = DMSTAG_RIGHT;
132:         pos[idxLoc].c   = 0;
133:         val[idxLoc]     = 3.0; /* fixed right BC */
134:         ++idxLoc;
135:       }
136:       PetscCall(DMStagVecSetValuesStencil(dm, *b, idxLoc, pos, val, INSERT_VALUES));
137:     }

139:     for (e = start; e < start + n; ++e) {
140:       if (isFirstRank && e == start) {
141:         DMStagStencil row;
142:         PetscScalar   val;

144:         row.i   = e;
145:         row.loc = DMSTAG_LEFT;
146:         row.c   = 0;
147:         val     = 1.0;
148:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
149:       } else {
150:         DMStagStencil row, col[5];
151:         PetscScalar   val[5];

153:         row.i   = e;
154:         row.loc = DMSTAG_LEFT;
155:         row.c   = 0;

157:         col[0].i   = e;
158:         col[0].loc = DMSTAG_ELEMENT;
159:         col[0].c   = 0;

161:         col[1].i   = e - 1;
162:         col[1].loc = DMSTAG_ELEMENT;
163:         col[1].c   = 0;

165:         val[0] = -1.0 / h;
166:         val[1] = 1.0 / h;

168:         col[2].i   = e;
169:         col[2].loc = DMSTAG_LEFT;
170:         col[2].c   = 0;
171:         val[2]     = -2.0 / (h * h);

173:         col[3].i   = e - 1;
174:         col[3].loc = DMSTAG_LEFT;
175:         col[3].c   = 0;
176:         val[3]     = 1.0 / (h * h);

178:         col[4].i   = e + 1;
179:         col[4].loc = DMSTAG_LEFT;
180:         col[4].c   = 0;
181:         val[4]     = 1.0 / (h * h);

183:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 5, col, val, INSERT_VALUES));
184:       }

186:       /* Additional velocity point (BC) on the right */
187:       if (isLastRank && e == start + n - 1) {
188:         DMStagStencil row;
189:         PetscScalar   val;

191:         row.i   = e;
192:         row.loc = DMSTAG_RIGHT;
193:         row.c   = 0;
194:         val     = 1.0;
195:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
196:       }

198:       /* Equation on pressure (element) variables */
199:       if (isFirstRank && e == 0) {
200:         /* Fix first pressure node to eliminate nullspace */
201:         DMStagStencil row;
202:         PetscScalar   val;

204:         row.i   = e;
205:         row.loc = DMSTAG_ELEMENT;
206:         row.c   = 0;
207:         val     = 1.0;

209:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
210:       } else {
211:         DMStagStencil row, col[2];
212:         PetscScalar   val[2];

214:         row.i   = e;
215:         row.loc = DMSTAG_ELEMENT;
216:         row.c   = 0;

218:         col[0].i   = e;
219:         col[0].loc = DMSTAG_LEFT;
220:         col[0].c   = 0;

222:         col[1].i   = e;
223:         col[1].loc = DMSTAG_RIGHT;
224:         col[1].c   = 0;

226:         val[0] = -1.0 / h;
227:         val[1] = 1.0 / h;

229:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 2, col, val, INSERT_VALUES));
230:       }
231:     }
232:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
233:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
234:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
235:   PetscCall(VecAssemblyBegin(*b));
236:   PetscCall(VecAssemblyEnd(*b));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: PetscErrorCode CreateSystem2d(DM dm, Mat *A, Vec *b)
241: {
242:   PetscInt  N[2];
243:   PetscBool isLastRankx, isLastRanky, isFirstRankx, isFirstRanky;
244:   PetscInt  ex, ey, startx, starty, nx, ny;
245:   PetscReal hx, hy, dv;

247:   PetscFunctionBeginUser;
248:   PetscCall(DMCreateMatrix(dm, A));
249:   PetscCall(DMCreateGlobalVector(dm, b));
250:   PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
251:   PetscCall(DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL));
252:   PetscCall(DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL));
253:   PetscCall(DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL));
254:   hx = 1.0 / N[0];
255:   hy = 1.0 / N[1];
256:   dv = hx * hy;

258:   for (ey = starty; ey < starty + ny; ++ey) {
259:     for (ex = startx; ex < startx + nx; ++ex) {
260:       if (ex == N[0] - 1) {
261:         /* Right Boundary velocity Dirichlet */
262:         DMStagStencil     row;
263:         PetscScalar       valRhs;
264:         const PetscScalar valA = 1.0;
265:         row.i                  = ex;
266:         row.j                  = ey;
267:         row.loc                = DMSTAG_RIGHT;
268:         row.c                  = 0;
269:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
270:         valRhs = 0.0; /* zero Dirichlet */
271:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
272:       }
273:       if (ey == N[1] - 1) {
274:         /* Top boundary velocity Dirichlet */
275:         DMStagStencil     row;
276:         PetscScalar       valRhs;
277:         const PetscScalar valA = 1.0;
278:         row.i                  = ex;
279:         row.j                  = ey;
280:         row.loc                = DMSTAG_UP;
281:         row.c                  = 0;
282:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
283:         valRhs = 0.0; /* zero Diri */
284:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
285:       }

287:       if (ey == 0) {
288:         /* Bottom boundary velocity Dirichlet */
289:         DMStagStencil     row;
290:         PetscScalar       valRhs;
291:         const PetscScalar valA = 1.0;
292:         row.i                  = ex;
293:         row.j                  = ey;
294:         row.loc                = DMSTAG_DOWN;
295:         row.c                  = 0;
296:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
297:         valRhs = 0.0; /* zero Diri */
298:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
299:       } else {
300:         /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
301:         DMStagStencil row, col[7];
302:         PetscScalar   valA[7], valRhs;
303:         PetscInt      nEntries;

305:         row.i   = ex;
306:         row.j   = ey;
307:         row.loc = DMSTAG_DOWN;
308:         row.c   = 0;
309:         if (ex == 0) {
310:           nEntries   = 6;
311:           col[0].i   = ex;
312:           col[0].j   = ey;
313:           col[0].loc = DMSTAG_DOWN;
314:           col[0].c   = 0;
315:           valA[0]    = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
316:           col[1].i   = ex;
317:           col[1].j   = ey - 1;
318:           col[1].loc = DMSTAG_DOWN;
319:           col[1].c   = 0;
320:           valA[1]    = dv * 1.0 / (hy * hy);
321:           col[2].i   = ex;
322:           col[2].j   = ey + 1;
323:           col[2].loc = DMSTAG_DOWN;
324:           col[2].c   = 0;
325:           valA[2]    = dv * 1.0 / (hy * hy);
326:           /* Missing left element */
327:           col[3].i   = ex + 1;
328:           col[3].j   = ey;
329:           col[3].loc = DMSTAG_DOWN;
330:           col[3].c   = 0;
331:           valA[3]    = dv * 1.0 / (hx * hx);
332:           col[4].i   = ex;
333:           col[4].j   = ey - 1;
334:           col[4].loc = DMSTAG_ELEMENT;
335:           col[4].c   = 0;
336:           valA[4]    = dv * 1.0 / hy;
337:           col[5].i   = ex;
338:           col[5].j   = ey;
339:           col[5].loc = DMSTAG_ELEMENT;
340:           col[5].c   = 0;
341:           valA[5]    = -dv * 1.0 / hy;
342:         } else if (ex == N[0] - 1) {
343:           /* Right boundary y velocity stencil */
344:           nEntries   = 6;
345:           col[0].i   = ex;
346:           col[0].j   = ey;
347:           col[0].loc = DMSTAG_DOWN;
348:           col[0].c   = 0;
349:           valA[0]    = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
350:           col[1].i   = ex;
351:           col[1].j   = ey - 1;
352:           col[1].loc = DMSTAG_DOWN;
353:           col[1].c   = 0;
354:           valA[1]    = dv * 1.0 / (hy * hy);
355:           col[2].i   = ex;
356:           col[2].j   = ey + 1;
357:           col[2].loc = DMSTAG_DOWN;
358:           col[2].c   = 0;
359:           valA[2]    = dv * 1.0 / (hy * hy);
360:           col[3].i   = ex - 1;
361:           col[3].j   = ey;
362:           col[3].loc = DMSTAG_DOWN;
363:           col[3].c   = 0;
364:           valA[3]    = dv * 1.0 / (hx * hx);
365:           /* Missing right element */
366:           col[4].i   = ex;
367:           col[4].j   = ey - 1;
368:           col[4].loc = DMSTAG_ELEMENT;
369:           col[4].c   = 0;
370:           valA[4]    = dv * 1.0 / hy;
371:           col[5].i   = ex;
372:           col[5].j   = ey;
373:           col[5].loc = DMSTAG_ELEMENT;
374:           col[5].c   = 0;
375:           valA[5]    = -dv * 1.0 / hy;
376:         } else {
377:           nEntries   = 7;
378:           col[0].i   = ex;
379:           col[0].j   = ey;
380:           col[0].loc = DMSTAG_DOWN;
381:           col[0].c   = 0;
382:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
383:           col[1].i   = ex;
384:           col[1].j   = ey - 1;
385:           col[1].loc = DMSTAG_DOWN;
386:           col[1].c   = 0;
387:           valA[1]    = dv * 1.0 / (hy * hy);
388:           col[2].i   = ex;
389:           col[2].j   = ey + 1;
390:           col[2].loc = DMSTAG_DOWN;
391:           col[2].c   = 0;
392:           valA[2]    = dv * 1.0 / (hy * hy);
393:           col[3].i   = ex - 1;
394:           col[3].j   = ey;
395:           col[3].loc = DMSTAG_DOWN;
396:           col[3].c   = 0;
397:           valA[3]    = dv * 1.0 / (hx * hx);
398:           col[4].i   = ex + 1;
399:           col[4].j   = ey;
400:           col[4].loc = DMSTAG_DOWN;
401:           col[4].c   = 0;
402:           valA[4]    = dv * 1.0 / (hx * hx);
403:           col[5].i   = ex;
404:           col[5].j   = ey - 1;
405:           col[5].loc = DMSTAG_ELEMENT;
406:           col[5].c   = 0;
407:           valA[5]    = dv * 1.0 / hy;
408:           col[6].i   = ex;
409:           col[6].j   = ey;
410:           col[6].loc = DMSTAG_ELEMENT;
411:           col[6].c   = 0;
412:           valA[6]    = -dv * 1.0 / hy;
413:         }
414:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, nEntries, col, valA, INSERT_VALUES));
415:         valRhs = dv * 1.0; /* non-zero */
416:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
417:       }

419:       if (ex == 0) {
420:         /* Left velocity Dirichlet */
421:         DMStagStencil     row;
422:         PetscScalar       valRhs;
423:         const PetscScalar valA = 1.0;
424:         row.i                  = ex;
425:         row.j                  = ey;
426:         row.loc                = DMSTAG_LEFT;
427:         row.c                  = 0;
428:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
429:         valRhs = 0.0; /* zero Diri */
430:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
431:       } else {
432:         /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
433:         DMStagStencil row, col[7];
434:         PetscScalar   valA[7], valRhs;
435:         PetscInt      nEntries;
436:         row.i   = ex;
437:         row.j   = ey;
438:         row.loc = DMSTAG_LEFT;
439:         row.c   = 0;

441:         if (ey == 0) {
442:           nEntries   = 6;
443:           col[0].i   = ex;
444:           col[0].j   = ey;
445:           col[0].loc = DMSTAG_LEFT;
446:           col[0].c   = 0;
447:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
448:           /* missing term from element below */
449:           col[1].i   = ex;
450:           col[1].j   = ey + 1;
451:           col[1].loc = DMSTAG_LEFT;
452:           col[1].c   = 0;
453:           valA[1]    = dv * 1.0 / (hy * hy);
454:           col[2].i   = ex - 1;
455:           col[2].j   = ey;
456:           col[2].loc = DMSTAG_LEFT;
457:           col[2].c   = 0;
458:           valA[2]    = dv * 1.0 / (hx * hx);
459:           col[3].i   = ex + 1;
460:           col[3].j   = ey;
461:           col[3].loc = DMSTAG_LEFT;
462:           col[3].c   = 0;
463:           valA[3]    = dv * 1.0 / (hx * hx);
464:           col[4].i   = ex - 1;
465:           col[4].j   = ey;
466:           col[4].loc = DMSTAG_ELEMENT;
467:           col[4].c   = 0;
468:           valA[4]    = dv * 1.0 / hx;
469:           col[5].i   = ex;
470:           col[5].j   = ey;
471:           col[5].loc = DMSTAG_ELEMENT;
472:           col[5].c   = 0;
473:           valA[5]    = -dv * 1.0 / hx;
474:         } else if (ey == N[1] - 1) {
475:           /* Top boundary x velocity stencil */
476:           nEntries   = 6;
477:           row.i      = ex;
478:           row.j      = ey;
479:           row.loc    = DMSTAG_LEFT;
480:           row.c      = 0;
481:           col[0].i   = ex;
482:           col[0].j   = ey;
483:           col[0].loc = DMSTAG_LEFT;
484:           col[0].c   = 0;
485:           valA[0]    = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
486:           col[1].i   = ex;
487:           col[1].j   = ey - 1;
488:           col[1].loc = DMSTAG_LEFT;
489:           col[1].c   = 0;
490:           valA[1]    = dv * 1.0 / (hy * hy);
491:           /* Missing element above term */
492:           col[2].i   = ex - 1;
493:           col[2].j   = ey;
494:           col[2].loc = DMSTAG_LEFT;
495:           col[2].c   = 0;
496:           valA[2]    = dv * 1.0 / (hx * hx);
497:           col[3].i   = ex + 1;
498:           col[3].j   = ey;
499:           col[3].loc = DMSTAG_LEFT;
500:           col[3].c   = 0;
501:           valA[3]    = dv * 1.0 / (hx * hx);
502:           col[4].i   = ex - 1;
503:           col[4].j   = ey;
504:           col[4].loc = DMSTAG_ELEMENT;
505:           col[4].c   = 0;
506:           valA[4]    = dv * 1.0 / hx;
507:           col[5].i   = ex;
508:           col[5].j   = ey;
509:           col[5].loc = DMSTAG_ELEMENT;
510:           col[5].c   = 0;
511:           valA[5]    = -dv * 1.0 / hx;
512:         } else {
513:           nEntries   = 7;
514:           col[0].i   = ex;
515:           col[0].j   = ey;
516:           col[0].loc = DMSTAG_LEFT;
517:           col[0].c   = 0;
518:           valA[0]    = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
519:           col[1].i   = ex;
520:           col[1].j   = ey - 1;
521:           col[1].loc = DMSTAG_LEFT;
522:           col[1].c   = 0;
523:           valA[1]    = dv * 1.0 / (hy * hy);
524:           col[2].i   = ex;
525:           col[2].j   = ey + 1;
526:           col[2].loc = DMSTAG_LEFT;
527:           col[2].c   = 0;
528:           valA[2]    = dv * 1.0 / (hy * hy);
529:           col[3].i   = ex - 1;
530:           col[3].j   = ey;
531:           col[3].loc = DMSTAG_LEFT;
532:           col[3].c   = 0;
533:           valA[3]    = dv * 1.0 / (hx * hx);
534:           col[4].i   = ex + 1;
535:           col[4].j   = ey;
536:           col[4].loc = DMSTAG_LEFT;
537:           col[4].c   = 0;
538:           valA[4]    = dv * 1.0 / (hx * hx);
539:           col[5].i   = ex - 1;
540:           col[5].j   = ey;
541:           col[5].loc = DMSTAG_ELEMENT;
542:           col[5].c   = 0;
543:           valA[5]    = dv * 1.0 / hx;
544:           col[6].i   = ex;
545:           col[6].j   = ey;
546:           col[6].loc = DMSTAG_ELEMENT;
547:           col[6].c   = 0;
548:           valA[6]    = -dv * 1.0 / hx;
549:         }
550:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, nEntries, col, valA, INSERT_VALUES));
551:         valRhs = dv * 0.0; /* zero */
552:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
553:       }

555:       /* P equation : u_x + v_y = g
556:          Note that this includes an explicit zero on the diagonal. This is only needed for
557:          direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
558:       if (ex == 0 && ey == 0) { /* Pin the first pressure node */
559:         DMStagStencil row;
560:         PetscScalar   valA, valRhs;
561:         row.i   = ex;
562:         row.j   = ey;
563:         row.loc = DMSTAG_ELEMENT;
564:         row.c   = 0;
565:         valA    = 1.0;
566:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
567:         valRhs = 0.0; /* zero pinned pressure */
568:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
569:       } else {
570:         DMStagStencil row, col[5];
571:         PetscScalar   valA[5], valRhs;

573:         /* Note: the scaling by dv here is not principled and likely suboptimal */
574:         row.i      = ex;
575:         row.j      = ey;
576:         row.loc    = DMSTAG_ELEMENT;
577:         row.c      = 0;
578:         col[0].i   = ex;
579:         col[0].j   = ey;
580:         col[0].loc = DMSTAG_LEFT;
581:         col[0].c   = 0;
582:         valA[0]    = -dv * 1.0 / hx;
583:         col[1].i   = ex;
584:         col[1].j   = ey;
585:         col[1].loc = DMSTAG_RIGHT;
586:         col[1].c   = 0;
587:         valA[1]    = dv * 1.0 / hx;
588:         col[2].i   = ex;
589:         col[2].j   = ey;
590:         col[2].loc = DMSTAG_DOWN;
591:         col[2].c   = 0;
592:         valA[2]    = -dv * 1.0 / hy;
593:         col[3].i   = ex;
594:         col[3].j   = ey;
595:         col[3].loc = DMSTAG_UP;
596:         col[3].c   = 0;
597:         valA[3]    = dv * 1.0 / hy;
598:         col[4]     = row;
599:         valA[4]    = 0.0;
600:         PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 5, col, valA, INSERT_VALUES));
601:         valRhs = dv * 0.0; /* zero pressure forcing */
602:         PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
603:       }
604:     }
605:   }
606:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
607:   PetscCall(VecAssemblyBegin(*b));
608:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
609:   PetscCall(VecAssemblyEnd(*b));
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

613: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b)
614: {
615:   PetscInt dim;

617:   PetscFunctionBeginUser;
618:   PetscCall(DMGetDimension(dm, &dim));
619:   if (dim == 1) {
620:     PetscCall(CreateSystem1d(dm, A, b));
621:   } else if (dim == 2) {
622:     PetscCall(CreateSystem2d(dm, A, b));
623:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
624:   PetscFunctionReturn(PETSC_SUCCESS);
625: }

627: /*TEST

629:    test:
630:       suffix: 1d_directsmooth_seq
631:       nsize: 1
632:       requires: suitesparse
633:       args: -dim 1 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type umfpack -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason

635:    test:
636:       suffix: 1d_directsmooth_par
637:       nsize: 4
638:       requires: mumps
639:       args: -dim 1 /ex15 -dim 1 -stag_grid_x 16 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type mumps -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason

641:    test:
642:       suffix: 1d_fssmooth_seq
643:       nsize: 1
644:       requires: suitesparse
645:       args: -dim 1 -stag_grid_x 256 -ksp_converged_reason -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point

647:    test:
648:       suffix: 1d_fssmooth_par
649:       nsize: 1
650:       requires: mumps
651:       args: -dim 1 -stag_grid_x 256 -ksp_converged_reason -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point

653:    test:
654:       suffix: 2d_directsmooth_seq
655:       nsize: 1
656:       requires: suitesparse
657:       args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type umfpack -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason

659:    test:
660:       suffix: 2d_directsmooth_par
661:       nsize: 4
662:       requires: mumps
663:       args: -dim 1 /ex15 -dim 1 -stag_grid_x 16 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type mumps -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason

665:    test:
666:       suffix: 2d_fssmooth_seq
667:       nsize: 1
668:       requires: suitesparse
669:       args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason

671:    test:
672:       suffix: 2d_fssmooth_par
673:       nsize: 1
674:       requires: mumps
675:       args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason

677:    test:
678:       suffix: 2d_mono_mg
679:       nsize: 1
680:       requires: suitesparse
681:       args: -dim 2 -pc_type mg -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_type SCHUR -ksp_monitor_true_residual -ksp_converged_reason -mg_levels_fieldsplit_element_pc_type jacobi -mg_levels_fieldsplit_face_pc_type jacobi -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -pc_mg_levels 3 -mg_levels_pc_fieldsplit_schur_fact_type UPPER  -mg_levels_fieldsplit_face_pc_type sor -mg_levels_fieldsplit_face_ksp_type richardson -mg_levels_fieldsplit_face_pc_sor_symmetric -mg_levels_fieldsplit_element_ksp_type richardson  -mg_levels_fieldsplit_element_pc_type none -ksp_type fgmres -stag_grid_x 48 -stag_grid_y 48 -mg_levels_fieldsplit_face_ksp_max_it 1 -mg_levels_ksp_max_it 4 -mg_levels_fieldsplit_element_ksp_type preonly -mg_levels_fieldsplit_element_pc_type none -pc_mg_cycle_type w -mg_levels_ksp_max_it 4 -mg_levels_2_ksp_max_it 8

683: TEST*/