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