Actual source code: ex3.c

      1: static char help[] = "Solve a toy 3D 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 3D domain.

      7:   u_{xx} + u_{yy} + u_{zz} - p_x = f^x
      8:   v_{xx} + v_{yy} + u_{zz} - p_y = f^y
      9:   w_{xx} + w_{yy} + w_{zz} - p_y = f^z
     10:   u_x    + v_y    + w_z          = g

     12:   g = 0 for the physical case.

     14:   Boundary conditions give prescribed flow perpendicular to the boundaries,
     15:   and zero derivative perpendicular to them (free slip). This involves
     16:   using a modified stencil at the boundaries. Another option would be to
     17:   use DM_BOUNDARY_GHOSTED in DMStagCreate3d() and a matrix-free operator (MATSHELL)
     18:   making use of the uniformly-available ghost layer.

     20:   Use the -pinpressure option to fix a pressure node, instead of providing
     21:   a constant-pressure nullspace. This allows use of direct solvers, e.g. to
     22:   use UMFPACK,

     24:      ./ex3 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack

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

     31: /* Shorter, more convenient names for DMStagStencilLocation entries */
     32: #define BACK_DOWN_LEFT   DMSTAG_BACK_DOWN_LEFT
     33: #define BACK_DOWN        DMSTAG_BACK_DOWN
     34: #define BACK_DOWN_RIGHT  DMSTAG_BACK_DOWN_RIGHT
     35: #define BACK_LEFT        DMSTAG_BACK_LEFT
     36: #define BACK             DMSTAG_BACK
     37: #define BACK_RIGHT       DMSTAG_BACK_RIGHT
     38: #define BACK_UP_LEFT     DMSTAG_BACK_UP_LEFT
     39: #define BACK_UP          DMSTAG_BACK_UP
     40: #define BACK_UP_RIGHT    DMSTAG_BACK_UP_RIGHT
     41: #define DOWN_LEFT        DMSTAG_DOWN_LEFT
     42: #define DOWN             DMSTAG_DOWN
     43: #define DOWN_RIGHT       DMSTAG_DOWN_RIGHT
     44: #define LEFT             DMSTAG_LEFT
     45: #define ELEMENT          DMSTAG_ELEMENT
     46: #define RIGHT            DMSTAG_RIGHT
     47: #define UP_LEFT          DMSTAG_UP_LEFT
     48: #define UP               DMSTAG_UP
     49: #define UP_RIGHT         DMSTAG_UP_RIGHT
     50: #define FRONT_DOWN_LEFT  DMSTAG_FRONT_DOWN_LEFT
     51: #define FRONT_DOWN       DMSTAG_FRONT_DOWN
     52: #define FRONT_DOWN_RIGHT DMSTAG_FRONT_DOWN_RIGHT
     53: #define FRONT_LEFT       DMSTAG_FRONT_LEFT
     54: #define FRONT            DMSTAG_FRONT
     55: #define FRONT_RIGHT      DMSTAG_FRONT_RIGHT
     56: #define FRONT_UP_LEFT    DMSTAG_FRONT_UP_LEFT
     57: #define FRONT_UP         DMSTAG_FRONT_UP
     58: #define FRONT_UP_RIGHT   DMSTAG_FRONT_UP_RIGHT

     60: static PetscErrorCode CreateReferenceSolution(DM, Vec *);
     61: static PetscErrorCode CreateSystem(DM, Mat *, Vec *, PetscBool);
     62: static PetscErrorCode AttachNullspace(DM, Mat);
     63: static PetscErrorCode CheckSolution(Vec, Vec);

     65: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
     66: and to have a zero derivative for flow parallel to the boundaries. That is,
     67: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
     68: and left boundaries.
     69: These expressions could be made more interesting with more z terms,
     70: and convergence could be confirmed.  */

     72: static PetscScalar uxRef(PetscScalar x, PetscScalar y, PetscScalar z)
     73: {
     74:   return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y + 0.0 * z;
     75: }
     76: static PetscScalar uyRef(PetscScalar x, PetscScalar y, PetscScalar z)
     77: {
     78:   return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y + 0.0 * z;
     79: }
     80: static PetscScalar uzRef(PetscScalar x, PetscScalar y, PetscScalar z)
     81: {
     82:   return 0.0 * x + 0.0 * y + 0.0 * z + 1.0;
     83: }
     84: static PetscScalar pRef(PetscScalar x, PetscScalar y, PetscScalar z)
     85: {
     86:   return -1.0 * (x - 0.5) + -3.0 / 2.0 * y * y + 0.5 - 2.0 * (z - 1.0);
     87: } /* zero integral */
     88: static PetscScalar fx(PetscScalar x, PetscScalar y, PetscScalar z)
     89: {
     90:   return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 0.0 * z + 1.0;
     91: }
     92: static PetscScalar fy(PetscScalar x, PetscScalar y, PetscScalar z)
     93: {
     94:   return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y + 0.0 * z;
     95: }
     96: static PetscScalar fz(PetscScalar x, PetscScalar y, PetscScalar z)
     97: {
     98:   return 0.0 * x + 0.0 * y + 0.0 * z + 2.0;
     99: }
    100: static PetscScalar g(PetscScalar x, PetscScalar y, PetscScalar z)
    101: {
    102:   return 0.0 * x + 0.0 * y + 0.0 * z + 0.0;
    103: }

    105: int main(int argc, char **argv)
    106: {
    107:   DM        dmSol;
    108:   Vec       sol, solRef, rhs;
    109:   Mat       A;
    110:   KSP       ksp;
    111:   PC        pc;
    112:   PetscBool pinPressure;

    114:   /* Initialize PETSc and process command line arguments */
    115:   PetscFunctionBeginUser;
    116:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    117:   pinPressure = PETSC_TRUE;
    118:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pinpressure", &pinPressure, NULL));

    120:   /* Create 3D DMStag for the solution, and set up. */
    121:   {
    122:     const PetscInt dof0 = 0, dof1 = 0, dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
    123:     const PetscInt stencilWidth = 1;
    124:     PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 5, 6, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, NULL, &dmSol));
    125:     PetscCall(DMSetFromOptions(dmSol));
    126:     PetscCall(DMSetUp(dmSol));
    127:     PetscCall(DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    128:     /* Note: also see ex2.c, where another, more efficient option is demonstrated,
    129:        using DMStagSetUniformCoordinatesProduct() */
    130:   }

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

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

    138:   /* Attach a constant-pressure nullspace to the operator */
    139:   if (!pinPressure) PetscCall(AttachNullspace(dmSol, A));

    141:   /* Solve */
    142:   PetscCall(DMCreateGlobalVector(dmSol, &sol));
    143:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    144:   PetscCall(KSPSetType(ksp, KSPFGMRES));
    145:   PetscCall(KSPSetOperators(ksp, A, A));
    146:   PetscCall(KSPGetPC(ksp, &pc));
    147:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
    148:   PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE));
    149:   PetscCall(KSPSetFromOptions(ksp));
    150:   PetscCall(KSPSolve(ksp, rhs, sol));

    152:   /* Check Solution */
    153:   PetscCall(CheckSolution(sol, solRef));

    155:   /* Clean up and finalize PETSc */
    156:   PetscCall(KSPDestroy(&ksp));
    157:   PetscCall(VecDestroy(&sol));
    158:   PetscCall(VecDestroy(&solRef));
    159:   PetscCall(VecDestroy(&rhs));
    160:   PetscCall(MatDestroy(&A));
    161:   PetscCall(DMDestroy(&dmSol));
    162:   PetscCall(PetscFinalize());
    163:   return 0;
    164: }

    166: static PetscErrorCode CreateSystem(DM dmSol, Mat *pA, Vec *pRhs, PetscBool pinPressure)
    167: {
    168:   Vec             rhs, coordLocal;
    169:   Mat             A;
    170:   PetscInt        startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
    171:   PetscInt        icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
    172:   PetscReal       hx, hy, hz;
    173:   DM              dmCoord;
    174:   PetscScalar ****arrCoord;

    176:   PetscFunctionBeginUser;
    177:   PetscCall(DMCreateMatrix(dmSol, pA));
    178:   A = *pA;
    179:   PetscCall(DMCreateGlobalVector(dmSol, pRhs));
    180:   rhs = *pRhs;

    182:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
    183:   PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
    184:   PetscCheck(N[0] >= 2 && N[1] >= 2 && N[2] >= 2, PetscObjectComm((PetscObject)dmSol), PETSC_ERR_ARG_SIZ, "This example requires at least two elements in each dimensions");
    185:   hx = 1.0 / N[0];
    186:   hy = 1.0 / N[1];
    187:   hz = 1.0 / N[2];
    188:   PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
    189:   PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
    190:   PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
    191:   for (d = 0; d < 3; ++d) {
    192:     PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
    193:     PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
    194:     PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
    195:     PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
    196:     PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
    197:     PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
    198:     PetscCall(DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]));
    199:   }

    201:   /* Loop over all local elements. Note that it may be more efficient in real
    202:      applications to loop over each boundary separately */
    203:   for (ez = startz; ez < startz + nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
    204:     for (ey = starty; ey < starty + ny; ++ey) {
    205:       for (ex = startx; ex < startx + nx; ++ex) {
    206:         if (ex == N[0] - 1) {
    207:           /* Right Boundary velocity Dirichlet */
    208:           DMStagStencil     row;
    209:           PetscScalar       valRhs;
    210:           const PetscScalar valA = 1.0;
    211:           row.i                  = ex;
    212:           row.j                  = ey;
    213:           row.k                  = ez;
    214:           row.loc                = RIGHT;
    215:           row.c                  = 0;
    216:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    217:           valRhs = uxRef(arrCoord[ez][ey][ex][icux_right[0]], arrCoord[ez][ey][ex][icux_right[1]], arrCoord[ez][ey][ex][icux_right[2]]);
    218:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    219:         }
    220:         if (ey == N[1] - 1) {
    221:           /* Top boundary velocity Dirichlet */
    222:           DMStagStencil     row;
    223:           PetscScalar       valRhs;
    224:           const PetscScalar valA = 1.0;
    225:           row.i                  = ex;
    226:           row.j                  = ey;
    227:           row.k                  = ez;
    228:           row.loc                = UP;
    229:           row.c                  = 0;
    230:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    231:           valRhs = uyRef(arrCoord[ez][ey][ex][icuy_up[0]], arrCoord[ez][ey][ex][icuy_up[1]], arrCoord[ez][ey][ex][icuy_up[2]]);
    232:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    233:         }
    234:         if (ez == N[2] - 1) {
    235:           /* Front boundary velocity Dirichlet */
    236:           DMStagStencil     row;
    237:           PetscScalar       valRhs;
    238:           const PetscScalar valA = 1.0;
    239:           row.i                  = ex;
    240:           row.j                  = ey;
    241:           row.k                  = ez;
    242:           row.loc                = FRONT;
    243:           row.c                  = 0;
    244:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    245:           valRhs = uzRef(arrCoord[ez][ey][ex][icuz_front[0]], arrCoord[ez][ey][ex][icuz_front[1]], arrCoord[ez][ey][ex][icuz_front[2]]);
    246:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    247:         }

    249:         /* Equation on left face of this element */
    250:         if (ex == 0) {
    251:           /* Left velocity Dirichlet */
    252:           DMStagStencil     row;
    253:           PetscScalar       valRhs;
    254:           const PetscScalar valA = 1.0;
    255:           row.i                  = ex;
    256:           row.j                  = ey;
    257:           row.k                  = ez;
    258:           row.loc                = LEFT;
    259:           row.c                  = 0;
    260:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    261:           valRhs = uxRef(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
    262:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    263:         } else {
    264:           /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
    265:           DMStagStencil row, col[9];
    266:           PetscScalar   valA[9], valRhs;
    267:           PetscInt      nEntries;

    269:           row.i   = ex;
    270:           row.j   = ey;
    271:           row.k   = ez;
    272:           row.loc = LEFT;
    273:           row.c   = 0;
    274:           if (ey == 0) {
    275:             if (ez == 0) {
    276:               nEntries   = 7;
    277:               col[0].i   = ex;
    278:               col[0].j   = ey;
    279:               col[0].k   = ez;
    280:               col[0].loc = LEFT;
    281:               col[0].c   = 0;
    282:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
    283:               /* Missing down term */
    284:               col[1].i   = ex;
    285:               col[1].j   = ey + 1;
    286:               col[1].k   = ez;
    287:               col[1].loc = LEFT;
    288:               col[1].c   = 0;
    289:               valA[1]    = 1.0 / (hy * hy);
    290:               col[2].i   = ex - 1;
    291:               col[2].j   = ey;
    292:               col[2].k   = ez;
    293:               col[2].loc = LEFT;
    294:               col[2].c   = 0;
    295:               valA[2]    = 1.0 / (hx * hx);
    296:               col[3].i   = ex + 1;
    297:               col[3].j   = ey;
    298:               col[3].k   = ez;
    299:               col[3].loc = LEFT;
    300:               col[3].c   = 0;
    301:               valA[3]    = 1.0 / (hx * hx);
    302:               /* Missing back term */
    303:               col[4].i   = ex;
    304:               col[4].j   = ey;
    305:               col[4].k   = ez + 1;
    306:               col[4].loc = LEFT;
    307:               col[4].c   = 0;
    308:               valA[4]    = 1.0 / (hz * hz);
    309:               col[5].i   = ex - 1;
    310:               col[5].j   = ey;
    311:               col[5].k   = ez;
    312:               col[5].loc = ELEMENT;
    313:               col[5].c   = 0;
    314:               valA[5]    = 1.0 / hx;
    315:               col[6].i   = ex;
    316:               col[6].j   = ey;
    317:               col[6].k   = ez;
    318:               col[6].loc = ELEMENT;
    319:               col[6].c   = 0;
    320:               valA[6]    = -1.0 / hx;
    321:             } else if (ez == N[2] - 1) {
    322:               nEntries   = 7;
    323:               col[0].i   = ex;
    324:               col[0].j   = ey;
    325:               col[0].k   = ez;
    326:               col[0].loc = LEFT;
    327:               col[0].c   = 0;
    328:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
    329:               /* Missing down term */
    330:               col[1].i   = ex;
    331:               col[1].j   = ey + 1;
    332:               col[1].k   = ez;
    333:               col[1].loc = LEFT;
    334:               col[1].c   = 0;
    335:               valA[1]    = 1.0 / (hy * hy);
    336:               col[2].i   = ex - 1;
    337:               col[2].j   = ey;
    338:               col[2].k   = ez;
    339:               col[2].loc = LEFT;
    340:               col[2].c   = 0;
    341:               valA[2]    = 1.0 / (hx * hx);
    342:               col[3].i   = ex + 1;
    343:               col[3].j   = ey;
    344:               col[3].k   = ez;
    345:               col[3].loc = LEFT;
    346:               col[3].c   = 0;
    347:               valA[3]    = 1.0 / (hx * hx);
    348:               col[4].i   = ex;
    349:               col[4].j   = ey;
    350:               col[4].k   = ez - 1;
    351:               col[4].loc = LEFT;
    352:               col[4].c   = 0;
    353:               valA[4]    = 1.0 / (hz * hz);
    354:               /* Missing front term */
    355:               col[5].i   = ex - 1;
    356:               col[5].j   = ey;
    357:               col[5].k   = ez;
    358:               col[5].loc = ELEMENT;
    359:               col[5].c   = 0;
    360:               valA[5]    = 1.0 / hx;
    361:               col[6].i   = ex;
    362:               col[6].j   = ey;
    363:               col[6].k   = ez;
    364:               col[6].loc = ELEMENT;
    365:               col[6].c   = 0;
    366:               valA[6]    = -1.0 / hx;
    367:             } else {
    368:               nEntries   = 8;
    369:               col[0].i   = ex;
    370:               col[0].j   = ey;
    371:               col[0].k   = ez;
    372:               col[0].loc = LEFT;
    373:               col[0].c   = 0;
    374:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    375:               /* Missing down term */
    376:               col[1].i   = ex;
    377:               col[1].j   = ey + 1;
    378:               col[1].k   = ez;
    379:               col[1].loc = LEFT;
    380:               col[1].c   = 0;
    381:               valA[1]    = 1.0 / (hy * hy);
    382:               col[2].i   = ex - 1;
    383:               col[2].j   = ey;
    384:               col[2].k   = ez;
    385:               col[2].loc = LEFT;
    386:               col[2].c   = 0;
    387:               valA[2]    = 1.0 / (hx * hx);
    388:               col[3].i   = ex + 1;
    389:               col[3].j   = ey;
    390:               col[3].k   = ez;
    391:               col[3].loc = LEFT;
    392:               col[3].c   = 0;
    393:               valA[3]    = 1.0 / (hx * hx);
    394:               col[4].i   = ex;
    395:               col[4].j   = ey;
    396:               col[4].k   = ez - 1;
    397:               col[4].loc = LEFT;
    398:               col[4].c   = 0;
    399:               valA[4]    = 1.0 / (hz * hz);
    400:               col[5].i   = ex;
    401:               col[5].j   = ey;
    402:               col[5].k   = ez + 1;
    403:               col[5].loc = LEFT;
    404:               col[5].c   = 0;
    405:               valA[5]    = 1.0 / (hz * hz);
    406:               col[6].i   = ex - 1;
    407:               col[6].j   = ey;
    408:               col[6].k   = ez;
    409:               col[6].loc = ELEMENT;
    410:               col[6].c   = 0;
    411:               valA[6]    = 1.0 / hx;
    412:               col[7].i   = ex;
    413:               col[7].j   = ey;
    414:               col[7].k   = ez;
    415:               col[7].loc = ELEMENT;
    416:               col[7].c   = 0;
    417:               valA[7]    = -1.0 / hx;
    418:             }
    419:           } else if (ey == N[1] - 1) {
    420:             if (ez == 0) {
    421:               nEntries   = 7;
    422:               col[0].i   = ex;
    423:               col[0].j   = ey;
    424:               col[0].k   = ez;
    425:               col[0].loc = LEFT;
    426:               col[0].c   = 0;
    427:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    428:               col[1].i   = ex;
    429:               col[1].j   = ey - 1;
    430:               col[1].k   = ez;
    431:               col[1].loc = LEFT;
    432:               col[1].c   = 0;
    433:               valA[1]    = 1.0 / (hy * hy);
    434:               /* Missing up term */
    435:               col[2].i   = ex - 1;
    436:               col[2].j   = ey;
    437:               col[2].k   = ez;
    438:               col[2].loc = LEFT;
    439:               col[2].c   = 0;
    440:               valA[2]    = 1.0 / (hx * hx);
    441:               col[3].i   = ex + 1;
    442:               col[3].j   = ey;
    443:               col[3].k   = ez;
    444:               col[3].loc = LEFT;
    445:               col[3].c   = 0;
    446:               valA[3]    = 1.0 / (hx * hx);
    447:               /* Missing back entry */
    448:               col[4].i   = ex;
    449:               col[4].j   = ey;
    450:               col[4].k   = ez + 1;
    451:               col[4].loc = LEFT;
    452:               col[4].c   = 0;
    453:               valA[4]    = 1.0 / (hz * hz);
    454:               col[5].i   = ex - 1;
    455:               col[5].j   = ey;
    456:               col[5].k   = ez;
    457:               col[5].loc = ELEMENT;
    458:               col[5].c   = 0;
    459:               valA[5]    = 1.0 / hx;
    460:               col[6].i   = ex;
    461:               col[6].j   = ey;
    462:               col[6].k   = ez;
    463:               col[6].loc = ELEMENT;
    464:               col[6].c   = 0;
    465:               valA[6]    = -1.0 / hx;
    466:             } else if (ez == N[2] - 1) {
    467:               nEntries   = 7;
    468:               col[0].i   = ex;
    469:               col[0].j   = ey;
    470:               col[0].k   = ez;
    471:               col[0].loc = LEFT;
    472:               col[0].c   = 0;
    473:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    474:               col[1].i   = ex;
    475:               col[1].j   = ey - 1;
    476:               col[1].k   = ez;
    477:               col[1].loc = LEFT;
    478:               col[1].c   = 0;
    479:               valA[1]    = 1.0 / (hy * hy);
    480:               /* Missing up term */
    481:               col[2].i   = ex - 1;
    482:               col[2].j   = ey;
    483:               col[2].k   = ez;
    484:               col[2].loc = LEFT;
    485:               col[2].c   = 0;
    486:               valA[2]    = 1.0 / (hx * hx);
    487:               col[3].i   = ex + 1;
    488:               col[3].j   = ey;
    489:               col[3].k   = ez;
    490:               col[3].loc = LEFT;
    491:               col[3].c   = 0;
    492:               valA[3]    = 1.0 / (hx * hx);
    493:               col[4].i   = ex;
    494:               col[4].j   = ey;
    495:               col[4].k   = ez - 1;
    496:               col[4].loc = LEFT;
    497:               col[4].c   = 0;
    498:               valA[4]    = 1.0 / (hz * hz);
    499:               /* Missing front term */
    500:               col[5].i   = ex - 1;
    501:               col[5].j   = ey;
    502:               col[5].k   = ez;
    503:               col[5].loc = ELEMENT;
    504:               col[5].c   = 0;
    505:               valA[5]    = 1.0 / hx;
    506:               col[6].i   = ex;
    507:               col[6].j   = ey;
    508:               col[6].k   = ez;
    509:               col[6].loc = ELEMENT;
    510:               col[6].c   = 0;
    511:               valA[6]    = -1.0 / hx;
    512:             } else {
    513:               nEntries   = 8;
    514:               col[0].i   = ex;
    515:               col[0].j   = ey;
    516:               col[0].k   = ez;
    517:               col[0].loc = LEFT;
    518:               col[0].c   = 0;
    519:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    520:               col[1].i   = ex;
    521:               col[1].j   = ey - 1;
    522:               col[1].k   = ez;
    523:               col[1].loc = LEFT;
    524:               col[1].c   = 0;
    525:               valA[1]    = 1.0 / (hy * hy);
    526:               /* Missing up term */
    527:               col[2].i   = ex - 1;
    528:               col[2].j   = ey;
    529:               col[2].k   = ez;
    530:               col[2].loc = LEFT;
    531:               col[2].c   = 0;
    532:               valA[2]    = 1.0 / (hx * hx);
    533:               col[3].i   = ex + 1;
    534:               col[3].j   = ey;
    535:               col[3].k   = ez;
    536:               col[3].loc = LEFT;
    537:               col[3].c   = 0;
    538:               valA[3]    = 1.0 / (hx * hx);
    539:               col[4].i   = ex;
    540:               col[4].j   = ey;
    541:               col[4].k   = ez - 1;
    542:               col[4].loc = LEFT;
    543:               col[4].c   = 0;
    544:               valA[4]    = 1.0 / (hz * hz);
    545:               col[5].i   = ex;
    546:               col[5].j   = ey;
    547:               col[5].k   = ez + 1;
    548:               col[5].loc = LEFT;
    549:               col[5].c   = 0;
    550:               valA[5]    = 1.0 / (hz * hz);
    551:               col[6].i   = ex - 1;
    552:               col[6].j   = ey;
    553:               col[6].k   = ez;
    554:               col[6].loc = ELEMENT;
    555:               col[6].c   = 0;
    556:               valA[6]    = 1.0 / hx;
    557:               col[7].i   = ex;
    558:               col[7].j   = ey;
    559:               col[7].k   = ez;
    560:               col[7].loc = ELEMENT;
    561:               col[7].c   = 0;
    562:               valA[7]    = -1.0 / hx;
    563:             }
    564:           } else if (ez == 0) {
    565:             nEntries   = 8;
    566:             col[0].i   = ex;
    567:             col[0].j   = ey;
    568:             col[0].k   = ez;
    569:             col[0].loc = LEFT;
    570:             col[0].c   = 0;
    571:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    572:             col[1].i   = ex;
    573:             col[1].j   = ey - 1;
    574:             col[1].k   = ez;
    575:             col[1].loc = LEFT;
    576:             col[1].c   = 0;
    577:             valA[1]    = 1.0 / (hy * hy);
    578:             col[2].i   = ex;
    579:             col[2].j   = ey + 1;
    580:             col[2].k   = ez;
    581:             col[2].loc = LEFT;
    582:             col[2].c   = 0;
    583:             valA[2]    = 1.0 / (hy * hy);
    584:             col[3].i   = ex - 1;
    585:             col[3].j   = ey;
    586:             col[3].k   = ez;
    587:             col[3].loc = LEFT;
    588:             col[3].c   = 0;
    589:             valA[3]    = 1.0 / (hx * hx);
    590:             col[4].i   = ex + 1;
    591:             col[4].j   = ey;
    592:             col[4].k   = ez;
    593:             col[4].loc = LEFT;
    594:             col[4].c   = 0;
    595:             valA[4]    = 1.0 / (hx * hx);
    596:             /* Missing back term */
    597:             col[5].i   = ex;
    598:             col[5].j   = ey;
    599:             col[5].k   = ez + 1;
    600:             col[5].loc = LEFT;
    601:             col[5].c   = 0;
    602:             valA[5]    = 1.0 / (hz * hz);
    603:             col[6].i   = ex - 1;
    604:             col[6].j   = ey;
    605:             col[6].k   = ez;
    606:             col[6].loc = ELEMENT;
    607:             col[6].c   = 0;
    608:             valA[6]    = 1.0 / hx;
    609:             col[7].i   = ex;
    610:             col[7].j   = ey;
    611:             col[7].k   = ez;
    612:             col[7].loc = ELEMENT;
    613:             col[7].c   = 0;
    614:             valA[7]    = -1.0 / hx;
    615:           } else if (ez == N[2] - 1) {
    616:             nEntries   = 8;
    617:             col[0].i   = ex;
    618:             col[0].j   = ey;
    619:             col[0].k   = ez;
    620:             col[0].loc = LEFT;
    621:             col[0].c   = 0;
    622:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    623:             col[1].i   = ex;
    624:             col[1].j   = ey - 1;
    625:             col[1].k   = ez;
    626:             col[1].loc = LEFT;
    627:             col[1].c   = 0;
    628:             valA[1]    = 1.0 / (hy * hy);
    629:             col[2].i   = ex;
    630:             col[2].j   = ey + 1;
    631:             col[2].k   = ez;
    632:             col[2].loc = LEFT;
    633:             col[2].c   = 0;
    634:             valA[2]    = 1.0 / (hy * hy);
    635:             col[3].i   = ex - 1;
    636:             col[3].j   = ey;
    637:             col[3].k   = ez;
    638:             col[3].loc = LEFT;
    639:             col[3].c   = 0;
    640:             valA[3]    = 1.0 / (hx * hx);
    641:             col[4].i   = ex + 1;
    642:             col[4].j   = ey;
    643:             col[4].k   = ez;
    644:             col[4].loc = LEFT;
    645:             col[4].c   = 0;
    646:             valA[4]    = 1.0 / (hx * hx);
    647:             col[5].i   = ex;
    648:             col[5].j   = ey;
    649:             col[5].k   = ez - 1;
    650:             col[5].loc = LEFT;
    651:             col[5].c   = 0;
    652:             valA[5]    = 1.0 / (hz * hz);
    653:             /* Missing front term */
    654:             col[6].i   = ex - 1;
    655:             col[6].j   = ey;
    656:             col[6].k   = ez;
    657:             col[6].loc = ELEMENT;
    658:             col[6].c   = 0;
    659:             valA[6]    = 1.0 / hx;
    660:             col[7].i   = ex;
    661:             col[7].j   = ey;
    662:             col[7].k   = ez;
    663:             col[7].loc = ELEMENT;
    664:             col[7].c   = 0;
    665:             valA[7]    = -1.0 / hx;
    666:           } else {
    667:             nEntries   = 9;
    668:             col[0].i   = ex;
    669:             col[0].j   = ey;
    670:             col[0].k   = ez;
    671:             col[0].loc = LEFT;
    672:             col[0].c   = 0;
    673:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    674:             col[1].i   = ex;
    675:             col[1].j   = ey - 1;
    676:             col[1].k   = ez;
    677:             col[1].loc = LEFT;
    678:             col[1].c   = 0;
    679:             valA[1]    = 1.0 / (hy * hy);
    680:             col[2].i   = ex;
    681:             col[2].j   = ey + 1;
    682:             col[2].k   = ez;
    683:             col[2].loc = LEFT;
    684:             col[2].c   = 0;
    685:             valA[2]    = 1.0 / (hy * hy);
    686:             col[3].i   = ex - 1;
    687:             col[3].j   = ey;
    688:             col[3].k   = ez;
    689:             col[3].loc = LEFT;
    690:             col[3].c   = 0;
    691:             valA[3]    = 1.0 / (hx * hx);
    692:             col[4].i   = ex + 1;
    693:             col[4].j   = ey;
    694:             col[4].k   = ez;
    695:             col[4].loc = LEFT;
    696:             col[4].c   = 0;
    697:             valA[4]    = 1.0 / (hx * hx);
    698:             col[5].i   = ex;
    699:             col[5].j   = ey;
    700:             col[5].k   = ez - 1;
    701:             col[5].loc = LEFT;
    702:             col[5].c   = 0;
    703:             valA[5]    = 1.0 / (hz * hz);
    704:             col[6].i   = ex;
    705:             col[6].j   = ey;
    706:             col[6].k   = ez + 1;
    707:             col[6].loc = LEFT;
    708:             col[6].c   = 0;
    709:             valA[6]    = 1.0 / (hz * hz);
    710:             col[7].i   = ex - 1;
    711:             col[7].j   = ey;
    712:             col[7].k   = ez;
    713:             col[7].loc = ELEMENT;
    714:             col[7].c   = 0;
    715:             valA[7]    = 1.0 / hx;
    716:             col[8].i   = ex;
    717:             col[8].j   = ey;
    718:             col[8].k   = ez;
    719:             col[8].loc = ELEMENT;
    720:             col[8].c   = 0;
    721:             valA[8]    = -1.0 / hx;
    722:           }
    723:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
    724:           valRhs = fx(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
    725:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    726:         }

    728:         /* Equation on bottom face of this element */
    729:         if (ey == 0) {
    730:           /* Bottom boundary velocity Dirichlet */
    731:           DMStagStencil     row;
    732:           PetscScalar       valRhs;
    733:           const PetscScalar valA = 1.0;
    734:           row.i                  = ex;
    735:           row.j                  = ey;
    736:           row.k                  = ez;
    737:           row.loc                = DOWN;
    738:           row.c                  = 0;
    739:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    740:           valRhs = uyRef(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
    741:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    742:         } else {
    743:           /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
    744:           DMStagStencil row, col[9];
    745:           PetscScalar   valA[9], valRhs;
    746:           PetscInt      nEntries;

    748:           row.i   = ex;
    749:           row.j   = ey;
    750:           row.k   = ez;
    751:           row.loc = DOWN;
    752:           row.c   = 0;
    753:           if (ex == 0) {
    754:             if (ez == 0) {
    755:               nEntries   = 7;
    756:               col[0].i   = ex;
    757:               col[0].j   = ey;
    758:               col[0].k   = ez;
    759:               col[0].loc = DOWN;
    760:               col[0].c   = 0;
    761:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    762:               col[1].i   = ex;
    763:               col[1].j   = ey - 1;
    764:               col[1].k   = ez;
    765:               col[1].loc = DOWN;
    766:               col[1].c   = 0;
    767:               valA[1]    = 1.0 / (hy * hy);
    768:               col[2].i   = ex;
    769:               col[2].j   = ey + 1;
    770:               col[2].k   = ez;
    771:               col[2].loc = DOWN;
    772:               col[2].c   = 0;
    773:               valA[2]    = 1.0 / (hy * hy);
    774:               /* Left term missing */
    775:               col[3].i   = ex + 1;
    776:               col[3].j   = ey;
    777:               col[3].k   = ez;
    778:               col[3].loc = DOWN;
    779:               col[3].c   = 0;
    780:               valA[3]    = 1.0 / (hx * hx);
    781:               /* Back term missing */
    782:               col[4].i   = ex;
    783:               col[4].j   = ey;
    784:               col[4].k   = ez + 1;
    785:               col[4].loc = DOWN;
    786:               col[4].c   = 0;
    787:               valA[4]    = 1.0 / (hz * hz);
    788:               col[5].i   = ex;
    789:               col[5].j   = ey - 1;
    790:               col[5].k   = ez;
    791:               col[5].loc = ELEMENT;
    792:               col[5].c   = 0;
    793:               valA[5]    = 1.0 / hy;
    794:               col[6].i   = ex;
    795:               col[6].j   = ey;
    796:               col[6].k   = ez;
    797:               col[6].loc = ELEMENT;
    798:               col[6].c   = 0;
    799:               valA[6]    = -1.0 / hy;
    800:             } else if (ez == N[2] - 1) {
    801:               nEntries   = 7;
    802:               col[0].i   = ex;
    803:               col[0].j   = ey;
    804:               col[0].k   = ez;
    805:               col[0].loc = DOWN;
    806:               col[0].c   = 0;
    807:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    808:               col[1].i   = ex;
    809:               col[1].j   = ey - 1;
    810:               col[1].k   = ez;
    811:               col[1].loc = DOWN;
    812:               col[1].c   = 0;
    813:               valA[1]    = 1.0 / (hy * hy);
    814:               col[2].i   = ex;
    815:               col[2].j   = ey + 1;
    816:               col[2].k   = ez;
    817:               col[2].loc = DOWN;
    818:               col[2].c   = 0;
    819:               valA[2]    = 1.0 / (hy * hy);
    820:               /* Left term missing */
    821:               col[3].i   = ex + 1;
    822:               col[3].j   = ey;
    823:               col[3].k   = ez;
    824:               col[3].loc = DOWN;
    825:               col[3].c   = 0;
    826:               valA[3]    = 1.0 / (hx * hx);
    827:               col[4].i   = ex;
    828:               col[4].j   = ey;
    829:               col[4].k   = ez - 1;
    830:               col[4].loc = DOWN;
    831:               col[4].c   = 0;
    832:               valA[4]    = 1.0 / (hz * hz);
    833:               /* Front term missing */
    834:               col[5].i   = ex;
    835:               col[5].j   = ey - 1;
    836:               col[5].k   = ez;
    837:               col[5].loc = ELEMENT;
    838:               col[5].c   = 0;
    839:               valA[5]    = 1.0 / hy;
    840:               col[6].i   = ex;
    841:               col[6].j   = ey;
    842:               col[6].k   = ez;
    843:               col[6].loc = ELEMENT;
    844:               col[6].c   = 0;
    845:               valA[6]    = -1.0 / hy;
    846:             } else {
    847:               nEntries   = 8;
    848:               col[0].i   = ex;
    849:               col[0].j   = ey;
    850:               col[0].k   = ez;
    851:               col[0].loc = DOWN;
    852:               col[0].c   = 0;
    853:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    854:               col[1].i   = ex;
    855:               col[1].j   = ey - 1;
    856:               col[1].k   = ez;
    857:               col[1].loc = DOWN;
    858:               col[1].c   = 0;
    859:               valA[1]    = 1.0 / (hy * hy);
    860:               col[2].i   = ex;
    861:               col[2].j   = ey + 1;
    862:               col[2].k   = ez;
    863:               col[2].loc = DOWN;
    864:               col[2].c   = 0;
    865:               valA[2]    = 1.0 / (hy * hy);
    866:               /* Left term missing */
    867:               col[3].i   = ex + 1;
    868:               col[3].j   = ey;
    869:               col[3].k   = ez;
    870:               col[3].loc = DOWN;
    871:               col[3].c   = 0;
    872:               valA[3]    = 1.0 / (hx * hx);
    873:               col[4].i   = ex;
    874:               col[4].j   = ey;
    875:               col[4].k   = ez - 1;
    876:               col[4].loc = DOWN;
    877:               col[4].c   = 0;
    878:               valA[4]    = 1.0 / (hz * hz);
    879:               col[5].i   = ex;
    880:               col[5].j   = ey;
    881:               col[5].k   = ez + 1;
    882:               col[5].loc = DOWN;
    883:               col[5].c   = 0;
    884:               valA[5]    = 1.0 / (hz * hz);
    885:               col[6].i   = ex;
    886:               col[6].j   = ey - 1;
    887:               col[6].k   = ez;
    888:               col[6].loc = ELEMENT;
    889:               col[6].c   = 0;
    890:               valA[6]    = 1.0 / hy;
    891:               col[7].i   = ex;
    892:               col[7].j   = ey;
    893:               col[7].k   = ez;
    894:               col[7].loc = ELEMENT;
    895:               col[7].c   = 0;
    896:               valA[7]    = -1.0 / hy;
    897:             }
    898:           } else if (ex == N[0] - 1) {
    899:             if (ez == 0) {
    900:               nEntries   = 7;
    901:               col[0].i   = ex;
    902:               col[0].j   = ey;
    903:               col[0].k   = ez;
    904:               col[0].loc = DOWN;
    905:               col[0].c   = 0;
    906:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    907:               col[1].i   = ex;
    908:               col[1].j   = ey - 1;
    909:               col[1].k   = ez;
    910:               col[1].loc = DOWN;
    911:               col[1].c   = 0;
    912:               valA[1]    = 1.0 / (hy * hy);
    913:               col[2].i   = ex;
    914:               col[2].j   = ey + 1;
    915:               col[2].k   = ez;
    916:               col[2].loc = DOWN;
    917:               col[2].c   = 0;
    918:               valA[2]    = 1.0 / (hy * hy);
    919:               col[3].i   = ex - 1;
    920:               col[3].j   = ey;
    921:               col[3].k   = ez;
    922:               col[3].loc = DOWN;
    923:               col[3].c   = 0;
    924:               valA[3]    = 1.0 / (hx * hx);
    925:               /* Right term missing */
    926:               /* Back term missing */
    927:               col[4].i   = ex;
    928:               col[4].j   = ey;
    929:               col[4].k   = ez + 1;
    930:               col[4].loc = DOWN;
    931:               col[4].c   = 0;
    932:               valA[4]    = 1.0 / (hz * hz);
    933:               col[5].i   = ex;
    934:               col[5].j   = ey - 1;
    935:               col[5].k   = ez;
    936:               col[5].loc = ELEMENT;
    937:               col[5].c   = 0;
    938:               valA[5]    = 1.0 / hy;
    939:               col[6].i   = ex;
    940:               col[6].j   = ey;
    941:               col[6].k   = ez;
    942:               col[6].loc = ELEMENT;
    943:               col[6].c   = 0;
    944:               valA[6]    = -1.0 / hy;
    945:             } else if (ez == N[2] - 1) {
    946:               nEntries   = 7;
    947:               col[0].i   = ex;
    948:               col[0].j   = ey;
    949:               col[0].k   = ez;
    950:               col[0].loc = DOWN;
    951:               col[0].c   = 0;
    952:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    953:               col[1].i   = ex;
    954:               col[1].j   = ey - 1;
    955:               col[1].k   = ez;
    956:               col[1].loc = DOWN;
    957:               col[1].c   = 0;
    958:               valA[1]    = 1.0 / (hy * hy);
    959:               col[2].i   = ex;
    960:               col[2].j   = ey + 1;
    961:               col[2].k   = ez;
    962:               col[2].loc = DOWN;
    963:               col[2].c   = 0;
    964:               valA[2]    = 1.0 / (hy * hy);
    965:               col[3].i   = ex - 1;
    966:               col[3].j   = ey;
    967:               col[3].k   = ez;
    968:               col[3].loc = DOWN;
    969:               col[3].c   = 0;
    970:               valA[3]    = 1.0 / (hx * hx);
    971:               /* Right term missing */
    972:               col[4].i   = ex;
    973:               col[4].j   = ey;
    974:               col[4].k   = ez - 1;
    975:               col[4].loc = DOWN;
    976:               col[4].c   = 0;
    977:               valA[4]    = 1.0 / (hz * hz);
    978:               /* Front term missing */
    979:               col[5].i   = ex;
    980:               col[5].j   = ey - 1;
    981:               col[5].k   = ez;
    982:               col[5].loc = ELEMENT;
    983:               col[5].c   = 0;
    984:               valA[5]    = 1.0 / hy;
    985:               col[6].i   = ex;
    986:               col[6].j   = ey;
    987:               col[6].k   = ez;
    988:               col[6].loc = ELEMENT;
    989:               col[6].c   = 0;
    990:               valA[6]    = -1.0 / hy;
    991:             } else {
    992:               nEntries   = 8;
    993:               col[0].i   = ex;
    994:               col[0].j   = ey;
    995:               col[0].k   = ez;
    996:               col[0].loc = DOWN;
    997:               col[0].c   = 0;
    998:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    999:               col[1].i   = ex;
    1000:               col[1].j   = ey - 1;
    1001:               col[1].k   = ez;
    1002:               col[1].loc = DOWN;
    1003:               col[1].c   = 0;
    1004:               valA[1]    = 1.0 / (hy * hy);
    1005:               col[2].i   = ex;
    1006:               col[2].j   = ey + 1;
    1007:               col[2].k   = ez;
    1008:               col[2].loc = DOWN;
    1009:               col[2].c   = 0;
    1010:               valA[2]    = 1.0 / (hy * hy);
    1011:               col[3].i   = ex - 1;
    1012:               col[3].j   = ey;
    1013:               col[3].k   = ez;
    1014:               col[3].loc = DOWN;
    1015:               col[3].c   = 0;
    1016:               valA[3]    = 1.0 / (hx * hx);
    1017:               /* Right term missing */
    1018:               col[4].i   = ex;
    1019:               col[4].j   = ey;
    1020:               col[4].k   = ez - 1;
    1021:               col[4].loc = DOWN;
    1022:               col[4].c   = 0;
    1023:               valA[4]    = 1.0 / (hz * hz);
    1024:               col[5].i   = ex;
    1025:               col[5].j   = ey;
    1026:               col[5].k   = ez + 1;
    1027:               col[5].loc = DOWN;
    1028:               col[5].c   = 0;
    1029:               valA[5]    = 1.0 / (hz * hz);
    1030:               col[6].i   = ex;
    1031:               col[6].j   = ey - 1;
    1032:               col[6].k   = ez;
    1033:               col[6].loc = ELEMENT;
    1034:               col[6].c   = 0;
    1035:               valA[6]    = 1.0 / hy;
    1036:               col[7].i   = ex;
    1037:               col[7].j   = ey;
    1038:               col[7].k   = ez;
    1039:               col[7].loc = ELEMENT;
    1040:               col[7].c   = 0;
    1041:               valA[7]    = -1.0 / hy;
    1042:             }
    1043:           } else if (ez == 0) {
    1044:             nEntries   = 8;
    1045:             col[0].i   = ex;
    1046:             col[0].j   = ey;
    1047:             col[0].k   = ez;
    1048:             col[0].loc = DOWN;
    1049:             col[0].c   = 0;
    1050:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    1051:             col[1].i   = ex;
    1052:             col[1].j   = ey - 1;
    1053:             col[1].k   = ez;
    1054:             col[1].loc = DOWN;
    1055:             col[1].c   = 0;
    1056:             valA[1]    = 1.0 / (hy * hy);
    1057:             col[2].i   = ex;
    1058:             col[2].j   = ey + 1;
    1059:             col[2].k   = ez;
    1060:             col[2].loc = DOWN;
    1061:             col[2].c   = 0;
    1062:             valA[2]    = 1.0 / (hy * hy);
    1063:             col[3].i   = ex - 1;
    1064:             col[3].j   = ey;
    1065:             col[3].k   = ez;
    1066:             col[3].loc = DOWN;
    1067:             col[3].c   = 0;
    1068:             valA[3]    = 1.0 / (hx * hx);
    1069:             col[4].i   = ex + 1;
    1070:             col[4].j   = ey;
    1071:             col[4].k   = ez;
    1072:             col[4].loc = DOWN;
    1073:             col[4].c   = 0;
    1074:             valA[4]    = 1.0 / (hx * hx);
    1075:             /* Back term missing */
    1076:             col[5].i   = ex;
    1077:             col[5].j   = ey;
    1078:             col[5].k   = ez + 1;
    1079:             col[5].loc = DOWN;
    1080:             col[5].c   = 0;
    1081:             valA[5]    = 1.0 / (hz * hz);
    1082:             col[6].i   = ex;
    1083:             col[6].j   = ey - 1;
    1084:             col[6].k   = ez;
    1085:             col[6].loc = ELEMENT;
    1086:             col[6].c   = 0;
    1087:             valA[6]    = 1.0 / hy;
    1088:             col[7].i   = ex;
    1089:             col[7].j   = ey;
    1090:             col[7].k   = ez;
    1091:             col[7].loc = ELEMENT;
    1092:             col[7].c   = 0;
    1093:             valA[7]    = -1.0 / hy;
    1094:           } else if (ez == N[2] - 1) {
    1095:             nEntries   = 8;
    1096:             col[0].i   = ex;
    1097:             col[0].j   = ey;
    1098:             col[0].k   = ez;
    1099:             col[0].loc = DOWN;
    1100:             col[0].c   = 0;
    1101:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
    1102:             col[1].i   = ex;
    1103:             col[1].j   = ey - 1;
    1104:             col[1].k   = ez;
    1105:             col[1].loc = DOWN;
    1106:             col[1].c   = 0;
    1107:             valA[1]    = 1.0 / (hy * hy);
    1108:             col[2].i   = ex;
    1109:             col[2].j   = ey + 1;
    1110:             col[2].k   = ez;
    1111:             col[2].loc = DOWN;
    1112:             col[2].c   = 0;
    1113:             valA[2]    = 1.0 / (hy * hy);
    1114:             col[3].i   = ex - 1;
    1115:             col[3].j   = ey;
    1116:             col[3].k   = ez;
    1117:             col[3].loc = DOWN;
    1118:             col[3].c   = 0;
    1119:             valA[3]    = 1.0 / (hx * hx);
    1120:             col[4].i   = ex + 1;
    1121:             col[4].j   = ey;
    1122:             col[4].k   = ez;
    1123:             col[4].loc = DOWN;
    1124:             col[4].c   = 0;
    1125:             valA[4]    = 1.0 / (hx * hx);
    1126:             col[5].i   = ex;
    1127:             col[5].j   = ey;
    1128:             col[5].k   = ez - 1;
    1129:             col[5].loc = DOWN;
    1130:             col[5].c   = 0;
    1131:             valA[5]    = 1.0 / (hz * hz);
    1132:             /* Front term missing */
    1133:             col[6].i   = ex;
    1134:             col[6].j   = ey - 1;
    1135:             col[6].k   = ez;
    1136:             col[6].loc = ELEMENT;
    1137:             col[6].c   = 0;
    1138:             valA[6]    = 1.0 / hy;
    1139:             col[7].i   = ex;
    1140:             col[7].j   = ey;
    1141:             col[7].k   = ez;
    1142:             col[7].loc = ELEMENT;
    1143:             col[7].c   = 0;
    1144:             valA[7]    = -1.0 / hy;
    1145:           } else {
    1146:             nEntries   = 9;
    1147:             col[0].i   = ex;
    1148:             col[0].j   = ey;
    1149:             col[0].k   = ez;
    1150:             col[0].loc = DOWN;
    1151:             col[0].c   = 0;
    1152:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    1153:             col[1].i   = ex;
    1154:             col[1].j   = ey - 1;
    1155:             col[1].k   = ez;
    1156:             col[1].loc = DOWN;
    1157:             col[1].c   = 0;
    1158:             valA[1]    = 1.0 / (hy * hy);
    1159:             col[2].i   = ex;
    1160:             col[2].j   = ey + 1;
    1161:             col[2].k   = ez;
    1162:             col[2].loc = DOWN;
    1163:             col[2].c   = 0;
    1164:             valA[2]    = 1.0 / (hy * hy);
    1165:             col[3].i   = ex - 1;
    1166:             col[3].j   = ey;
    1167:             col[3].k   = ez;
    1168:             col[3].loc = DOWN;
    1169:             col[3].c   = 0;
    1170:             valA[3]    = 1.0 / (hx * hx);
    1171:             col[4].i   = ex + 1;
    1172:             col[4].j   = ey;
    1173:             col[4].k   = ez;
    1174:             col[4].loc = DOWN;
    1175:             col[4].c   = 0;
    1176:             valA[4]    = 1.0 / (hx * hx);
    1177:             col[5].i   = ex;
    1178:             col[5].j   = ey;
    1179:             col[5].k   = ez - 1;
    1180:             col[5].loc = DOWN;
    1181:             col[5].c   = 0;
    1182:             valA[5]    = 1.0 / (hz * hz);
    1183:             col[6].i   = ex;
    1184:             col[6].j   = ey;
    1185:             col[6].k   = ez + 1;
    1186:             col[6].loc = DOWN;
    1187:             col[6].c   = 0;
    1188:             valA[6]    = 1.0 / (hz * hz);
    1189:             col[7].i   = ex;
    1190:             col[7].j   = ey - 1;
    1191:             col[7].k   = ez;
    1192:             col[7].loc = ELEMENT;
    1193:             col[7].c   = 0;
    1194:             valA[7]    = 1.0 / hy;
    1195:             col[8].i   = ex;
    1196:             col[8].j   = ey;
    1197:             col[8].k   = ez;
    1198:             col[8].loc = ELEMENT;
    1199:             col[8].c   = 0;
    1200:             valA[8]    = -1.0 / hy;
    1201:           }
    1202:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
    1203:           valRhs = fy(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
    1204:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    1205:         }

    1207:         /* Equation on back face of this element */
    1208:         if (ez == 0) {
    1209:           /* Back boundary velocity Dirichlet */
    1210:           DMStagStencil     row;
    1211:           PetscScalar       valRhs;
    1212:           const PetscScalar valA = 1.0;
    1213:           row.i                  = ex;
    1214:           row.j                  = ey;
    1215:           row.k                  = ez;
    1216:           row.loc                = BACK;
    1217:           row.c                  = 0;
    1218:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    1219:           valRhs = uzRef(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
    1220:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    1221:         } else {
    1222:           /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
    1223:           DMStagStencil row, col[9];
    1224:           PetscScalar   valA[9], valRhs;
    1225:           PetscInt      nEntries;

    1227:           row.i   = ex;
    1228:           row.j   = ey;
    1229:           row.k   = ez;
    1230:           row.loc = BACK;
    1231:           row.c   = 0;
    1232:           if (ex == 0) {
    1233:             if (ey == 0) {
    1234:               nEntries   = 7;
    1235:               col[0].i   = ex;
    1236:               col[0].j   = ey;
    1237:               col[0].k   = ez;
    1238:               col[0].loc = BACK;
    1239:               col[0].c   = 0;
    1240:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1241:               /* Down term missing */
    1242:               col[1].i   = ex;
    1243:               col[1].j   = ey + 1;
    1244:               col[1].k   = ez;
    1245:               col[1].loc = BACK;
    1246:               col[1].c   = 0;
    1247:               valA[1]    = 1.0 / (hy * hy);
    1248:               /* Left term missing */
    1249:               col[2].i   = ex + 1;
    1250:               col[2].j   = ey;
    1251:               col[2].k   = ez;
    1252:               col[2].loc = BACK;
    1253:               col[2].c   = 0;
    1254:               valA[2]    = 1.0 / (hx * hx);
    1255:               col[3].i   = ex;
    1256:               col[3].j   = ey;
    1257:               col[3].k   = ez - 1;
    1258:               col[3].loc = BACK;
    1259:               col[3].c   = 0;
    1260:               valA[3]    = 1.0 / (hz * hz);
    1261:               col[4].i   = ex;
    1262:               col[4].j   = ey;
    1263:               col[4].k   = ez + 1;
    1264:               col[4].loc = BACK;
    1265:               col[4].c   = 0;
    1266:               valA[4]    = 1.0 / (hz * hz);
    1267:               col[5].i   = ex;
    1268:               col[5].j   = ey;
    1269:               col[5].k   = ez - 1;
    1270:               col[5].loc = ELEMENT;
    1271:               col[5].c   = 0;
    1272:               valA[5]    = 1.0 / hz;
    1273:               col[6].i   = ex;
    1274:               col[6].j   = ey;
    1275:               col[6].k   = ez;
    1276:               col[6].loc = ELEMENT;
    1277:               col[6].c   = 0;
    1278:               valA[6]    = -1.0 / hz;
    1279:             } else if (ey == N[1] - 1) {
    1280:               nEntries   = 7;
    1281:               col[0].i   = ex;
    1282:               col[0].j   = ey;
    1283:               col[0].k   = ez;
    1284:               col[0].loc = BACK;
    1285:               col[0].c   = 0;
    1286:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1287:               col[1].i   = ex;
    1288:               col[1].j   = ey - 1;
    1289:               col[1].k   = ez;
    1290:               col[1].loc = BACK;
    1291:               col[1].c   = 0;
    1292:               valA[1]    = 1.0 / (hy * hy);
    1293:               /* Up term missing */
    1294:               /* Left term missing */
    1295:               col[2].i   = ex + 1;
    1296:               col[2].j   = ey;
    1297:               col[2].k   = ez;
    1298:               col[2].loc = BACK;
    1299:               col[2].c   = 0;
    1300:               valA[2]    = 1.0 / (hx * hx);
    1301:               col[3].i   = ex;
    1302:               col[3].j   = ey;
    1303:               col[3].k   = ez - 1;
    1304:               col[3].loc = BACK;
    1305:               col[3].c   = 0;
    1306:               valA[3]    = 1.0 / (hz * hz);
    1307:               col[4].i   = ex;
    1308:               col[4].j   = ey;
    1309:               col[4].k   = ez + 1;
    1310:               col[4].loc = BACK;
    1311:               col[4].c   = 0;
    1312:               valA[4]    = 1.0 / (hz * hz);
    1313:               col[5].i   = ex;
    1314:               col[5].j   = ey;
    1315:               col[5].k   = ez - 1;
    1316:               col[5].loc = ELEMENT;
    1317:               col[5].c   = 0;
    1318:               valA[5]    = 1.0 / hz;
    1319:               col[6].i   = ex;
    1320:               col[6].j   = ey;
    1321:               col[6].k   = ez;
    1322:               col[6].loc = ELEMENT;
    1323:               col[6].c   = 0;
    1324:               valA[6]    = -1.0 / hz;
    1325:             } else {
    1326:               nEntries   = 8;
    1327:               col[0].i   = ex;
    1328:               col[0].j   = ey;
    1329:               col[0].k   = ez;
    1330:               col[0].loc = BACK;
    1331:               col[0].c   = 0;
    1332:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    1333:               col[1].i   = ex;
    1334:               col[1].j   = ey - 1;
    1335:               col[1].k   = ez;
    1336:               col[1].loc = BACK;
    1337:               col[1].c   = 0;
    1338:               valA[1]    = 1.0 / (hy * hy);
    1339:               col[2].i   = ex;
    1340:               col[2].j   = ey + 1;
    1341:               col[2].k   = ez;
    1342:               col[2].loc = BACK;
    1343:               col[2].c   = 0;
    1344:               valA[2]    = 1.0 / (hy * hy);
    1345:               /* Left term missing */
    1346:               col[3].i   = ex + 1;
    1347:               col[3].j   = ey;
    1348:               col[3].k   = ez;
    1349:               col[3].loc = BACK;
    1350:               col[3].c   = 0;
    1351:               valA[3]    = 1.0 / (hx * hx);
    1352:               col[4].i   = ex;
    1353:               col[4].j   = ey;
    1354:               col[4].k   = ez - 1;
    1355:               col[4].loc = BACK;
    1356:               col[4].c   = 0;
    1357:               valA[4]    = 1.0 / (hz * hz);
    1358:               col[5].i   = ex;
    1359:               col[5].j   = ey;
    1360:               col[5].k   = ez + 1;
    1361:               col[5].loc = BACK;
    1362:               col[5].c   = 0;
    1363:               valA[5]    = 1.0 / (hz * hz);
    1364:               col[6].i   = ex;
    1365:               col[6].j   = ey;
    1366:               col[6].k   = ez - 1;
    1367:               col[6].loc = ELEMENT;
    1368:               col[6].c   = 0;
    1369:               valA[6]    = 1.0 / hz;
    1370:               col[7].i   = ex;
    1371:               col[7].j   = ey;
    1372:               col[7].k   = ez;
    1373:               col[7].loc = ELEMENT;
    1374:               col[7].c   = 0;
    1375:               valA[7]    = -1.0 / hz;
    1376:             }
    1377:           } else if (ex == N[0] - 1) {
    1378:             if (ey == 0) {
    1379:               nEntries   = 7;
    1380:               col[0].i   = ex;
    1381:               col[0].j   = ey;
    1382:               col[0].k   = ez;
    1383:               col[0].loc = BACK;
    1384:               col[0].c   = 0;
    1385:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1386:               /* Down term missing */
    1387:               col[1].i   = ex;
    1388:               col[1].j   = ey + 1;
    1389:               col[1].k   = ez;
    1390:               col[1].loc = BACK;
    1391:               col[1].c   = 0;
    1392:               valA[1]    = 1.0 / (hy * hy);
    1393:               col[2].i   = ex - 1;
    1394:               col[2].j   = ey;
    1395:               col[2].k   = ez;
    1396:               col[2].loc = BACK;
    1397:               col[2].c   = 0;
    1398:               valA[2]    = 1.0 / (hx * hx);
    1399:               /* Right term missing */
    1400:               col[3].i   = ex;
    1401:               col[3].j   = ey;
    1402:               col[3].k   = ez - 1;
    1403:               col[3].loc = BACK;
    1404:               col[3].c   = 0;
    1405:               valA[3]    = 1.0 / (hz * hz);
    1406:               col[4].i   = ex;
    1407:               col[4].j   = ey;
    1408:               col[4].k   = ez + 1;
    1409:               col[4].loc = BACK;
    1410:               col[4].c   = 0;
    1411:               valA[4]    = 1.0 / (hz * hz);
    1412:               col[5].i   = ex;
    1413:               col[5].j   = ey;
    1414:               col[5].k   = ez - 1;
    1415:               col[5].loc = ELEMENT;
    1416:               col[5].c   = 0;
    1417:               valA[5]    = 1.0 / hz;
    1418:               col[6].i   = ex;
    1419:               col[6].j   = ey;
    1420:               col[6].k   = ez;
    1421:               col[6].loc = ELEMENT;
    1422:               col[6].c   = 0;
    1423:               valA[6]    = -1.0 / hz;
    1424:             } else if (ey == N[1] - 1) {
    1425:               nEntries   = 7;
    1426:               col[0].i   = ex;
    1427:               col[0].j   = ey;
    1428:               col[0].k   = ez;
    1429:               col[0].loc = BACK;
    1430:               col[0].c   = 0;
    1431:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1432:               col[1].i   = ex;
    1433:               col[1].j   = ey - 1;
    1434:               col[1].k   = ez;
    1435:               col[1].loc = BACK;
    1436:               col[1].c   = 0;
    1437:               valA[1]    = 1.0 / (hy * hy);
    1438:               /* Up term missing */
    1439:               col[2].i   = ex - 1;
    1440:               col[2].j   = ey;
    1441:               col[2].k   = ez;
    1442:               col[2].loc = BACK;
    1443:               col[2].c   = 0;
    1444:               valA[2]    = 1.0 / (hx * hx);
    1445:               /* Right term missing */
    1446:               col[3].i   = ex;
    1447:               col[3].j   = ey;
    1448:               col[3].k   = ez - 1;
    1449:               col[3].loc = BACK;
    1450:               col[3].c   = 0;
    1451:               valA[3]    = 1.0 / (hz * hz);
    1452:               col[4].i   = ex;
    1453:               col[4].j   = ey;
    1454:               col[4].k   = ez + 1;
    1455:               col[4].loc = BACK;
    1456:               col[4].c   = 0;
    1457:               valA[4]    = 1.0 / (hz * hz);
    1458:               col[5].i   = ex;
    1459:               col[5].j   = ey;
    1460:               col[5].k   = ez - 1;
    1461:               col[5].loc = ELEMENT;
    1462:               col[5].c   = 0;
    1463:               valA[5]    = 1.0 / hz;
    1464:               col[6].i   = ex;
    1465:               col[6].j   = ey;
    1466:               col[6].k   = ez;
    1467:               col[6].loc = ELEMENT;
    1468:               col[6].c   = 0;
    1469:               valA[6]    = -1.0 / hz;
    1470:             } else {
    1471:               nEntries   = 8;
    1472:               col[0].i   = ex;
    1473:               col[0].j   = ey;
    1474:               col[0].k   = ez;
    1475:               col[0].loc = BACK;
    1476:               col[0].c   = 0;
    1477:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    1478:               col[1].i   = ex;
    1479:               col[1].j   = ey - 1;
    1480:               col[1].k   = ez;
    1481:               col[1].loc = BACK;
    1482:               col[1].c   = 0;
    1483:               valA[1]    = 1.0 / (hy * hy);
    1484:               col[2].i   = ex;
    1485:               col[2].j   = ey + 1;
    1486:               col[2].k   = ez;
    1487:               col[2].loc = BACK;
    1488:               col[2].c   = 0;
    1489:               valA[2]    = 1.0 / (hy * hy);
    1490:               col[3].i   = ex - 1;
    1491:               col[3].j   = ey;
    1492:               col[3].k   = ez;
    1493:               col[3].loc = BACK;
    1494:               col[3].c   = 0;
    1495:               valA[3]    = 1.0 / (hx * hx);
    1496:               /* Right term missing */
    1497:               col[4].i   = ex;
    1498:               col[4].j   = ey;
    1499:               col[4].k   = ez - 1;
    1500:               col[4].loc = BACK;
    1501:               col[4].c   = 0;
    1502:               valA[4]    = 1.0 / (hz * hz);
    1503:               col[5].i   = ex;
    1504:               col[5].j   = ey;
    1505:               col[5].k   = ez + 1;
    1506:               col[5].loc = BACK;
    1507:               col[5].c   = 0;
    1508:               valA[5]    = 1.0 / (hz * hz);
    1509:               col[6].i   = ex;
    1510:               col[6].j   = ey;
    1511:               col[6].k   = ez - 1;
    1512:               col[6].loc = ELEMENT;
    1513:               col[6].c   = 0;
    1514:               valA[6]    = 1.0 / hz;
    1515:               col[7].i   = ex;
    1516:               col[7].j   = ey;
    1517:               col[7].k   = ez;
    1518:               col[7].loc = ELEMENT;
    1519:               col[7].c   = 0;
    1520:               valA[7]    = -1.0 / hz;
    1521:             }
    1522:           } else if (ey == 0) {
    1523:             nEntries   = 8;
    1524:             col[0].i   = ex;
    1525:             col[0].j   = ey;
    1526:             col[0].k   = ez;
    1527:             col[0].loc = BACK;
    1528:             col[0].c   = 0;
    1529:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1530:             /* Down term missing */
    1531:             col[1].i   = ex;
    1532:             col[1].j   = ey + 1;
    1533:             col[1].k   = ez;
    1534:             col[1].loc = BACK;
    1535:             col[1].c   = 0;
    1536:             valA[1]    = 1.0 / (hy * hy);
    1537:             col[2].i   = ex - 1;
    1538:             col[2].j   = ey;
    1539:             col[2].k   = ez;
    1540:             col[2].loc = BACK;
    1541:             col[2].c   = 0;
    1542:             valA[2]    = 1.0 / (hx * hx);
    1543:             col[3].i   = ex + 1;
    1544:             col[3].j   = ey;
    1545:             col[3].k   = ez;
    1546:             col[3].loc = BACK;
    1547:             col[3].c   = 0;
    1548:             valA[3]    = 1.0 / (hx * hx);
    1549:             col[4].i   = ex;
    1550:             col[4].j   = ey;
    1551:             col[4].k   = ez - 1;
    1552:             col[4].loc = BACK;
    1553:             col[4].c   = 0;
    1554:             valA[4]    = 1.0 / (hz * hz);
    1555:             col[5].i   = ex;
    1556:             col[5].j   = ey;
    1557:             col[5].k   = ez + 1;
    1558:             col[5].loc = BACK;
    1559:             col[5].c   = 0;
    1560:             valA[5]    = 1.0 / (hz * hz);
    1561:             col[6].i   = ex;
    1562:             col[6].j   = ey;
    1563:             col[6].k   = ez - 1;
    1564:             col[6].loc = ELEMENT;
    1565:             col[6].c   = 0;
    1566:             valA[6]    = 1.0 / hz;
    1567:             col[7].i   = ex;
    1568:             col[7].j   = ey;
    1569:             col[7].k   = ez;
    1570:             col[7].loc = ELEMENT;
    1571:             col[7].c   = 0;
    1572:             valA[7]    = -1.0 / hz;
    1573:           } else if (ey == N[1] - 1) {
    1574:             nEntries   = 8;
    1575:             col[0].i   = ex;
    1576:             col[0].j   = ey;
    1577:             col[0].k   = ez;
    1578:             col[0].loc = BACK;
    1579:             col[0].c   = 0;
    1580:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
    1581:             col[1].i   = ex;
    1582:             col[1].j   = ey - 1;
    1583:             col[1].k   = ez;
    1584:             col[1].loc = BACK;
    1585:             col[1].c   = 0;
    1586:             valA[1]    = 1.0 / (hy * hy);
    1587:             /* Up term missing */
    1588:             col[2].i   = ex - 1;
    1589:             col[2].j   = ey;
    1590:             col[2].k   = ez;
    1591:             col[2].loc = BACK;
    1592:             col[2].c   = 0;
    1593:             valA[2]    = 1.0 / (hx * hx);
    1594:             col[3].i   = ex + 1;
    1595:             col[3].j   = ey;
    1596:             col[3].k   = ez;
    1597:             col[3].loc = BACK;
    1598:             col[3].c   = 0;
    1599:             valA[3]    = 1.0 / (hx * hx);
    1600:             col[4].i   = ex;
    1601:             col[4].j   = ey;
    1602:             col[4].k   = ez - 1;
    1603:             col[4].loc = BACK;
    1604:             col[4].c   = 0;
    1605:             valA[4]    = 1.0 / (hz * hz);
    1606:             col[5].i   = ex;
    1607:             col[5].j   = ey;
    1608:             col[5].k   = ez + 1;
    1609:             col[5].loc = BACK;
    1610:             col[5].c   = 0;
    1611:             valA[5]    = 1.0 / (hz * hz);
    1612:             col[6].i   = ex;
    1613:             col[6].j   = ey;
    1614:             col[6].k   = ez - 1;
    1615:             col[6].loc = ELEMENT;
    1616:             col[6].c   = 0;
    1617:             valA[6]    = 1.0 / hz;
    1618:             col[7].i   = ex;
    1619:             col[7].j   = ey;
    1620:             col[7].k   = ez;
    1621:             col[7].loc = ELEMENT;
    1622:             col[7].c   = 0;
    1623:             valA[7]    = -1.0 / hz;
    1624:           } else {
    1625:             nEntries   = 9;
    1626:             col[0].i   = ex;
    1627:             col[0].j   = ey;
    1628:             col[0].k   = ez;
    1629:             col[0].loc = BACK;
    1630:             col[0].c   = 0;
    1631:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
    1632:             col[1].i   = ex;
    1633:             col[1].j   = ey - 1;
    1634:             col[1].k   = ez;
    1635:             col[1].loc = BACK;
    1636:             col[1].c   = 0;
    1637:             valA[1]    = 1.0 / (hy * hy);
    1638:             col[2].i   = ex;
    1639:             col[2].j   = ey + 1;
    1640:             col[2].k   = ez;
    1641:             col[2].loc = BACK;
    1642:             col[2].c   = 0;
    1643:             valA[2]    = 1.0 / (hy * hy);
    1644:             col[3].i   = ex - 1;
    1645:             col[3].j   = ey;
    1646:             col[3].k   = ez;
    1647:             col[3].loc = BACK;
    1648:             col[3].c   = 0;
    1649:             valA[3]    = 1.0 / (hx * hx);
    1650:             col[4].i   = ex + 1;
    1651:             col[4].j   = ey;
    1652:             col[4].k   = ez;
    1653:             col[4].loc = BACK;
    1654:             col[4].c   = 0;
    1655:             valA[4]    = 1.0 / (hx * hx);
    1656:             col[5].i   = ex;
    1657:             col[5].j   = ey;
    1658:             col[5].k   = ez - 1;
    1659:             col[5].loc = BACK;
    1660:             col[5].c   = 0;
    1661:             valA[5]    = 1.0 / (hz * hz);
    1662:             col[6].i   = ex;
    1663:             col[6].j   = ey;
    1664:             col[6].k   = ez + 1;
    1665:             col[6].loc = BACK;
    1666:             col[6].c   = 0;
    1667:             valA[6]    = 1.0 / (hz * hz);
    1668:             col[7].i   = ex;
    1669:             col[7].j   = ey;
    1670:             col[7].k   = ez - 1;
    1671:             col[7].loc = ELEMENT;
    1672:             col[7].c   = 0;
    1673:             valA[7]    = 1.0 / hz;
    1674:             col[8].i   = ex;
    1675:             col[8].j   = ey;
    1676:             col[8].k   = ez;
    1677:             col[8].loc = ELEMENT;
    1678:             col[8].c   = 0;
    1679:             valA[8]    = -1.0 / hz;
    1680:           }
    1681:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
    1682:           valRhs = fz(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
    1683:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    1684:         }

    1686:         /* P equation : u_x + v_y + w_z = g
    1687:            Note that this includes an explicit zero on the diagonal. This is only needed for
    1688:            direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
    1689:         if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
    1690:           DMStagStencil row;
    1691:           PetscScalar   valA, valRhs;
    1692:           row.i   = ex;
    1693:           row.j   = ey;
    1694:           row.k   = ez;
    1695:           row.loc = ELEMENT;
    1696:           row.c   = 0;
    1697:           valA    = 1.0;
    1698:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
    1699:           valRhs = pRef(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
    1700:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    1701:         } else {
    1702:           DMStagStencil row, col[7];
    1703:           PetscScalar   valA[7], valRhs;

    1705:           row.i      = ex;
    1706:           row.j      = ey;
    1707:           row.k      = ez;
    1708:           row.loc    = ELEMENT;
    1709:           row.c      = 0;
    1710:           col[0].i   = ex;
    1711:           col[0].j   = ey;
    1712:           col[0].k   = ez;
    1713:           col[0].loc = LEFT;
    1714:           col[0].c   = 0;
    1715:           valA[0]    = -1.0 / hx;
    1716:           col[1].i   = ex;
    1717:           col[1].j   = ey;
    1718:           col[1].k   = ez;
    1719:           col[1].loc = RIGHT;
    1720:           col[1].c   = 0;
    1721:           valA[1]    = 1.0 / hx;
    1722:           col[2].i   = ex;
    1723:           col[2].j   = ey;
    1724:           col[2].k   = ez;
    1725:           col[2].loc = DOWN;
    1726:           col[2].c   = 0;
    1727:           valA[2]    = -1.0 / hy;
    1728:           col[3].i   = ex;
    1729:           col[3].j   = ey;
    1730:           col[3].k   = ez;
    1731:           col[3].loc = UP;
    1732:           col[3].c   = 0;
    1733:           valA[3]    = 1.0 / hy;
    1734:           col[4].i   = ex;
    1735:           col[4].j   = ey;
    1736:           col[4].k   = ez;
    1737:           col[4].loc = BACK;
    1738:           col[4].c   = 0;
    1739:           valA[4]    = -1.0 / hz;
    1740:           col[5].i   = ex;
    1741:           col[5].j   = ey;
    1742:           col[5].k   = ez;
    1743:           col[5].loc = FRONT;
    1744:           col[5].c   = 0;
    1745:           valA[5]    = 1.0 / hz;
    1746:           col[6]     = row;
    1747:           valA[6]    = 0.0;
    1748:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 7, col, valA, INSERT_VALUES));
    1749:           valRhs = g(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
    1750:           PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
    1751:         }
    1752:       }
    1753:     }
    1754:   }
    1755:   PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
    1756:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1757:   PetscCall(VecAssemblyBegin(rhs));
    1758:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    1759:   PetscCall(VecAssemblyEnd(rhs));
    1760:   PetscFunctionReturn(PETSC_SUCCESS);
    1761: }

    1763: /* Create a pressure-only DMStag and use it to generate a nullspace vector
    1764:    - Create a compatible DMStag with one dof per element (and nothing else).
    1765:    - Create a constant vector and normalize it
    1766:    - Migrate it to a vector on the original dmSol (making use of the fact
    1767:    that this will fill in zeros for "extra" dof)
    1768:    - Set the nullspace for the operator
    1769:    - Destroy everything (the operator keeps the references it needs) */
    1770: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
    1771: {
    1772:   DM           dmPressure;
    1773:   Vec          constantPressure, basis;
    1774:   PetscReal    nrm;
    1775:   MatNullSpace matNullSpace;

    1777:   PetscFunctionBeginUser;
    1778:   PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 0, 0, 1, &dmPressure));
    1779:   PetscCall(DMGetGlobalVector(dmPressure, &constantPressure));
    1780:   PetscCall(VecSet(constantPressure, 1.0));
    1781:   PetscCall(VecNorm(constantPressure, NORM_2, &nrm));
    1782:   PetscCall(VecScale(constantPressure, 1.0 / nrm));
    1783:   PetscCall(DMCreateGlobalVector(dmSol, &basis));
    1784:   PetscCall(DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis));
    1785:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace));
    1786:   PetscCall(VecDestroy(&basis));
    1787:   PetscCall(VecDestroy(&constantPressure));
    1788:   PetscCall(MatSetNullSpace(A, matNullSpace));
    1789:   PetscCall(MatNullSpaceDestroy(&matNullSpace));
    1790:   PetscFunctionReturn(PETSC_SUCCESS);
    1791: }

    1793: static PetscErrorCode CreateReferenceSolution(DM dmSol, Vec *pSolRef)
    1794: {
    1795:   PetscInt        start[3], n[3], nExtra[3], ex, ey, ez, d;
    1796:   PetscInt        ip, iux, iuy, iuz, icp[3], icux[3], icuy[3], icuz[3];
    1797:   Vec             solRef, solRefLocal, coord, coordLocal;
    1798:   DM              dmCoord;
    1799:   PetscScalar ****arrSol, ****arrCoord;

    1801:   PetscFunctionBeginUser;
    1802:   PetscCall(DMCreateGlobalVector(dmSol, pSolRef));
    1803:   solRef = *pSolRef;
    1804:   PetscCall(DMStagGetCorners(dmSol, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]));
    1805:   PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
    1806:   PetscCall(DMGetCoordinates(dmSol, &coord));
    1807:   PetscCall(DMGetLocalVector(dmCoord, &coordLocal));
    1808:   PetscCall(DMGlobalToLocal(dmCoord, coord, INSERT_VALUES, coordLocal));
    1809:   PetscCall(DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip));
    1810:   PetscCall(DMStagGetLocationSlot(dmSol, LEFT, 0, &iux));
    1811:   PetscCall(DMStagGetLocationSlot(dmSol, DOWN, 0, &iuy));
    1812:   PetscCall(DMStagGetLocationSlot(dmSol, BACK, 0, &iuz));
    1813:   for (d = 0; d < 3; ++d) {
    1814:     PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
    1815:     PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
    1816:     PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
    1817:     PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
    1818:   }
    1819:   PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
    1820:   PetscCall(DMGetLocalVector(dmSol, &solRefLocal));
    1821:   PetscCall(DMStagVecGetArray(dmSol, solRefLocal, &arrSol));
    1822:   for (ez = start[2]; ez < start[2] + n[2] + nExtra[2]; ++ez) {
    1823:     for (ey = start[1]; ey < start[1] + n[1] + nExtra[1]; ++ey) {
    1824:       for (ex = start[0]; ex < start[0] + n[0] + nExtra[0]; ++ex) {
    1825:         if (ex < start[0] + n[0] && ey < start[1] + n[1]) arrSol[ez][ey][ex][iuz] = uzRef(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
    1826:         if (ex < start[0] + n[0] && ey < start[2] + n[2]) arrSol[ez][ey][ex][iuy] = uyRef(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
    1827:         if (ex < start[1] + n[1] && ey < start[2] + n[2]) arrSol[ez][ey][ex][iux] = uxRef(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
    1828:         if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) arrSol[ez][ey][ex][ip] = pRef(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
    1829:       }
    1830:     }
    1831:   }
    1832:   PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
    1833:   PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
    1834:   PetscCall(DMLocalToGlobal(dmSol, solRefLocal, INSERT_VALUES, solRef));
    1835:   PetscCall(DMRestoreLocalVector(dmCoord, &coordLocal));
    1836:   PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
    1837:   PetscFunctionReturn(PETSC_SUCCESS);
    1838: }

    1840: static PetscErrorCode CheckSolution(Vec sol, Vec solRef)
    1841: {
    1842:   Vec       diff;
    1843:   PetscReal normsolRef, errAbs, errRel;

    1845:   PetscFunctionBeginUser;
    1846:   PetscCall(VecDuplicate(sol, &diff));
    1847:   PetscCall(VecCopy(sol, diff));
    1848:   PetscCall(VecAXPY(diff, -1.0, solRef));
    1849:   PetscCall(VecNorm(diff, NORM_2, &errAbs));
    1850:   PetscCall(VecNorm(solRef, NORM_2, &normsolRef));
    1851:   errRel = errAbs / normsolRef;
    1852:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel));
    1853:   PetscCall(VecDestroy(&diff));
    1854:   PetscFunctionReturn(PETSC_SUCCESS);
    1855: }

    1857: /*TEST

    1859:    test:
    1860:       suffix: 1
    1861:       requires: mumps
    1862:       nsize: 27
    1863:       args: -ksp_monitor_short -ksp_converged_reason -stag_ranks_x 3 -stag_ranks_y 3 -stag_ranks_z 3 -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type lu -fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20

    1865:    test:
    1866:       suffix: 2
    1867:       requires: !single
    1868:       nsize: 4
    1869:       args: -ksp_monitor_short -ksp_converged_reason -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type gamg -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_gamg_esteig_ksp_max_it 5

    1871:    test:
    1872:       suffix: direct_umfpack
    1873:       requires: suitesparse
    1874:       nsize: 1
    1875:       args: -pinpressure 1 -stag_grid_x 5 -stag_grid_y 3 -stag_grid_z 4 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack

    1877: TEST*/