Actual source code: ex4.c

      1: static char help[] = "Solves the incompressible, variable-viscosity Stokes equation in 2D or 3D, driven by buoyancy variations.\n"
      2:                      "-dim: set dimension (2 or 3)\n"
      3:                      "-nondimensional: replace dimensional domain and coefficients with nondimensional ones\n"
      4:                      "-isoviscous: use constant viscosity\n"
      5:                      "-levels: number of grids to create, by coarsening\n"
      6:                      "-rediscretize: create operators for all grids and set up a FieldSplit/MG solver\n"
      7:                      "-dump_solution: dump VTK files\n\n";

      9: #include <petscdm.h>
     10: #include <petscksp.h>
     11: #include <petscdmstag.h>
     12: #include <petscdmda.h>

     14: /* Application context - grid-based data*/
     15: typedef struct LevelCtxData_ {
     16:   DM          dm_stokes, dm_coefficients, dm_faces;
     17:   Vec         coeff;
     18:   PetscInt    cells_x, cells_y, cells_z; /* redundant with DMs */
     19:   PetscReal   hx_characteristic, hy_characteristic, hz_characteristic;
     20:   PetscScalar K_bound, K_cont;
     21: } LevelCtxData;
     22: typedef LevelCtxData *LevelCtx;

     24: /* Application context - problem and grid(s) (but not solver-specific data) */
     25: typedef struct CtxData_ {
     26:   MPI_Comm    comm;
     27:   PetscInt    dim;                       /* redundant with DMs */
     28:   PetscInt    cells_x, cells_y, cells_z; /* Redundant with finest DMs */
     29:   PetscReal   xmax, ymax, xmin, ymin, zmin, zmax;
     30:   PetscScalar eta1, eta2, rho1, rho2, gy, eta_characteristic;
     31:   PetscBool   pin_pressure;
     32:   PetscScalar (*GetEta)(struct CtxData_ *, PetscScalar, PetscScalar, PetscScalar);
     33:   PetscScalar (*GetRho)(struct CtxData_ *, PetscScalar, PetscScalar, PetscScalar);
     34:   PetscInt  n_levels;
     35:   LevelCtx *levels;
     36: } CtxData;
     37: typedef CtxData *Ctx;

     39: /* Helper to pass system-creation parameters */
     40: typedef struct SystemParameters_ {
     41:   Ctx       ctx;
     42:   PetscInt  level;
     43:   PetscBool include_inverse_visc, faces_only;
     44: } SystemParametersData;
     45: typedef SystemParametersData *SystemParameters;

     47: /* Main logic */
     48: static PetscErrorCode AttachNullspace(DM, Mat);
     49: static PetscErrorCode CreateAuxiliaryOperator(Ctx, PetscInt, Mat *);
     50: static PetscErrorCode CreateSystem(SystemParameters, Mat *, Vec *);
     51: static PetscErrorCode CtxCreateAndSetFromOptions(Ctx *);
     52: static PetscErrorCode CtxDestroy(Ctx *);
     53: static PetscErrorCode DumpSolution(Ctx, PetscInt, Vec);
     54: static PetscErrorCode OperatorInsertInverseViscosityPressureTerms(DM, DM, Vec, PetscScalar, Mat);
     55: static PetscErrorCode PopulateCoefficientData(Ctx, PetscInt);
     56: static PetscErrorCode SystemParametersCreate(SystemParameters *, Ctx);
     57: static PetscErrorCode SystemParametersDestroy(SystemParameters *);

     59: int main(int argc, char **argv)
     60: {
     61:   Ctx       ctx;
     62:   Mat       A, *A_faces = NULL, S_hat = NULL, P = NULL;
     63:   Vec       x, b;
     64:   KSP       ksp;
     65:   DM        dm_stokes, dm_coefficients;
     66:   PetscBool dump_solution, build_auxiliary_operator, rediscretize, custom_pc_mat;

     68:   PetscFunctionBeginUser;
     69:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     71:   /* Accept options for program behavior */
     72:   dump_solution = PETSC_FALSE;
     73:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dump_solution", &dump_solution, NULL));
     74:   rediscretize = PETSC_FALSE;
     75:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-rediscretize", &rediscretize, NULL));
     76:   build_auxiliary_operator = rediscretize;
     77:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-build_auxiliary_operator", &build_auxiliary_operator, NULL));
     78:   custom_pc_mat = PETSC_FALSE;
     79:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-custom_pc_mat", &custom_pc_mat, NULL));

     81:   /* Populate application context */
     82:   PetscCall(CtxCreateAndSetFromOptions(&ctx));

     84:   /* Create two DMStag objects, corresponding to the same domain and parallel
     85:      decomposition ("topology"). Each defines a different set of fields on
     86:      the domain ("section"); the first the solution to the Stokes equations
     87:      (x- and y-velocities and scalar pressure), and the second holds coefficients
     88:      (viscosities on elements and viscosities+densities on corners/edges in 2d/3d) */
     89:   if (ctx->dim == 2) {
     90:     PetscCall(DMStagCreate2d(ctx->comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, ctx->cells_x, ctx->cells_y, /* Global element counts */
     91:                              PETSC_DECIDE, PETSC_DECIDE,                                                /* Determine parallel decomposition automatically */
     92:                              0, 1, 1,                                                                   /* dof: 0 per vertex, 1 per edge, 1 per face/element */
     93:                              DMSTAG_STENCIL_BOX, 1,                                                     /* elementwise stencil width */
     94:                              NULL, NULL, &ctx->levels[ctx->n_levels - 1]->dm_stokes));
     95:   } else if (ctx->dim == 3) {
     96:     PetscCall(DMStagCreate3d(ctx->comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, ctx->cells_x, ctx->cells_y, ctx->cells_z, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 0, 0, 1, 1, /* dof: 1 per face, 1 per element */
     97:                              DMSTAG_STENCIL_BOX, 1, NULL, NULL, NULL, &ctx->levels[ctx->n_levels - 1]->dm_stokes));
     98:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension: %" PetscInt_FMT, ctx->dim);
     99:   dm_stokes = ctx->levels[ctx->n_levels - 1]->dm_stokes;
    100:   PetscCall(DMSetFromOptions(dm_stokes));
    101:   PetscCall(DMStagGetGlobalSizes(dm_stokes, &ctx->cells_x, &ctx->cells_y, &ctx->cells_z));
    102:   PetscCall(DMSetUp(dm_stokes));
    103:   PetscCall(DMStagSetUniformCoordinatesProduct(dm_stokes, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax, ctx->zmin, ctx->zmax));

    105:   if (ctx->dim == 2) PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 2, 0, 1, 0, &ctx->levels[ctx->n_levels - 1]->dm_coefficients));
    106:   else if (ctx->dim == 3) PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 2, 0, 1, &ctx->levels[ctx->n_levels - 1]->dm_coefficients));
    107:   else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension: %" PetscInt_FMT, ctx->dim);
    108:   dm_coefficients = ctx->levels[ctx->n_levels - 1]->dm_coefficients;
    109:   PetscCall(DMSetUp(dm_coefficients));
    110:   PetscCall(DMStagSetUniformCoordinatesProduct(dm_coefficients, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax, ctx->zmin, ctx->zmax));

    112:   /* Create additional DMs by coarsening. 0 is the coarsest level */
    113:   for (PetscInt level = ctx->n_levels - 1; level > 0; --level) {
    114:     PetscCall(DMCoarsen(ctx->levels[level]->dm_stokes, MPI_COMM_NULL, &ctx->levels[level - 1]->dm_stokes));
    115:     PetscCall(DMCoarsen(ctx->levels[level]->dm_coefficients, MPI_COMM_NULL, &ctx->levels[level - 1]->dm_coefficients));
    116:   }

    118:   /* Compute scaling constants, knowing grid spacing */
    119:   ctx->eta_characteristic = PetscMin(PetscRealPart(ctx->eta1), PetscRealPart(ctx->eta2));
    120:   for (PetscInt level = 0; level < ctx->n_levels; ++level) {
    121:     PetscInt  N[3];
    122:     PetscReal hx_avg_inv;

    124:     PetscCall(DMStagGetGlobalSizes(ctx->levels[level]->dm_stokes, &N[0], &N[1], &N[2]));
    125:     ctx->levels[level]->hx_characteristic = (ctx->xmax - ctx->xmin) / N[0];
    126:     ctx->levels[level]->hy_characteristic = (ctx->ymax - ctx->ymin) / N[1];
    127:     if (N[2]) ctx->levels[level]->hz_characteristic = (ctx->zmax - ctx->zmin) / N[2];
    128:     if (ctx->dim == 2) {
    129:       hx_avg_inv = 2.0 / (ctx->levels[level]->hx_characteristic + ctx->levels[level]->hy_characteristic);
    130:     } else if (ctx->dim == 3) {
    131:       hx_avg_inv = 3.0 / (ctx->levels[level]->hx_characteristic + ctx->levels[level]->hy_characteristic + ctx->levels[level]->hz_characteristic);
    132:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not Implemented for dimension %" PetscInt_FMT, ctx->dim);
    133:     ctx->levels[level]->K_cont  = ctx->eta_characteristic * hx_avg_inv;
    134:     ctx->levels[level]->K_bound = ctx->eta_characteristic * hx_avg_inv * hx_avg_inv;
    135:   }

    137:   /* Populate coefficient data */
    138:   for (PetscInt level = 0; level < ctx->n_levels; ++level) PetscCall(PopulateCoefficientData(ctx, level));

    140:   /* Construct main system */
    141:   {
    142:     SystemParameters system_parameters;

    144:     PetscCall(SystemParametersCreate(&system_parameters, ctx));
    145:     PetscCall(CreateSystem(system_parameters, &A, &b));
    146:     PetscCall(SystemParametersDestroy(&system_parameters));
    147:   }

    149:   /* Attach a constant-pressure nullspace to the fine-level operator */
    150:   if (!ctx->pin_pressure) PetscCall(AttachNullspace(dm_stokes, A));

    152:   /* Set up solver */
    153:   PetscCall(KSPCreate(ctx->comm, &ksp));
    154:   PetscCall(KSPSetType(ksp, KSPFGMRES));
    155:   {
    156:     /* Default to a direct solver, if a package is available */
    157:     PetscMPIInt size;

    159:     PetscCallMPI(MPI_Comm_size(ctx->comm, &size));
    160:     if (PetscDefined(HAVE_SUITESPARSE) && size == 1) {
    161:       PC pc;

    163:       PetscCall(KSPGetPC(ksp, &pc));
    164:       PetscCall(PCSetType(pc, PCLU));
    165:       PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK)); /* A default, requires SuiteSparse */
    166:     }
    167:     if (PetscDefined(HAVE_MUMPS) && size > 1) {
    168:       PC pc;

    170:       PetscCall(KSPGetPC(ksp, &pc));
    171:       PetscCall(PCSetType(pc, PCLU));
    172:       PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); /* A default, requires MUMPS */
    173:     }
    174:   }

    176:   /* Create and set a custom preconditioning matrix */
    177:   if (custom_pc_mat) {
    178:     SystemParameters system_parameters;

    180:     PetscCall(SystemParametersCreate(&system_parameters, ctx));
    181:     system_parameters->include_inverse_visc = PETSC_TRUE;
    182:     PetscCall(CreateSystem(system_parameters, &P, NULL));
    183:     PetscCall(SystemParametersDestroy(&system_parameters));
    184:     PetscCall(KSPSetOperators(ksp, A, P));
    185:   } else {
    186:     PetscCall(KSPSetOperators(ksp, A, A));
    187:   }

    189:   PetscCall(KSPSetDM(ksp, dm_stokes));
    190:   PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));

    192:   /* Finish setting up solver (can override options set above) */
    193:   PetscCall(KSPSetFromOptions(ksp));

    195:   /* Additional solver configuration that can involve setting up and CANNOT be
    196:      overridden from the command line */

    198:   /* Construct an auxiliary operator for use a Schur complement preconditioner,
    199:      and tell PCFieldSplit to use it (which has no effect if not using that PC) */
    200:   if (build_auxiliary_operator && !rediscretize) {
    201:     PC pc;

    203:     PetscCall(KSPGetPC(ksp, &pc));
    204:     PetscCall(CreateAuxiliaryOperator(ctx, ctx->n_levels - 1, &S_hat));
    205:     PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S_hat));
    206:     PetscCall(KSPSetFromOptions(ksp));
    207:   }

    209:   if (rediscretize) {
    210:     /* Set up an ABF solver with rediscretized geometric multigrid on the velocity-velocity block */
    211:     PC  pc, pc_faces;
    212:     KSP ksp_faces;

    214:     if (ctx->n_levels < 2) PetscCall(PetscPrintf(ctx->comm, "Warning: not using multiple levels!\n"));

    216:     PetscCall(KSPGetPC(ksp, &pc));
    217:     PetscCall(PCSetType(pc, PCFIELDSPLIT));
    218:     PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
    219:     PetscCall(PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_UPPER));
    220:     if (build_auxiliary_operator) {
    221:       PetscCall(CreateAuxiliaryOperator(ctx, ctx->n_levels - 1, &S_hat));
    222:       PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S_hat));
    223:     }

    225:     /* Create rediscretized velocity-only DMs and operators */
    226:     PetscCall(PetscMalloc1(ctx->n_levels, &A_faces));
    227:     for (PetscInt level = 0; level < ctx->n_levels; ++level) {
    228:       if (ctx->dim == 2) {
    229:         PetscCall(DMStagCreateCompatibleDMStag(ctx->levels[level]->dm_stokes, 0, 1, 0, 0, &ctx->levels[level]->dm_faces));
    230:       } else if (ctx->dim == 3) {
    231:         PetscCall(DMStagCreateCompatibleDMStag(ctx->levels[level]->dm_stokes, 0, 0, 1, 0, &ctx->levels[level]->dm_faces));
    232:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not Implemented for dimension %" PetscInt_FMT, ctx->dim);
    233:       {
    234:         SystemParameters system_parameters;

    236:         PetscCall(SystemParametersCreate(&system_parameters, ctx));
    237:         system_parameters->faces_only = PETSC_TRUE;
    238:         system_parameters->level      = level;
    239:         PetscCall(CreateSystem(system_parameters, &A_faces[level], NULL));
    240:         PetscCall(SystemParametersDestroy(&system_parameters));
    241:       }
    242:     }

    244:     /* Set up to populate enough to define the sub-solver */
    245:     PetscCall(KSPSetUp(ksp));

    247:     /* Set multigrid components and settings */
    248:     {
    249:       KSP *sub_ksp;

    251:       PetscCall(PCFieldSplitSchurGetSubKSP(pc, NULL, &sub_ksp));
    252:       ksp_faces = sub_ksp[0];
    253:       PetscCall(PetscFree(sub_ksp));
    254:     }
    255:     PetscCall(KSPSetType(ksp_faces, KSPGCR));
    256:     PetscCall(KSPGetPC(ksp_faces, &pc_faces));
    257:     PetscCall(PCSetType(pc_faces, PCMG));
    258:     PetscCall(PCMGSetLevels(pc_faces, ctx->n_levels, NULL));
    259:     for (PetscInt level = 0; level < ctx->n_levels; ++level) {
    260:       KSP ksp_level;
    261:       PC  pc_level;

    263:       /* Smoothers */
    264:       PetscCall(PCMGGetSmoother(pc_faces, level, &ksp_level));
    265:       PetscCall(KSPGetPC(ksp_level, &pc_level));
    266:       PetscCall(KSPSetOperators(ksp_level, A_faces[level], A_faces[level]));
    267:       if (level > 0) PetscCall(PCSetType(pc_level, PCJACOBI));

    269:       /* Transfer Operators */
    270:       if (level > 0) {
    271:         Mat restriction, interpolation;
    272:         DM  dm_level   = ctx->levels[level]->dm_faces;
    273:         DM  dm_coarser = ctx->levels[level - 1]->dm_faces;

    275:         PetscCall(DMCreateInterpolation(dm_coarser, dm_level, &interpolation, NULL));
    276:         PetscCall(PCMGSetInterpolation(pc_faces, level, interpolation));
    277:         PetscCall(MatDestroy(&interpolation));
    278:         PetscCall(DMCreateRestriction(dm_coarser, dm_level, &restriction));
    279:         PetscCall(PCMGSetRestriction(pc_faces, level, restriction));
    280:         PetscCall(MatDestroy(&restriction));
    281:       }
    282:     }
    283:   }

    285:   /* Solve */
    286:   PetscCall(VecDuplicate(b, &x));
    287:   PetscCall(KSPSolve(ksp, b, x));
    288:   {
    289:     KSPConvergedReason reason;
    290:     PetscCall(KSPGetConvergedReason(ksp, &reason));
    291:     PetscCheck(reason >= 0, PETSC_COMM_WORLD, PETSC_ERR_CONV_FAILED, "Linear solve failed");
    292:   }

    294:   /* Dump solution by converting to DMDAs and dumping */
    295:   if (dump_solution) PetscCall(DumpSolution(ctx, ctx->n_levels - 1, x));

    297:   /* Destroy PETSc objects and finalize */
    298:   PetscCall(MatDestroy(&A));
    299:   PetscCall(PetscFree(A));
    300:   if (A_faces) {
    301:     for (PetscInt level = 0; level < ctx->n_levels; ++level) {
    302:       if (A_faces[level]) PetscCall(MatDestroy(&A_faces[level]));
    303:     }
    304:     PetscCall(PetscFree(A_faces));
    305:   }
    306:   if (P) PetscCall(MatDestroy(&P));
    307:   PetscCall(VecDestroy(&x));
    308:   PetscCall(VecDestroy(&b));
    309:   PetscCall(MatDestroy(&S_hat));
    310:   PetscCall(KSPDestroy(&ksp));
    311:   PetscCall(CtxDestroy(&ctx));
    312:   PetscCall(PetscFinalize());
    313:   return 0;
    314: }

    316: static PetscScalar GetEta_constant(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    317: {
    318:   (void)ctx;
    319:   (void)x;
    320:   (void)y;
    321:   (void)z;
    322:   return 1.0;
    323: }

    325: static PetscScalar GetRho_layers(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    326: {
    327:   (void)y;
    328:   (void)z;
    329:   return PetscRealPart(x) < (ctx->xmax - ctx->xmin) / 2.0 ? ctx->rho1 : ctx->rho2;
    330: }

    332: static PetscScalar GetEta_layers(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    333: {
    334:   (void)y;
    335:   (void)z;
    336:   return PetscRealPart(x) < (ctx->xmax - ctx->xmin) / 2.0 ? ctx->eta1 : ctx->eta2;
    337: }

    339: static PetscScalar GetRho_sinker_box2(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    340: {
    341:   (void)z;
    342:   const PetscReal d  = ctx->xmax - ctx->xmin;
    343:   const PetscReal xx = PetscRealPart(x) / d - 0.5;
    344:   const PetscReal yy = PetscRealPart(y) / d - 0.5;
    345:   return (xx * xx > 0.15 * 0.15 || yy * yy > 0.15 * 0.15) ? ctx->rho1 : ctx->rho2;
    346: }

    348: static PetscScalar GetEta_sinker_box2(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    349: {
    350:   (void)z;
    351:   const PetscReal d  = ctx->xmax - ctx->xmin;
    352:   const PetscReal xx = PetscRealPart(x) / d - 0.5;
    353:   const PetscReal yy = PetscRealPart(y) / d - 0.5;
    354:   return (xx * xx > 0.15 * 0.15 || yy * yy > 0.15 * 0.15) ? ctx->eta1 : ctx->eta2;
    355: }

    357: static PetscScalar GetRho_sinker_box3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    358: {
    359:   const PetscReal d          = ctx->xmax - ctx->xmin;
    360:   const PetscReal xx         = PetscRealPart(x) / d - 0.5;
    361:   const PetscReal yy         = PetscRealPart(y) / d - 0.5;
    362:   const PetscReal zz         = PetscRealPart(z) / d - 0.5;
    363:   const PetscReal half_width = 0.15;
    364:   return (PetscAbsReal(xx) > half_width || PetscAbsReal(yy) > half_width || PetscAbsReal(zz) > half_width) ? ctx->rho1 : ctx->rho2;
    365: }

    367: static PetscScalar GetEta_sinker_box3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    368: {
    369:   const PetscReal d          = ctx->xmax - ctx->xmin;
    370:   const PetscReal xx         = PetscRealPart(x) / d - 0.5;
    371:   const PetscReal yy         = PetscRealPart(y) / d - 0.5;
    372:   const PetscReal zz         = PetscRealPart(z) / d - 0.5;
    373:   const PetscReal half_width = 0.15;
    374:   return (PetscAbsReal(xx) > half_width || PetscAbsReal(yy) > half_width || PetscAbsReal(zz) > half_width) ? ctx->eta1 : ctx->eta2;
    375: }

    377: static PetscScalar GetRho_sinker_sphere3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    378: {
    379:   const PetscReal d          = ctx->xmax - ctx->xmin;
    380:   const PetscReal xx         = PetscRealPart(x) / d - 0.5;
    381:   const PetscReal yy         = PetscRealPart(y) / d - 0.5;
    382:   const PetscReal zz         = PetscRealPart(z) / d - 0.5;
    383:   const PetscReal half_width = 0.3;
    384:   return (xx * xx + yy * yy + zz * zz > half_width * half_width) ? ctx->rho1 : ctx->rho2;
    385: }

    387: static PetscScalar GetEta_sinker_sphere3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    388: {
    389:   const PetscReal d          = ctx->xmax - ctx->xmin;
    390:   const PetscReal xx         = PetscRealPart(x) / d - 0.5;
    391:   const PetscReal yy         = PetscRealPart(y) / d - 0.5;
    392:   const PetscReal zz         = PetscRealPart(z) / d - 0.5;
    393:   const PetscReal half_width = 0.3;
    394:   return (xx * xx + yy * yy + zz * zz > half_width * half_width) ? ctx->eta1 : ctx->eta2;
    395: }

    397: static PetscScalar GetEta_blob3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    398: {
    399:   const PetscReal d  = ctx->xmax - ctx->xmin;
    400:   const PetscReal xx = PetscRealPart(x) / d - 0.5;
    401:   const PetscReal yy = PetscRealPart(y) / d - 0.5;
    402:   const PetscReal zz = PetscRealPart(z) / d - 0.5;
    403:   return ctx->eta1 + ctx->eta2 * PetscExpScalar(-20.0 * (xx * xx + yy * yy + zz * zz));
    404: }

    406: static PetscScalar GetRho_blob3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
    407: {
    408:   const PetscReal d  = ctx->xmax - ctx->xmin;
    409:   const PetscReal xx = PetscRealPart(x) / d - 0.5;
    410:   const PetscReal yy = PetscRealPart(y) / d - 0.5;
    411:   const PetscReal zz = PetscRealPart(z) / d - 0.5;
    412:   return ctx->rho1 + ctx->rho2 * PetscExpScalar(-20.0 * (xx * xx + yy * yy + zz * zz));
    413: }

    415: static PetscErrorCode LevelCtxCreate(LevelCtx *p_level_ctx)
    416: {
    417:   LevelCtx level_ctx;

    419:   PetscFunctionBeginUser;
    420:   PetscCall(PetscMalloc1(1, p_level_ctx));
    421:   level_ctx                  = *p_level_ctx;
    422:   level_ctx->dm_stokes       = NULL;
    423:   level_ctx->dm_coefficients = NULL;
    424:   level_ctx->dm_faces        = NULL;
    425:   PetscFunctionReturn(PETSC_SUCCESS);
    426: }

    428: static PetscErrorCode LevelCtxDestroy(LevelCtx *p_level_ctx)
    429: {
    430:   LevelCtx level_ctx;

    432:   PetscFunctionBeginUser;
    433:   level_ctx = *p_level_ctx;
    434:   if (level_ctx->dm_stokes) PetscCall(DMDestroy(&level_ctx->dm_stokes));
    435:   if (level_ctx->dm_coefficients) PetscCall(DMDestroy(&level_ctx->dm_coefficients));
    436:   if (level_ctx->dm_faces) PetscCall(DMDestroy(&level_ctx->dm_faces));
    437:   if (level_ctx->coeff) PetscCall(VecDestroy(&level_ctx->coeff));
    438:   PetscCall(PetscFree(*p_level_ctx));
    439:   PetscFunctionReturn(PETSC_SUCCESS);
    440: }

    442: static PetscErrorCode CtxCreateAndSetFromOptions(Ctx *p_ctx)
    443: {
    444:   Ctx ctx;

    446:   PetscFunctionBeginUser;
    447:   PetscCall(PetscMalloc1(1, p_ctx));
    448:   ctx = *p_ctx;

    450:   ctx->comm         = PETSC_COMM_WORLD;
    451:   ctx->pin_pressure = PETSC_FALSE;
    452:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pin_pressure", &ctx->pin_pressure, NULL));
    453:   ctx->dim = 3;
    454:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &ctx->dim, NULL));
    455:   if (ctx->dim <= 2) {
    456:     ctx->cells_x = 32;
    457:   } else {
    458:     ctx->cells_x = 16;
    459:   }
    460:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-s", &ctx->cells_x, NULL)); /* shortcut. Usually, use -stag_grid_x etc. */
    461:   ctx->cells_z = ctx->cells_y = ctx->cells_x;
    462:   ctx->xmin = ctx->ymin = ctx->zmin = 0.0;
    463:   {
    464:     PetscBool nondimensional = PETSC_TRUE;

    466:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-nondimensional", &nondimensional, NULL));
    467:     if (nondimensional) {
    468:       ctx->xmax = ctx->ymax = ctx->zmax = 1.0;
    469:       ctx->rho1                         = 0.0;
    470:       ctx->rho2                         = 1.0;
    471:       ctx->eta1                         = 1.0;
    472:       ctx->eta2                         = 1e2;
    473:       ctx->gy                           = -1.0; /* downwards */
    474:     } else {
    475:       ctx->xmax = 1e6;
    476:       ctx->ymax = 1.5e6;
    477:       ctx->zmax = 1e6;
    478:       ctx->rho1 = 3200;
    479:       ctx->rho2 = 3300;
    480:       ctx->eta1 = 1e20;
    481:       ctx->eta2 = 1e22;
    482:       ctx->gy   = -10.0; /* downwards */
    483:     }
    484:   }
    485:   {
    486:     PetscBool isoviscous;

    488:     isoviscous = PETSC_FALSE;
    489:     PetscCall(PetscOptionsGetScalar(NULL, NULL, "-eta1", &ctx->eta1, NULL));
    490:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-isoviscous", &isoviscous, NULL));
    491:     if (isoviscous) {
    492:       ctx->eta2   = ctx->eta1;
    493:       ctx->GetEta = GetEta_constant; /* override */
    494:     } else {
    495:       PetscCall(PetscOptionsGetScalar(NULL, NULL, "-eta2", &ctx->eta2, NULL));
    496:     }
    497:   }
    498:   {
    499:     char      mode[1024] = "sinker";
    500:     PetscBool is_layers, is_blob, is_sinker_box, is_sinker_sphere;

    502:     PetscCall(PetscOptionsGetString(NULL, NULL, "-coefficients", mode, sizeof(mode), NULL));
    503:     PetscCall(PetscStrncmp(mode, "layers", sizeof(mode), &is_layers));
    504:     PetscCall(PetscStrncmp(mode, "sinker", sizeof(mode), &is_sinker_box));
    505:     if (!is_sinker_box) PetscCall(PetscStrncmp(mode, "sinker_box", sizeof(mode), &is_sinker_box));
    506:     PetscCall(PetscStrncmp(mode, "sinker_sphere", sizeof(mode), &is_sinker_sphere));
    507:     PetscCall(PetscStrncmp(mode, "blob", sizeof(mode), &is_blob));

    509:     if (is_layers) {
    510:       ctx->GetRho = GetRho_layers;
    511:       ctx->GetEta = GetEta_layers;
    512:     }
    513:     if (is_blob) {
    514:       if (ctx->dim == 3) {
    515:         ctx->GetRho = GetRho_blob3;
    516:         ctx->GetEta = GetEta_blob3;
    517:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
    518:     }
    519:     if (is_sinker_box) {
    520:       if (ctx->dim == 2) {
    521:         ctx->GetRho = GetRho_sinker_box2;
    522:         ctx->GetEta = GetEta_sinker_box2;
    523:       } else if (ctx->dim == 3) {
    524:         ctx->GetRho = GetRho_sinker_box3;
    525:         ctx->GetEta = GetEta_sinker_box3;
    526:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
    527:     }
    528:     if (is_sinker_sphere) {
    529:       if (ctx->dim == 3) {
    530:         ctx->GetRho = GetRho_sinker_sphere3;
    531:         ctx->GetEta = GetEta_sinker_sphere3;
    532:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
    533:     }
    534:   }

    536:   /* Per-level data */
    537:   ctx->n_levels = 1;
    538:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-levels", &ctx->n_levels, NULL));
    539:   PetscCall(PetscMalloc1(ctx->n_levels, &ctx->levels));
    540:   for (PetscInt i = 0; i < ctx->n_levels; ++i) PetscCall(LevelCtxCreate(&ctx->levels[i]));
    541:   PetscFunctionReturn(PETSC_SUCCESS);
    542: }

    544: static PetscErrorCode CtxDestroy(Ctx *p_ctx)
    545: {
    546:   Ctx ctx;

    548:   PetscFunctionBeginUser;
    549:   ctx = *p_ctx;
    550:   for (PetscInt i = 0; i < ctx->n_levels; ++i) PetscCall(LevelCtxDestroy(&ctx->levels[i]));
    551:   PetscCall(PetscFree(ctx->levels));
    552:   PetscCall(PetscFree(*p_ctx));
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: static PetscErrorCode SystemParametersCreate(SystemParameters *parameters, Ctx ctx)
    557: {
    558:   PetscFunctionBeginUser;
    559:   PetscCall(PetscMalloc1(1, parameters));
    560:   (*parameters)->ctx                  = ctx;
    561:   (*parameters)->level                = ctx->n_levels - 1;
    562:   (*parameters)->include_inverse_visc = PETSC_FALSE;
    563:   (*parameters)->faces_only           = PETSC_FALSE;
    564:   PetscFunctionReturn(PETSC_SUCCESS);
    565: }

    567: static PetscErrorCode SystemParametersDestroy(SystemParameters *parameters)
    568: {
    569:   PetscFunctionBeginUser;
    570:   PetscCall(PetscFree(*parameters));
    571:   *parameters = NULL;
    572:   PetscFunctionReturn(PETSC_SUCCESS);
    573: }

    575: static PetscErrorCode CreateSystem2d(SystemParameters parameters, Mat *pA, Vec *pRhs)
    576: {
    577:   PetscInt    N[2];
    578:   PetscInt    ex, ey, startx, starty, nx, ny;
    579:   Mat         A;
    580:   Vec         rhs;
    581:   PetscReal   hx, hy, dv;
    582:   Vec         coefficients_local;
    583:   PetscBool   build_rhs;
    584:   DM          dm_main, dm_coefficients;
    585:   PetscScalar K_cont, K_bound;
    586:   Ctx         ctx   = parameters->ctx;
    587:   PetscInt    level = parameters->level;

    589:   PetscFunctionBeginUser;
    590:   if (parameters->faces_only) {
    591:     dm_main = ctx->levels[level]->dm_faces;
    592:   } else {
    593:     dm_main = ctx->levels[level]->dm_stokes;
    594:   }
    595:   dm_coefficients = ctx->levels[level]->dm_coefficients;
    596:   K_cont          = ctx->levels[level]->K_cont;
    597:   K_bound         = ctx->levels[level]->K_bound;
    598:   PetscCall(DMCreateMatrix(dm_main, pA));
    599:   A         = *pA;
    600:   build_rhs = (PetscBool)(pRhs != NULL);
    601:   PetscCheck(!(parameters->faces_only && build_rhs), PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "RHS for faces-only not supported");
    602:   if (build_rhs) {
    603:     PetscCall(DMCreateGlobalVector(dm_main, pRhs));
    604:     rhs = *pRhs;
    605:   } else {
    606:     rhs = NULL;
    607:   }
    608:   PetscCall(DMStagGetCorners(dm_main, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
    609:   PetscCall(DMStagGetGlobalSizes(dm_main, &N[0], &N[1], NULL));
    610:   hx = ctx->levels[level]->hx_characteristic;
    611:   hy = ctx->levels[level]->hy_characteristic;
    612:   dv = hx * hy;
    613:   PetscCall(DMGetLocalVector(dm_coefficients, &coefficients_local));
    614:   PetscCall(DMGlobalToLocal(dm_coefficients, ctx->levels[level]->coeff, INSERT_VALUES, coefficients_local));

    616:   /* Loop over all local elements */
    617:   for (ey = starty; ey < starty + ny; ++ey) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
    618:     for (ex = startx; ex < startx + nx; ++ex) {
    619:       const PetscBool left_boundary   = (PetscBool)(ex == 0);
    620:       const PetscBool right_boundary  = (PetscBool)(ex == N[0] - 1);
    621:       const PetscBool bottom_boundary = (PetscBool)(ey == 0);
    622:       const PetscBool top_boundary    = (PetscBool)(ey == N[1] - 1);

    624:       if (ey == N[1] - 1) {
    625:         /* Top boundary velocity Dirichlet */
    626:         DMStagStencil     row;
    627:         const PetscScalar val_rhs = 0.0;
    628:         const PetscScalar val_A   = K_bound;

    630:         row.i   = ex;
    631:         row.j   = ey;
    632:         row.loc = DMSTAG_UP;
    633:         row.c   = 0;
    634:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    635:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    636:       }

    638:       if (ey == 0) {
    639:         /* Bottom boundary velocity Dirichlet */
    640:         DMStagStencil     row;
    641:         const PetscScalar val_rhs = 0.0;
    642:         const PetscScalar val_A   = K_bound;

    644:         row.i   = ex;
    645:         row.j   = ey;
    646:         row.loc = DMSTAG_DOWN;
    647:         row.c   = 0;
    648:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    649:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    650:       } else {
    651:         /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y
    652:            includes non-zero forcing and free-slip boundary conditions */
    653:         PetscInt      count;
    654:         DMStagStencil row, col[11];
    655:         PetscScalar   val_A[11];
    656:         DMStagStencil rhoPoint[2];
    657:         PetscScalar   rho[2], val_rhs;
    658:         DMStagStencil etaPoint[4];
    659:         PetscScalar   eta[4], eta_left, eta_right, eta_up, eta_down;

    661:         row.i   = ex;
    662:         row.j   = ey;
    663:         row.loc = DMSTAG_DOWN;
    664:         row.c   = 0;

    666:         /* get rho values  and compute rhs value*/
    667:         rhoPoint[0].i   = ex;
    668:         rhoPoint[0].j   = ey;
    669:         rhoPoint[0].loc = DMSTAG_DOWN_LEFT;
    670:         rhoPoint[0].c   = 1;
    671:         rhoPoint[1].i   = ex;
    672:         rhoPoint[1].j   = ey;
    673:         rhoPoint[1].loc = DMSTAG_DOWN_RIGHT;
    674:         rhoPoint[1].c   = 1;
    675:         PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 2, rhoPoint, rho));
    676:         val_rhs = -ctx->gy * dv * 0.5 * (rho[0] + rho[1]);

    678:         /* Get eta values */
    679:         etaPoint[0].i   = ex;
    680:         etaPoint[0].j   = ey;
    681:         etaPoint[0].loc = DMSTAG_DOWN_LEFT;
    682:         etaPoint[0].c   = 0; /* Left  */
    683:         etaPoint[1].i   = ex;
    684:         etaPoint[1].j   = ey;
    685:         etaPoint[1].loc = DMSTAG_DOWN_RIGHT;
    686:         etaPoint[1].c   = 0; /* Right */
    687:         etaPoint[2].i   = ex;
    688:         etaPoint[2].j   = ey;
    689:         etaPoint[2].loc = DMSTAG_ELEMENT;
    690:         etaPoint[2].c   = 0; /* Up    */
    691:         etaPoint[3].i   = ex;
    692:         etaPoint[3].j   = ey - 1;
    693:         etaPoint[3].loc = DMSTAG_ELEMENT;
    694:         etaPoint[3].c   = 0; /* Down  */
    695:         PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 4, etaPoint, eta));
    696:         eta_left  = eta[0];
    697:         eta_right = eta[1];
    698:         eta_up    = eta[2];
    699:         eta_down  = eta[3];

    701:         count = 0;

    703:         col[count]   = row;
    704:         val_A[count] = -2.0 * dv * (eta_down + eta_up) / (hy * hy);
    705:         if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
    706:         if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
    707:         ++count;

    709:         col[count].i   = ex;
    710:         col[count].j   = ey - 1;
    711:         col[count].loc = DMSTAG_DOWN;
    712:         col[count].c   = 0;
    713:         val_A[count]   = 2.0 * dv * eta_down / (hy * hy);
    714:         ++count;
    715:         col[count].i   = ex;
    716:         col[count].j   = ey + 1;
    717:         col[count].loc = DMSTAG_DOWN;
    718:         col[count].c   = 0;
    719:         val_A[count]   = 2.0 * dv * eta_up / (hy * hy);
    720:         ++count;
    721:         if (!left_boundary) {
    722:           col[count].i   = ex - 1;
    723:           col[count].j   = ey;
    724:           col[count].loc = DMSTAG_DOWN;
    725:           col[count].c   = 0;
    726:           val_A[count]   = dv * eta_left / (hx * hx);
    727:           ++count;
    728:         }
    729:         if (!right_boundary) {
    730:           col[count].i   = ex + 1;
    731:           col[count].j   = ey;
    732:           col[count].loc = DMSTAG_DOWN;
    733:           col[count].c   = 0;
    734:           val_A[count]   = dv * eta_right / (hx * hx);
    735:           ++count;
    736:         }
    737:         col[count].i   = ex;
    738:         col[count].j   = ey - 1;
    739:         col[count].loc = DMSTAG_LEFT;
    740:         col[count].c   = 0;
    741:         val_A[count]   = dv * eta_left / (hx * hy);
    742:         ++count; /* down left x edge */
    743:         col[count].i   = ex;
    744:         col[count].j   = ey - 1;
    745:         col[count].loc = DMSTAG_RIGHT;
    746:         col[count].c   = 0;
    747:         val_A[count]   = -1.0 * dv * eta_right / (hx * hy);
    748:         ++count; /* down right x edge */
    749:         col[count].i   = ex;
    750:         col[count].j   = ey;
    751:         col[count].loc = DMSTAG_LEFT;
    752:         col[count].c   = 0;
    753:         val_A[count]   = -1.0 * dv * eta_left / (hx * hy);
    754:         ++count; /* up left x edge */
    755:         col[count].i   = ex;
    756:         col[count].j   = ey;
    757:         col[count].loc = DMSTAG_RIGHT;
    758:         col[count].c   = 0;
    759:         val_A[count]   = dv * eta_right / (hx * hy);
    760:         ++count; /* up right x edge */
    761:         if (!parameters->faces_only) {
    762:           col[count].i   = ex;
    763:           col[count].j   = ey - 1;
    764:           col[count].loc = DMSTAG_ELEMENT;
    765:           col[count].c   = 0;
    766:           val_A[count]   = K_cont * dv / hy;
    767:           ++count;
    768:           col[count].i   = ex;
    769:           col[count].j   = ey;
    770:           col[count].loc = DMSTAG_ELEMENT;
    771:           col[count].c   = 0;
    772:           val_A[count]   = -1.0 * K_cont * dv / hy;
    773:           ++count;
    774:         }

    776:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
    777:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    778:       }

    780:       if (ex == N[0] - 1) {
    781:         /* Right Boundary velocity Dirichlet */
    782:         /* Redundant in the corner */
    783:         DMStagStencil     row;
    784:         const PetscScalar val_rhs = 0.0;
    785:         const PetscScalar val_A   = K_bound;

    787:         row.i   = ex;
    788:         row.j   = ey;
    789:         row.loc = DMSTAG_RIGHT;
    790:         row.c   = 0;
    791:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    792:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    793:       }
    794:       if (ex == 0) {
    795:         /* Left velocity Dirichlet */
    796:         DMStagStencil     row;
    797:         const PetscScalar val_rhs = 0.0;
    798:         const PetscScalar val_A   = K_bound;

    800:         row.i   = ex;
    801:         row.j   = ey;
    802:         row.loc = DMSTAG_LEFT;
    803:         row.c   = 0;
    804:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    805:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    806:       } else {
    807:         /* X-momentum equation : (u_xx + u_yy) - p_x = f^x
    808:           zero RHS, including free-slip boundary conditions */
    809:         PetscInt          count;
    810:         DMStagStencil     row, col[11];
    811:         PetscScalar       val_A[11];
    812:         DMStagStencil     etaPoint[4];
    813:         PetscScalar       eta[4], eta_left, eta_right, eta_up, eta_down;
    814:         const PetscScalar val_rhs = 0.0;

    816:         row.i   = ex;
    817:         row.j   = ey;
    818:         row.loc = DMSTAG_LEFT;
    819:         row.c   = 0;

    821:         /* Get eta values */
    822:         etaPoint[0].i   = ex - 1;
    823:         etaPoint[0].j   = ey;
    824:         etaPoint[0].loc = DMSTAG_ELEMENT;
    825:         etaPoint[0].c   = 0; /* Left  */
    826:         etaPoint[1].i   = ex;
    827:         etaPoint[1].j   = ey;
    828:         etaPoint[1].loc = DMSTAG_ELEMENT;
    829:         etaPoint[1].c   = 0; /* Right */
    830:         etaPoint[2].i   = ex;
    831:         etaPoint[2].j   = ey;
    832:         etaPoint[2].loc = DMSTAG_UP_LEFT;
    833:         etaPoint[2].c   = 0; /* Up    */
    834:         etaPoint[3].i   = ex;
    835:         etaPoint[3].j   = ey;
    836:         etaPoint[3].loc = DMSTAG_DOWN_LEFT;
    837:         etaPoint[3].c   = 0; /* Down  */
    838:         PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 4, etaPoint, eta));
    839:         eta_left  = eta[0];
    840:         eta_right = eta[1];
    841:         eta_up    = eta[2];
    842:         eta_down  = eta[3];

    844:         count        = 0;
    845:         col[count]   = row;
    846:         val_A[count] = -2.0 * dv * (eta_left + eta_right) / (hx * hx);
    847:         if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
    848:         if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
    849:         ++count;

    851:         if (!bottom_boundary) {
    852:           col[count].i   = ex;
    853:           col[count].j   = ey - 1;
    854:           col[count].loc = DMSTAG_LEFT;
    855:           col[count].c   = 0;
    856:           val_A[count]   = dv * eta_down / (hy * hy);
    857:           ++count;
    858:         }
    859:         if (!top_boundary) {
    860:           col[count].i   = ex;
    861:           col[count].j   = ey + 1;
    862:           col[count].loc = DMSTAG_LEFT;
    863:           col[count].c   = 0;
    864:           val_A[count]   = dv * eta_up / (hy * hy);
    865:           ++count;
    866:         }
    867:         col[count].i   = ex - 1;
    868:         col[count].j   = ey;
    869:         col[count].loc = DMSTAG_LEFT;
    870:         col[count].c   = 0;
    871:         val_A[count]   = 2.0 * dv * eta_left / (hx * hx);
    872:         ++count;
    873:         col[count].i   = ex + 1;
    874:         col[count].j   = ey;
    875:         col[count].loc = DMSTAG_LEFT;
    876:         col[count].c   = 0;
    877:         val_A[count]   = 2.0 * dv * eta_right / (hx * hx);
    878:         ++count;
    879:         col[count].i   = ex - 1;
    880:         col[count].j   = ey;
    881:         col[count].loc = DMSTAG_DOWN;
    882:         col[count].c   = 0;
    883:         val_A[count]   = dv * eta_down / (hx * hy);
    884:         ++count; /* down left */
    885:         col[count].i   = ex;
    886:         col[count].j   = ey;
    887:         col[count].loc = DMSTAG_DOWN;
    888:         col[count].c   = 0;
    889:         val_A[count]   = -1.0 * dv * eta_down / (hx * hy);
    890:         ++count; /* down right */
    891:         col[count].i   = ex - 1;
    892:         col[count].j   = ey;
    893:         col[count].loc = DMSTAG_UP;
    894:         col[count].c   = 0;
    895:         val_A[count]   = -1.0 * dv * eta_up / (hx * hy);
    896:         ++count; /* up left */
    897:         col[count].i   = ex;
    898:         col[count].j   = ey;
    899:         col[count].loc = DMSTAG_UP;
    900:         col[count].c   = 0;
    901:         val_A[count]   = dv * eta_up / (hx * hy);
    902:         ++count; /* up right */
    903:         if (!parameters->faces_only) {
    904:           col[count].i   = ex - 1;
    905:           col[count].j   = ey;
    906:           col[count].loc = DMSTAG_ELEMENT;
    907:           col[count].c   = 0;
    908:           val_A[count]   = K_cont * dv / hx;
    909:           ++count;
    910:           col[count].i   = ex;
    911:           col[count].j   = ey;
    912:           col[count].loc = DMSTAG_ELEMENT;
    913:           col[count].c   = 0;
    914:           val_A[count]   = -1.0 * K_cont * dv / hx;
    915:           ++count;
    916:         }

    918:         PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
    919:         if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    920:       }

    922:       /* P equation : u_x + v_y = 0

    924:          Note that this includes an explicit zero on the diagonal. This is only needed for
    925:          direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace)

    927:         Note: the scaling by dv is not chosen in a principled way and is likely sub-optimal
    928:        */
    929:       if (!parameters->faces_only) {
    930:         if (ctx->pin_pressure && ex == 0 && ey == 0) { /* Pin the first pressure node to zero, if requested */
    931:           DMStagStencil     row;
    932:           const PetscScalar val_A   = K_bound;
    933:           const PetscScalar val_rhs = 0.0;

    935:           row.i   = ex;
    936:           row.j   = ey;
    937:           row.loc = DMSTAG_ELEMENT;
    938:           row.c   = 0;
    939:           PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    940:           if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    941:         } else {
    942:           DMStagStencil     row, col[5];
    943:           PetscScalar       val_A[5];
    944:           const PetscScalar val_rhs = 0.0;

    946:           row.i      = ex;
    947:           row.j      = ey;
    948:           row.loc    = DMSTAG_ELEMENT;
    949:           row.c      = 0;
    950:           col[0].i   = ex;
    951:           col[0].j   = ey;
    952:           col[0].loc = DMSTAG_LEFT;
    953:           col[0].c   = 0;
    954:           val_A[0]   = -1.0 * K_cont * dv / hx;
    955:           col[1].i   = ex;
    956:           col[1].j   = ey;
    957:           col[1].loc = DMSTAG_RIGHT;
    958:           col[1].c   = 0;
    959:           val_A[1]   = K_cont * dv / hx;
    960:           col[2].i   = ex;
    961:           col[2].j   = ey;
    962:           col[2].loc = DMSTAG_DOWN;
    963:           col[2].c   = 0;
    964:           val_A[2]   = -1.0 * K_cont * dv / hy;
    965:           col[3].i   = ex;
    966:           col[3].j   = ey;
    967:           col[3].loc = DMSTAG_UP;
    968:           col[3].c   = 0;
    969:           val_A[3]   = K_cont * dv / hy;
    970:           col[4]     = row;
    971:           val_A[4]   = 0.0;
    972:           PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 5, col, val_A, INSERT_VALUES));
    973:           if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    974:         }
    975:       }
    976:     }
    977:   }
    978:   PetscCall(DMRestoreLocalVector(dm_coefficients, &coefficients_local));

    980:   /* Add additional inverse viscosity terms (for use in building a preconditioning matrix) */
    981:   if (parameters->include_inverse_visc) {
    982:     PetscCheck(!parameters->faces_only, PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "Does not make sense with faces only");
    983:     PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_main, dm_coefficients, ctx->levels[level]->coeff, 1.0, A));
    984:   }

    986:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    987:   if (build_rhs) PetscCall(VecAssemblyBegin(rhs));
    988:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    989:   if (build_rhs) PetscCall(VecAssemblyEnd(rhs));
    990:   PetscFunctionReturn(PETSC_SUCCESS);
    991: }

    993: static PetscErrorCode CreateSystem3d(SystemParameters parameters, Mat *pA, Vec *pRhs)
    994: {
    995:   PetscInt    N[3];
    996:   PetscInt    ex, ey, ez, startx, starty, startz, nx, ny, nz;
    997:   Mat         A;
    998:   PetscReal   hx, hy, hz, dv;
    999:   PetscInt    pinx, piny, pinz;
    1000:   Vec         coeff_local, rhs;
    1001:   PetscBool   build_rhs;
    1002:   DM          dm_main, dm_coefficients;
    1003:   PetscScalar K_cont, K_bound;
    1004:   Ctx         ctx   = parameters->ctx;
    1005:   PetscInt    level = parameters->level;

    1007:   PetscFunctionBeginUser;
    1008:   if (parameters->faces_only) {
    1009:     dm_main = ctx->levels[level]->dm_faces;
    1010:   } else {
    1011:     dm_main = ctx->levels[level]->dm_stokes;
    1012:   }
    1013:   dm_coefficients = ctx->levels[level]->dm_coefficients;
    1014:   K_cont          = ctx->levels[level]->K_cont;
    1015:   K_bound         = ctx->levels[level]->K_bound;
    1016:   PetscCall(DMCreateMatrix(dm_main, pA));
    1017:   A         = *pA;
    1018:   build_rhs = (PetscBool)(pRhs != NULL);
    1019:   if (build_rhs) {
    1020:     PetscCall(DMCreateGlobalVector(dm_main, pRhs));
    1021:     rhs = *pRhs;
    1022:   } else {
    1023:     rhs = NULL;
    1024:   }
    1025:   PetscCall(DMStagGetCorners(dm_main, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
    1026:   PetscCall(DMStagGetGlobalSizes(dm_main, &N[0], &N[1], &N[2]));
    1027:   hx = ctx->levels[level]->hx_characteristic;
    1028:   hy = ctx->levels[level]->hy_characteristic;
    1029:   hz = ctx->levels[level]->hz_characteristic;
    1030:   dv = hx * hy * hz;
    1031:   PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
    1032:   PetscCall(DMGlobalToLocal(dm_coefficients, ctx->levels[level]->coeff, INSERT_VALUES, coeff_local));

    1034:   PetscCheck(N[0] >= 2 && N[1] >= 2 && N[2] >= 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for less than 2 elements in any direction");
    1035:   pinx = 1;
    1036:   piny = 0;
    1037:   pinz = 0; /* Depends on assertion above that there are at least two element in the x direction */

    1039:   /* Loop over all local elements.

    1041:      For each element, fill 4-7 rows of the matrix, corresponding to
    1042:      - the pressure degree of freedom (dof), centered on the element
    1043:      - the 3 velocity dofs on left, bottom, and back faces of the element
    1044:      - velocity dof on the right, top, and front faces of the element (only on domain boundaries)

    1046:    */
    1047:   for (ez = startz; ez < startz + nz; ++ez) {
    1048:     for (ey = starty; ey < starty + ny; ++ey) {
    1049:       for (ex = startx; ex < startx + nx; ++ex) {
    1050:         const PetscBool left_boundary   = (PetscBool)(ex == 0);
    1051:         const PetscBool right_boundary  = (PetscBool)(ex == N[0] - 1);
    1052:         const PetscBool bottom_boundary = (PetscBool)(ey == 0);
    1053:         const PetscBool top_boundary    = (PetscBool)(ey == N[1] - 1);
    1054:         const PetscBool back_boundary   = (PetscBool)(ez == 0);
    1055:         const PetscBool front_boundary  = (PetscBool)(ez == N[2] - 1);

    1057:         /* Note that below, we depend on the check above that there is never one
    1058:            element (globally) in a given direction.  Thus, for example, an
    1059:            element is never both on the left and right boundary */

    1061:         /* X-faces - right boundary */
    1062:         if (right_boundary) {
    1063:           /* Right x-velocity Dirichlet */
    1064:           DMStagStencil     row;
    1065:           const PetscScalar val_rhs = 0.0;
    1066:           const PetscScalar val_A   = K_bound;

    1068:           row.i   = ex;
    1069:           row.j   = ey;
    1070:           row.k   = ez;
    1071:           row.loc = DMSTAG_RIGHT;
    1072:           row.c   = 0;
    1073:           PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1074:           if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1075:         }

    1077:         /* X faces - left*/
    1078:         {
    1079:           DMStagStencil row;

    1081:           row.i   = ex;
    1082:           row.j   = ey;
    1083:           row.k   = ez;
    1084:           row.loc = DMSTAG_LEFT;
    1085:           row.c   = 0;

    1087:           if (left_boundary) {
    1088:             /* Left x-velocity Dirichlet */
    1089:             const PetscScalar val_rhs = 0.0;
    1090:             const PetscScalar val_A   = K_bound;

    1092:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1093:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1094:           } else {
    1095:             /* X-momentum equation */
    1096:             PetscInt          count;
    1097:             DMStagStencil     col[17];
    1098:             PetscScalar       val_A[17];
    1099:             DMStagStencil     eta_point[6];
    1100:             PetscScalar       eta[6], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the left face */
    1101:             const PetscScalar val_rhs = 0.0;

    1103:             /* Get eta values */
    1104:             eta_point[0].i   = ex - 1;
    1105:             eta_point[0].j   = ey;
    1106:             eta_point[0].k   = ez;
    1107:             eta_point[0].loc = DMSTAG_ELEMENT;
    1108:             eta_point[0].c   = 0; /* Left  */
    1109:             eta_point[1].i   = ex;
    1110:             eta_point[1].j   = ey;
    1111:             eta_point[1].k   = ez;
    1112:             eta_point[1].loc = DMSTAG_ELEMENT;
    1113:             eta_point[1].c   = 0; /* Right */
    1114:             eta_point[2].i   = ex;
    1115:             eta_point[2].j   = ey;
    1116:             eta_point[2].k   = ez;
    1117:             eta_point[2].loc = DMSTAG_DOWN_LEFT;
    1118:             eta_point[2].c   = 0; /* Down  */
    1119:             eta_point[3].i   = ex;
    1120:             eta_point[3].j   = ey;
    1121:             eta_point[3].k   = ez;
    1122:             eta_point[3].loc = DMSTAG_UP_LEFT;
    1123:             eta_point[3].c   = 0; /* Up    */
    1124:             eta_point[4].i   = ex;
    1125:             eta_point[4].j   = ey;
    1126:             eta_point[4].k   = ez;
    1127:             eta_point[4].loc = DMSTAG_BACK_LEFT;
    1128:             eta_point[4].c   = 0; /* Back  */
    1129:             eta_point[5].i   = ex;
    1130:             eta_point[5].j   = ey;
    1131:             eta_point[5].k   = ez;
    1132:             eta_point[5].loc = DMSTAG_FRONT_LEFT;
    1133:             eta_point[5].c   = 0; /* Front  */
    1134:             PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
    1135:             eta_left  = eta[0];
    1136:             eta_right = eta[1];
    1137:             eta_down  = eta[2];
    1138:             eta_up    = eta[3];
    1139:             eta_back  = eta[4];
    1140:             eta_front = eta[5];

    1142:             count = 0;

    1144:             col[count]   = row;
    1145:             val_A[count] = -2.0 * dv * (eta_left + eta_right) / (hx * hx);
    1146:             if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
    1147:             if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
    1148:             if (!back_boundary) val_A[count] += -1.0 * dv * eta_back / (hz * hz);
    1149:             if (!front_boundary) val_A[count] += -1.0 * dv * eta_front / (hz * hz);
    1150:             ++count;

    1152:             col[count].i   = ex - 1;
    1153:             col[count].j   = ey;
    1154:             col[count].k   = ez;
    1155:             col[count].loc = DMSTAG_LEFT;
    1156:             col[count].c   = 0;
    1157:             val_A[count]   = 2.0 * dv * eta_left / (hx * hx);
    1158:             ++count;
    1159:             col[count].i   = ex + 1;
    1160:             col[count].j   = ey;
    1161:             col[count].k   = ez;
    1162:             col[count].loc = DMSTAG_LEFT;
    1163:             col[count].c   = 0;
    1164:             val_A[count]   = 2.0 * dv * eta_right / (hx * hx);
    1165:             ++count;
    1166:             if (!bottom_boundary) {
    1167:               col[count].i   = ex;
    1168:               col[count].j   = ey - 1;
    1169:               col[count].k   = ez;
    1170:               col[count].loc = DMSTAG_LEFT;
    1171:               col[count].c   = 0;
    1172:               val_A[count]   = dv * eta_down / (hy * hy);
    1173:               ++count;
    1174:             }
    1175:             if (!top_boundary) {
    1176:               col[count].i   = ex;
    1177:               col[count].j   = ey + 1;
    1178:               col[count].k   = ez;
    1179:               col[count].loc = DMSTAG_LEFT;
    1180:               col[count].c   = 0;
    1181:               val_A[count]   = dv * eta_up / (hy * hy);
    1182:               ++count;
    1183:             }
    1184:             if (!back_boundary) {
    1185:               col[count].i   = ex;
    1186:               col[count].j   = ey;
    1187:               col[count].k   = ez - 1;
    1188:               col[count].loc = DMSTAG_LEFT;
    1189:               col[count].c   = 0;
    1190:               val_A[count]   = dv * eta_back / (hz * hz);
    1191:               ++count;
    1192:             }
    1193:             if (!front_boundary) {
    1194:               col[count].i   = ex;
    1195:               col[count].j   = ey;
    1196:               col[count].k   = ez + 1;
    1197:               col[count].loc = DMSTAG_LEFT;
    1198:               col[count].c   = 0;
    1199:               val_A[count]   = dv * eta_front / (hz * hz);
    1200:               ++count;
    1201:             }

    1203:             col[count].i   = ex - 1;
    1204:             col[count].j   = ey;
    1205:             col[count].k   = ez;
    1206:             col[count].loc = DMSTAG_DOWN;
    1207:             col[count].c   = 0;
    1208:             val_A[count]   = dv * eta_down / (hx * hy);
    1209:             ++count; /* down left */
    1210:             col[count].i   = ex;
    1211:             col[count].j   = ey;
    1212:             col[count].k   = ez;
    1213:             col[count].loc = DMSTAG_DOWN;
    1214:             col[count].c   = 0;
    1215:             val_A[count]   = -1.0 * dv * eta_down / (hx * hy);
    1216:             ++count; /* down right */

    1218:             col[count].i   = ex - 1;
    1219:             col[count].j   = ey;
    1220:             col[count].k   = ez;
    1221:             col[count].loc = DMSTAG_UP;
    1222:             col[count].c   = 0;
    1223:             val_A[count]   = -1.0 * dv * eta_up / (hx * hy);
    1224:             ++count; /* up left */
    1225:             col[count].i   = ex;
    1226:             col[count].j   = ey;
    1227:             col[count].k   = ez;
    1228:             col[count].loc = DMSTAG_UP;
    1229:             col[count].c   = 0;
    1230:             val_A[count]   = dv * eta_up / (hx * hy);
    1231:             ++count; /* up right */

    1233:             col[count].i   = ex - 1;
    1234:             col[count].j   = ey;
    1235:             col[count].k   = ez;
    1236:             col[count].loc = DMSTAG_BACK;
    1237:             col[count].c   = 0;
    1238:             val_A[count]   = dv * eta_back / (hx * hz);
    1239:             ++count; /* back left */
    1240:             col[count].i   = ex;
    1241:             col[count].j   = ey;
    1242:             col[count].k   = ez;
    1243:             col[count].loc = DMSTAG_BACK;
    1244:             col[count].c   = 0;
    1245:             val_A[count]   = -1.0 * dv * eta_back / (hx * hz);
    1246:             ++count; /* back right */

    1248:             col[count].i   = ex - 1;
    1249:             col[count].j   = ey;
    1250:             col[count].k   = ez;
    1251:             col[count].loc = DMSTAG_FRONT;
    1252:             col[count].c   = 0;
    1253:             val_A[count]   = -1.0 * dv * eta_front / (hx * hz);
    1254:             ++count; /* front left */
    1255:             col[count].i   = ex;
    1256:             col[count].j   = ey;
    1257:             col[count].k   = ez;
    1258:             col[count].loc = DMSTAG_FRONT;
    1259:             col[count].c   = 0;
    1260:             val_A[count]   = dv * eta_front / (hx * hz);
    1261:             ++count; /* front right */

    1263:             if (!parameters->faces_only) {
    1264:               col[count].i   = ex - 1;
    1265:               col[count].j   = ey;
    1266:               col[count].k   = ez;
    1267:               col[count].loc = DMSTAG_ELEMENT;
    1268:               col[count].c   = 0;
    1269:               val_A[count]   = K_cont * dv / hx;
    1270:               ++count;
    1271:               col[count].i   = ex;
    1272:               col[count].j   = ey;
    1273:               col[count].k   = ez;
    1274:               col[count].loc = DMSTAG_ELEMENT;
    1275:               col[count].c   = 0;
    1276:               val_A[count]   = -1.0 * K_cont * dv / hx;
    1277:               ++count;
    1278:             }

    1280:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
    1281:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1282:           }
    1283:         }

    1285:         /* Y faces - top boundary */
    1286:         if (top_boundary) {
    1287:           /* Top y-velocity Dirichlet */
    1288:           DMStagStencil     row;
    1289:           const PetscScalar val_rhs = 0.0;
    1290:           const PetscScalar val_A   = K_bound;

    1292:           row.i   = ex;
    1293:           row.j   = ey;
    1294:           row.k   = ez;
    1295:           row.loc = DMSTAG_UP;
    1296:           row.c   = 0;
    1297:           PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1298:           if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1299:         }

    1301:         /* Y faces - down */
    1302:         {
    1303:           DMStagStencil row;

    1305:           row.i   = ex;
    1306:           row.j   = ey;
    1307:           row.k   = ez;
    1308:           row.loc = DMSTAG_DOWN;
    1309:           row.c   = 0;

    1311:           if (bottom_boundary) {
    1312:             /* Bottom y-velocity Dirichlet */
    1313:             const PetscScalar val_rhs = 0.0;
    1314:             const PetscScalar val_A   = K_bound;

    1316:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1317:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1318:           } else {
    1319:             /* Y-momentum equation (including non-zero forcing) */
    1320:             PetscInt      count;
    1321:             DMStagStencil col[17];
    1322:             PetscScalar   val_rhs, val_A[17];
    1323:             DMStagStencil eta_point[6], rho_point[4];
    1324:             PetscScalar   eta[6], rho[4], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the bottom face */

    1326:             if (build_rhs) {
    1327:               /* get rho values  (note .c = 1) */
    1328:               /* Note that we have rho at perhaps strange points (edges not corners) */
    1329:               rho_point[0].i   = ex;
    1330:               rho_point[0].j   = ey;
    1331:               rho_point[0].k   = ez;
    1332:               rho_point[0].loc = DMSTAG_DOWN_LEFT;
    1333:               rho_point[0].c   = 1;
    1334:               rho_point[1].i   = ex;
    1335:               rho_point[1].j   = ey;
    1336:               rho_point[1].k   = ez;
    1337:               rho_point[1].loc = DMSTAG_DOWN_RIGHT;
    1338:               rho_point[1].c   = 1;
    1339:               rho_point[2].i   = ex;
    1340:               rho_point[2].j   = ey;
    1341:               rho_point[2].k   = ez;
    1342:               rho_point[2].loc = DMSTAG_BACK_DOWN;
    1343:               rho_point[2].c   = 1;
    1344:               rho_point[3].i   = ex;
    1345:               rho_point[3].j   = ey;
    1346:               rho_point[3].k   = ez;
    1347:               rho_point[3].loc = DMSTAG_FRONT_DOWN;
    1348:               rho_point[3].c   = 1;
    1349:               PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 4, rho_point, rho));

    1351:               /* Compute forcing */
    1352:               val_rhs = ctx->gy * dv * (0.25 * (rho[0] + rho[1] + rho[2] + rho[3]));
    1353:             }

    1355:             /* Get eta values */
    1356:             eta_point[0].i   = ex;
    1357:             eta_point[0].j   = ey;
    1358:             eta_point[0].k   = ez;
    1359:             eta_point[0].loc = DMSTAG_DOWN_LEFT;
    1360:             eta_point[0].c   = 0; /* Left  */
    1361:             eta_point[1].i   = ex;
    1362:             eta_point[1].j   = ey;
    1363:             eta_point[1].k   = ez;
    1364:             eta_point[1].loc = DMSTAG_DOWN_RIGHT;
    1365:             eta_point[1].c   = 0; /* Right */
    1366:             eta_point[2].i   = ex;
    1367:             eta_point[2].j   = ey - 1;
    1368:             eta_point[2].k   = ez;
    1369:             eta_point[2].loc = DMSTAG_ELEMENT;
    1370:             eta_point[2].c   = 0; /* Down  */
    1371:             eta_point[3].i   = ex;
    1372:             eta_point[3].j   = ey;
    1373:             eta_point[3].k   = ez;
    1374:             eta_point[3].loc = DMSTAG_ELEMENT;
    1375:             eta_point[3].c   = 0; /* Up    */
    1376:             eta_point[4].i   = ex;
    1377:             eta_point[4].j   = ey;
    1378:             eta_point[4].k   = ez;
    1379:             eta_point[4].loc = DMSTAG_BACK_DOWN;
    1380:             eta_point[4].c   = 0; /* Back  */
    1381:             eta_point[5].i   = ex;
    1382:             eta_point[5].j   = ey;
    1383:             eta_point[5].k   = ez;
    1384:             eta_point[5].loc = DMSTAG_FRONT_DOWN;
    1385:             eta_point[5].c   = 0; /* Front  */
    1386:             PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
    1387:             eta_left  = eta[0];
    1388:             eta_right = eta[1];
    1389:             eta_down  = eta[2];
    1390:             eta_up    = eta[3];
    1391:             eta_back  = eta[4];
    1392:             eta_front = eta[5];

    1394:             count = 0;

    1396:             col[count]   = row;
    1397:             val_A[count] = -2.0 * dv * (eta_up + eta_down) / (hy * hy);
    1398:             if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
    1399:             if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
    1400:             if (!back_boundary) val_A[count] += -1.0 * dv * eta_back / (hz * hz);
    1401:             if (!front_boundary) val_A[count] += -1.0 * dv * eta_front / (hz * hz);
    1402:             ++count;

    1404:             col[count].i   = ex;
    1405:             col[count].j   = ey - 1;
    1406:             col[count].k   = ez;
    1407:             col[count].loc = DMSTAG_DOWN;
    1408:             col[count].c   = 0;
    1409:             val_A[count]   = 2.0 * dv * eta_down / (hy * hy);
    1410:             ++count;
    1411:             col[count].i   = ex;
    1412:             col[count].j   = ey + 1;
    1413:             col[count].k   = ez;
    1414:             col[count].loc = DMSTAG_DOWN;
    1415:             col[count].c   = 0;
    1416:             val_A[count]   = 2.0 * dv * eta_up / (hy * hy);
    1417:             ++count;

    1419:             if (!left_boundary) {
    1420:               col[count].i   = ex - 1;
    1421:               col[count].j   = ey;
    1422:               col[count].k   = ez;
    1423:               col[count].loc = DMSTAG_DOWN;
    1424:               col[count].c   = 0;
    1425:               val_A[count]   = dv * eta_left / (hx * hx);
    1426:               ++count;
    1427:             }
    1428:             if (!right_boundary) {
    1429:               col[count].i   = ex + 1;
    1430:               col[count].j   = ey;
    1431:               col[count].k   = ez;
    1432:               col[count].loc = DMSTAG_DOWN;
    1433:               col[count].c   = 0;
    1434:               val_A[count]   = dv * eta_right / (hx * hx);
    1435:               ++count;
    1436:             }
    1437:             if (!back_boundary) {
    1438:               col[count].i   = ex;
    1439:               col[count].j   = ey;
    1440:               col[count].k   = ez - 1;
    1441:               col[count].loc = DMSTAG_DOWN;
    1442:               col[count].c   = 0;
    1443:               val_A[count]   = dv * eta_back / (hz * hz);
    1444:               ++count;
    1445:             }
    1446:             if (!front_boundary) {
    1447:               col[count].i   = ex;
    1448:               col[count].j   = ey;
    1449:               col[count].k   = ez + 1;
    1450:               col[count].loc = DMSTAG_DOWN;
    1451:               col[count].c   = 0;
    1452:               val_A[count]   = dv * eta_front / (hz * hz);
    1453:               ++count;
    1454:             }

    1456:             col[count].i   = ex;
    1457:             col[count].j   = ey - 1;
    1458:             col[count].k   = ez;
    1459:             col[count].loc = DMSTAG_LEFT;
    1460:             col[count].c   = 0;
    1461:             val_A[count]   = dv * eta_left / (hx * hy);
    1462:             ++count; /* down left*/
    1463:             col[count].i   = ex;
    1464:             col[count].j   = ey;
    1465:             col[count].k   = ez;
    1466:             col[count].loc = DMSTAG_LEFT;
    1467:             col[count].c   = 0;
    1468:             val_A[count]   = -1.0 * dv * eta_left / (hx * hy);
    1469:             ++count; /* up left*/

    1471:             col[count].i   = ex;
    1472:             col[count].j   = ey - 1;
    1473:             col[count].k   = ez;
    1474:             col[count].loc = DMSTAG_RIGHT;
    1475:             col[count].c   = 0;
    1476:             val_A[count]   = -1.0 * dv * eta_right / (hx * hy);
    1477:             ++count; /* down right*/
    1478:             col[count].i   = ex;
    1479:             col[count].j   = ey;
    1480:             col[count].k   = ez;
    1481:             col[count].loc = DMSTAG_RIGHT;
    1482:             col[count].c   = 0;
    1483:             val_A[count]   = dv * eta_right / (hx * hy);
    1484:             ++count; /* up right*/

    1486:             col[count].i   = ex;
    1487:             col[count].j   = ey - 1;
    1488:             col[count].k   = ez;
    1489:             col[count].loc = DMSTAG_BACK;
    1490:             col[count].c   = 0;
    1491:             val_A[count]   = dv * eta_back / (hy * hz);
    1492:             ++count; /* back down */
    1493:             col[count].i   = ex;
    1494:             col[count].j   = ey;
    1495:             col[count].k   = ez;
    1496:             col[count].loc = DMSTAG_BACK;
    1497:             col[count].c   = 0;
    1498:             val_A[count]   = -1.0 * dv * eta_back / (hy * hz);
    1499:             ++count; /* back up */

    1501:             col[count].i   = ex;
    1502:             col[count].j   = ey - 1;
    1503:             col[count].k   = ez;
    1504:             col[count].loc = DMSTAG_FRONT;
    1505:             col[count].c   = 0;
    1506:             val_A[count]   = -1.0 * dv * eta_front / (hy * hz);
    1507:             ++count; /* front down */
    1508:             col[count].i   = ex;
    1509:             col[count].j   = ey;
    1510:             col[count].k   = ez;
    1511:             col[count].loc = DMSTAG_FRONT;
    1512:             col[count].c   = 0;
    1513:             val_A[count]   = dv * eta_front / (hy * hz);
    1514:             ++count; /* front up */

    1516:             if (!parameters->faces_only) {
    1517:               col[count].i   = ex;
    1518:               col[count].j   = ey - 1;
    1519:               col[count].k   = ez;
    1520:               col[count].loc = DMSTAG_ELEMENT;
    1521:               col[count].c   = 0;
    1522:               val_A[count]   = K_cont * dv / hy;
    1523:               ++count;
    1524:               col[count].i   = ex;
    1525:               col[count].j   = ey;
    1526:               col[count].k   = ez;
    1527:               col[count].loc = DMSTAG_ELEMENT;
    1528:               col[count].c   = 0;
    1529:               val_A[count]   = -1.0 * K_cont * dv / hy;
    1530:               ++count;
    1531:             }

    1533:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
    1534:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1535:           }
    1536:         }

    1538:         if (front_boundary) {
    1539:           /* Front z-velocity Dirichlet */
    1540:           DMStagStencil     row;
    1541:           const PetscScalar val_rhs = 0.0;
    1542:           const PetscScalar val_A   = K_bound;

    1544:           row.i   = ex;
    1545:           row.j   = ey;
    1546:           row.k   = ez;
    1547:           row.loc = DMSTAG_FRONT;
    1548:           row.c   = 0;
    1549:           PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1550:           if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1551:         }

    1553:         /* Z faces - back */
    1554:         {
    1555:           DMStagStencil row;

    1557:           row.i   = ex;
    1558:           row.j   = ey;
    1559:           row.k   = ez;
    1560:           row.loc = DMSTAG_BACK;
    1561:           row.c   = 0;

    1563:           if (back_boundary) {
    1564:             /* Back z-velocity Dirichlet */
    1565:             const PetscScalar val_rhs = 0.0;
    1566:             const PetscScalar val_A   = K_bound;

    1568:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1569:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1570:           } else {
    1571:             /* Z-momentum equation */
    1572:             PetscInt          count;
    1573:             DMStagStencil     col[17];
    1574:             PetscScalar       val_A[17];
    1575:             DMStagStencil     eta_point[6];
    1576:             PetscScalar       eta[6], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the back face */
    1577:             const PetscScalar val_rhs = 0.0;

    1579:             /* Get eta values */
    1580:             eta_point[0].i   = ex;
    1581:             eta_point[0].j   = ey;
    1582:             eta_point[0].k   = ez;
    1583:             eta_point[0].loc = DMSTAG_BACK_LEFT;
    1584:             eta_point[0].c   = 0; /* Left  */
    1585:             eta_point[1].i   = ex;
    1586:             eta_point[1].j   = ey;
    1587:             eta_point[1].k   = ez;
    1588:             eta_point[1].loc = DMSTAG_BACK_RIGHT;
    1589:             eta_point[1].c   = 0; /* Right */
    1590:             eta_point[2].i   = ex;
    1591:             eta_point[2].j   = ey;
    1592:             eta_point[2].k   = ez;
    1593:             eta_point[2].loc = DMSTAG_BACK_DOWN;
    1594:             eta_point[2].c   = 0; /* Down  */
    1595:             eta_point[3].i   = ex;
    1596:             eta_point[3].j   = ey;
    1597:             eta_point[3].k   = ez;
    1598:             eta_point[3].loc = DMSTAG_BACK_UP;
    1599:             eta_point[3].c   = 0; /* Up    */
    1600:             eta_point[4].i   = ex;
    1601:             eta_point[4].j   = ey;
    1602:             eta_point[4].k   = ez - 1;
    1603:             eta_point[4].loc = DMSTAG_ELEMENT;
    1604:             eta_point[4].c   = 0; /* Back  */
    1605:             eta_point[5].i   = ex;
    1606:             eta_point[5].j   = ey;
    1607:             eta_point[5].k   = ez;
    1608:             eta_point[5].loc = DMSTAG_ELEMENT;
    1609:             eta_point[5].c   = 0; /* Front  */
    1610:             PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
    1611:             eta_left  = eta[0];
    1612:             eta_right = eta[1];
    1613:             eta_down  = eta[2];
    1614:             eta_up    = eta[3];
    1615:             eta_back  = eta[4];
    1616:             eta_front = eta[5];

    1618:             count = 0;

    1620:             col[count]   = row;
    1621:             val_A[count] = -2.0 * dv * (eta_back + eta_front) / (hz * hz);
    1622:             if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
    1623:             if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
    1624:             if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
    1625:             if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
    1626:             ++count;

    1628:             col[count].i   = ex;
    1629:             col[count].j   = ey;
    1630:             col[count].k   = ez - 1;
    1631:             col[count].loc = DMSTAG_BACK;
    1632:             col[count].c   = 0;
    1633:             val_A[count]   = 2.0 * dv * eta_back / (hz * hz);
    1634:             ++count;
    1635:             col[count].i   = ex;
    1636:             col[count].j   = ey;
    1637:             col[count].k   = ez + 1;
    1638:             col[count].loc = DMSTAG_BACK;
    1639:             col[count].c   = 0;
    1640:             val_A[count]   = 2.0 * dv * eta_front / (hz * hz);
    1641:             ++count;

    1643:             if (!left_boundary) {
    1644:               col[count].i   = ex - 1;
    1645:               col[count].j   = ey;
    1646:               col[count].k   = ez;
    1647:               col[count].loc = DMSTAG_BACK;
    1648:               col[count].c   = 0;
    1649:               val_A[count]   = dv * eta_left / (hx * hx);
    1650:               ++count;
    1651:             }
    1652:             if (!right_boundary) {
    1653:               col[count].i   = ex + 1;
    1654:               col[count].j   = ey;
    1655:               col[count].k   = ez;
    1656:               col[count].loc = DMSTAG_BACK;
    1657:               col[count].c   = 0;
    1658:               val_A[count]   = dv * eta_right / (hx * hx);
    1659:               ++count;
    1660:             }
    1661:             if (!bottom_boundary) {
    1662:               col[count].i   = ex;
    1663:               col[count].j   = ey - 1;
    1664:               col[count].k   = ez;
    1665:               col[count].loc = DMSTAG_BACK;
    1666:               col[count].c   = 0;
    1667:               val_A[count]   = dv * eta_down / (hy * hy);
    1668:               ++count;
    1669:             }
    1670:             if (!top_boundary) {
    1671:               col[count].i   = ex;
    1672:               col[count].j   = ey + 1;
    1673:               col[count].k   = ez;
    1674:               col[count].loc = DMSTAG_BACK;
    1675:               col[count].c   = 0;
    1676:               val_A[count]   = dv * eta_up / (hy * hy);
    1677:               ++count;
    1678:             }

    1680:             col[count].i   = ex;
    1681:             col[count].j   = ey;
    1682:             col[count].k   = ez - 1;
    1683:             col[count].loc = DMSTAG_LEFT;
    1684:             col[count].c   = 0;
    1685:             val_A[count]   = dv * eta_left / (hx * hz);
    1686:             ++count; /* back left*/
    1687:             col[count].i   = ex;
    1688:             col[count].j   = ey;
    1689:             col[count].k   = ez;
    1690:             col[count].loc = DMSTAG_LEFT;
    1691:             col[count].c   = 0;
    1692:             val_A[count]   = -1.0 * dv * eta_left / (hx * hz);
    1693:             ++count; /* front left*/

    1695:             col[count].i   = ex;
    1696:             col[count].j   = ey;
    1697:             col[count].k   = ez - 1;
    1698:             col[count].loc = DMSTAG_RIGHT;
    1699:             col[count].c   = 0;
    1700:             val_A[count]   = -1.0 * dv * eta_right / (hx * hz);
    1701:             ++count; /* back right */
    1702:             col[count].i   = ex;
    1703:             col[count].j   = ey;
    1704:             col[count].k   = ez;
    1705:             col[count].loc = DMSTAG_RIGHT;
    1706:             col[count].c   = 0;
    1707:             val_A[count]   = dv * eta_right / (hx * hz);
    1708:             ++count; /* front right*/

    1710:             col[count].i   = ex;
    1711:             col[count].j   = ey;
    1712:             col[count].k   = ez - 1;
    1713:             col[count].loc = DMSTAG_DOWN;
    1714:             col[count].c   = 0;
    1715:             val_A[count]   = dv * eta_down / (hy * hz);
    1716:             ++count; /* back down */
    1717:             col[count].i   = ex;
    1718:             col[count].j   = ey;
    1719:             col[count].k   = ez;
    1720:             col[count].loc = DMSTAG_DOWN;
    1721:             col[count].c   = 0;
    1722:             val_A[count]   = -1.0 * dv * eta_down / (hy * hz);
    1723:             ++count; /* back down */

    1725:             col[count].i   = ex;
    1726:             col[count].j   = ey;
    1727:             col[count].k   = ez - 1;
    1728:             col[count].loc = DMSTAG_UP;
    1729:             col[count].c   = 0;
    1730:             val_A[count]   = -1.0 * dv * eta_up / (hy * hz);
    1731:             ++count; /* back up */
    1732:             col[count].i   = ex;
    1733:             col[count].j   = ey;
    1734:             col[count].k   = ez;
    1735:             col[count].loc = DMSTAG_UP;
    1736:             col[count].c   = 0;
    1737:             val_A[count]   = dv * eta_up / (hy * hz);
    1738:             ++count; /* back up */

    1740:             if (!parameters->faces_only) {
    1741:               col[count].i   = ex;
    1742:               col[count].j   = ey;
    1743:               col[count].k   = ez - 1;
    1744:               col[count].loc = DMSTAG_ELEMENT;
    1745:               col[count].c   = 0;
    1746:               val_A[count]   = K_cont * dv / hz;
    1747:               ++count;
    1748:               col[count].i   = ex;
    1749:               col[count].j   = ey;
    1750:               col[count].k   = ez;
    1751:               col[count].loc = DMSTAG_ELEMENT;
    1752:               col[count].c   = 0;
    1753:               val_A[count]   = -1.0 * K_cont * dv / hz;
    1754:               ++count;
    1755:             }

    1757:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
    1758:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1759:           }
    1760:         }

    1762:         /* Elements */
    1763:         if (!parameters->faces_only) {
    1764:           DMStagStencil row;

    1766:           row.i   = ex;
    1767:           row.j   = ey;
    1768:           row.k   = ez;
    1769:           row.loc = DMSTAG_ELEMENT;
    1770:           row.c   = 0;

    1772:           if (ctx->pin_pressure && ex == pinx && ey == piny && ez == pinz) {
    1773:             /* Pin a pressure node to zero, if requested */
    1774:             const PetscScalar val_A   = K_bound;
    1775:             const PetscScalar val_rhs = 0.0;

    1777:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
    1778:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1779:           } else {
    1780:             /* Continuity equation */
    1781:             /* Note that this includes an explicit zero on the diagonal. This is only needed for
    1782:                some direct solvers (not required if using an iterative solver and setting a constant-pressure nullspace) */
    1783:             /* Note: the scaling by dv is not chosen in a principled way and is likely sub-optimal */
    1784:             DMStagStencil     col[7];
    1785:             PetscScalar       val_A[7];
    1786:             const PetscScalar val_rhs = 0.0;

    1788:             col[0].i   = ex;
    1789:             col[0].j   = ey;
    1790:             col[0].k   = ez;
    1791:             col[0].loc = DMSTAG_LEFT;
    1792:             col[0].c   = 0;
    1793:             val_A[0]   = -1.0 * K_cont * dv / hx;
    1794:             col[1].i   = ex;
    1795:             col[1].j   = ey;
    1796:             col[1].k   = ez;
    1797:             col[1].loc = DMSTAG_RIGHT;
    1798:             col[1].c   = 0;
    1799:             val_A[1]   = K_cont * dv / hx;
    1800:             col[2].i   = ex;
    1801:             col[2].j   = ey;
    1802:             col[2].k   = ez;
    1803:             col[2].loc = DMSTAG_DOWN;
    1804:             col[2].c   = 0;
    1805:             val_A[2]   = -1.0 * K_cont * dv / hy;
    1806:             col[3].i   = ex;
    1807:             col[3].j   = ey;
    1808:             col[3].k   = ez;
    1809:             col[3].loc = DMSTAG_UP;
    1810:             col[3].c   = 0;
    1811:             val_A[3]   = K_cont * dv / hy;
    1812:             col[4].i   = ex;
    1813:             col[4].j   = ey;
    1814:             col[4].k   = ez;
    1815:             col[4].loc = DMSTAG_BACK;
    1816:             col[4].c   = 0;
    1817:             val_A[4]   = -1.0 * K_cont * dv / hz;
    1818:             col[5].i   = ex;
    1819:             col[5].j   = ey;
    1820:             col[5].k   = ez;
    1821:             col[5].loc = DMSTAG_FRONT;
    1822:             col[5].c   = 0;
    1823:             val_A[5]   = K_cont * dv / hz;
    1824:             col[6]     = row;
    1825:             val_A[6]   = 0.0;
    1826:             PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 7, col, val_A, INSERT_VALUES));
    1827:             if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
    1828:           }
    1829:         }
    1830:       }
    1831:     }
    1832:   }
    1833:   PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));

    1835:   /* Add additional inverse viscosity terms (for use in building a preconditioning matrix) */
    1836:   if (parameters->include_inverse_visc) {
    1837:     PetscCheck(!parameters->faces_only, PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "Does not make sense with faces only");
    1838:     PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_main, dm_coefficients, ctx->levels[level]->coeff, 1.0, A));
    1839:   }

    1841:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1842:   if (build_rhs) PetscCall(VecAssemblyBegin(rhs));
    1843:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    1844:   if (build_rhs) PetscCall(VecAssemblyEnd(rhs));
    1845:   PetscFunctionReturn(PETSC_SUCCESS);
    1846: }

    1848: static PetscErrorCode CreateSystem(SystemParameters parameters, Mat *pA, Vec *pRhs)
    1849: {
    1850:   PetscFunctionBeginUser;
    1851:   if (parameters->ctx->dim == 2) {
    1852:     PetscCall(CreateSystem2d(parameters, pA, pRhs));
    1853:   } else if (parameters->ctx->dim == 3) {
    1854:     PetscCall(CreateSystem3d(parameters, pA, pRhs));
    1855:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, parameters->ctx->dim);
    1856:   PetscFunctionReturn(PETSC_SUCCESS);
    1857: }

    1859: PetscErrorCode PopulateCoefficientData(Ctx ctx, PetscInt level)
    1860: {
    1861:   PetscInt    dim;
    1862:   PetscInt    N[3];
    1863:   PetscInt    ex, ey, ez, startx, starty, startz, nx, ny, nz;
    1864:   PetscInt    slot_prev, slot_center;
    1865:   PetscInt    slot_rho_downleft, slot_rho_backleft, slot_rho_backdown, slot_eta_element, slot_eta_downleft, slot_eta_backleft, slot_eta_backdown;
    1866:   Vec         coeff_local;
    1867:   PetscReal **arr_coordinates_x, **arr_coordinates_y, **arr_coordinates_z;
    1868:   DM          dm_coefficients;
    1869:   Vec         coeff;

    1871:   PetscFunctionBeginUser;
    1872:   dm_coefficients = ctx->levels[level]->dm_coefficients;
    1873:   PetscCall(DMGetDimension(dm_coefficients, &dim));

    1875:   /* Create global coefficient vector */
    1876:   PetscCall(DMCreateGlobalVector(dm_coefficients, &ctx->levels[level]->coeff));
    1877:   coeff = ctx->levels[level]->coeff;

    1879:   /* Get temporary access to a local representation of the coefficient data */
    1880:   PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
    1881:   PetscCall(DMGlobalToLocal(dm_coefficients, coeff, INSERT_VALUES, coeff_local));

    1883:   /* Use direct array access to coefficient and coordinate arrays, to popoulate coefficient data */
    1884:   PetscCall(DMStagGetGhostCorners(dm_coefficients, &startx, &starty, &startz, &nx, &ny, &nz));
    1885:   PetscCall(DMStagGetGlobalSizes(dm_coefficients, &N[0], &N[1], &N[2]));
    1886:   PetscCall(DMStagGetProductCoordinateArraysRead(dm_coefficients, &arr_coordinates_x, &arr_coordinates_y, &arr_coordinates_z));
    1887:   PetscCall(DMStagGetProductCoordinateLocationSlot(dm_coefficients, DMSTAG_ELEMENT, &slot_center));
    1888:   PetscCall(DMStagGetProductCoordinateLocationSlot(dm_coefficients, DMSTAG_LEFT, &slot_prev));
    1889:   PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_ELEMENT, 0, &slot_eta_element));
    1890:   PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_DOWN_LEFT, 0, &slot_eta_downleft));
    1891:   PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_DOWN_LEFT, 1, &slot_rho_downleft));
    1892:   if (dim == 2) {
    1893:     PetscScalar ***arr_coefficients;

    1895:     PetscCall(DMStagVecGetArray(dm_coefficients, coeff_local, &arr_coefficients));
    1896:     /* Note that these ranges are with respect to the local representation */
    1897:     for (ey = starty; ey < starty + ny; ++ey) {
    1898:       for (ex = startx; ex < startx + nx; ++ex) {
    1899:         arr_coefficients[ey][ex][slot_eta_element]  = ctx->GetEta(ctx, arr_coordinates_x[ex][slot_center], arr_coordinates_y[ey][slot_center], 0.0);
    1900:         arr_coefficients[ey][ex][slot_eta_downleft] = ctx->GetEta(ctx, arr_coordinates_x[ex][slot_prev], arr_coordinates_y[ey][slot_prev], 0.0);
    1901:         arr_coefficients[ey][ex][slot_rho_downleft] = ctx->GetRho(ctx, arr_coordinates_x[ex][slot_prev], arr_coordinates_y[ey][slot_prev], 0.0);
    1902:       }
    1903:     }
    1904:     PetscCall(DMStagVecRestoreArray(dm_coefficients, coeff_local, &arr_coefficients));
    1905:   } else if (dim == 3) {
    1906:     PetscScalar ****arr_coefficients;

    1908:     PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_LEFT, 0, &slot_eta_backleft));
    1909:     PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_LEFT, 1, &slot_rho_backleft));
    1910:     PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_DOWN, 0, &slot_eta_backdown));
    1911:     PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_DOWN, 1, &slot_rho_backdown));
    1912:     PetscCall(DMStagVecGetArray(dm_coefficients, coeff_local, &arr_coefficients));
    1913:     /* Note that these are with respect to the entire local representation, including ghosts */
    1914:     for (ez = startz; ez < startz + nz; ++ez) {
    1915:       for (ey = starty; ey < starty + ny; ++ey) {
    1916:         for (ex = startx; ex < startx + nx; ++ex) {
    1917:           const PetscScalar x_prev   = arr_coordinates_x[ex][slot_prev];
    1918:           const PetscScalar y_prev   = arr_coordinates_y[ey][slot_prev];
    1919:           const PetscScalar z_prev   = arr_coordinates_z[ez][slot_prev];
    1920:           const PetscScalar x_center = arr_coordinates_x[ex][slot_center];
    1921:           const PetscScalar y_center = arr_coordinates_y[ey][slot_center];
    1922:           const PetscScalar z_center = arr_coordinates_z[ez][slot_center];

    1924:           arr_coefficients[ez][ey][ex][slot_eta_element]  = ctx->GetEta(ctx, x_center, y_center, z_center);
    1925:           arr_coefficients[ez][ey][ex][slot_eta_downleft] = ctx->GetEta(ctx, x_prev, y_prev, z_center);
    1926:           arr_coefficients[ez][ey][ex][slot_rho_downleft] = ctx->GetRho(ctx, x_prev, y_prev, z_center);
    1927:           arr_coefficients[ez][ey][ex][slot_eta_backleft] = ctx->GetEta(ctx, x_prev, y_center, z_prev);
    1928:           arr_coefficients[ez][ey][ex][slot_rho_backleft] = ctx->GetRho(ctx, x_prev, y_center, z_prev);
    1929:           arr_coefficients[ez][ey][ex][slot_eta_backdown] = ctx->GetEta(ctx, x_center, y_prev, z_prev);
    1930:           arr_coefficients[ez][ey][ex][slot_rho_backdown] = ctx->GetRho(ctx, x_center, y_prev, z_prev);
    1931:         }
    1932:       }
    1933:     }
    1934:     PetscCall(DMStagVecRestoreArray(dm_coefficients, coeff_local, &arr_coefficients));
    1935:   } else SETERRQ(PetscObjectComm((PetscObject)dm_coefficients), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
    1936:   PetscCall(DMStagRestoreProductCoordinateArraysRead(dm_coefficients, &arr_coordinates_x, &arr_coordinates_y, &arr_coordinates_z));
    1937:   PetscCall(DMLocalToGlobal(dm_coefficients, coeff_local, INSERT_VALUES, coeff));
    1938:   PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));
    1939:   PetscFunctionReturn(PETSC_SUCCESS);
    1940: }

    1942: static PetscErrorCode CreateAuxiliaryOperator(Ctx ctx, PetscInt level, Mat *p_S_hat)
    1943: {
    1944:   DM  dm_element;
    1945:   Mat S_hat;
    1946:   DM  dm_stokes, dm_coefficients;
    1947:   Vec coeff;

    1949:   PetscFunctionBeginUser;
    1950:   dm_stokes       = ctx->levels[level]->dm_stokes;
    1951:   dm_coefficients = ctx->levels[level]->dm_coefficients;
    1952:   coeff           = ctx->levels[level]->coeff;
    1953:   if (ctx->dim == 2) {
    1954:     PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 1, 0, &dm_element));
    1955:   } else if (ctx->dim == 3) {
    1956:     PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 0, 1, &dm_element));
    1957:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
    1958:   PetscCall(DMCreateMatrix(dm_element, p_S_hat));
    1959:   S_hat = *p_S_hat;
    1960:   PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_element, dm_coefficients, coeff, 1.0, S_hat));
    1961:   PetscCall(MatAssemblyBegin(S_hat, MAT_FINAL_ASSEMBLY));
    1962:   PetscCall(MatAssemblyEnd(S_hat, MAT_FINAL_ASSEMBLY));
    1963:   PetscCall(DMDestroy(&dm_element));
    1964:   PetscFunctionReturn(PETSC_SUCCESS);
    1965: }

    1967: static PetscErrorCode OperatorInsertInverseViscosityPressureTerms(DM dm, DM dm_coefficients, Vec coefficients, PetscScalar scale, Mat mat)
    1968: {
    1969:   PetscInt dim, ex, ey, ez, startx, starty, startz, nx, ny, nz;
    1970:   Vec      coeff_local;

    1972:   PetscFunctionBeginUser;
    1973:   PetscCall(DMGetDimension(dm, &dim));
    1974:   PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
    1975:   PetscCall(DMGlobalToLocal(dm_coefficients, coefficients, INSERT_VALUES, coeff_local));
    1976:   PetscCall(DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
    1977:   if (dim == 2) { /* Trick to have one loop nest */
    1978:     startz = 0;
    1979:     nz     = 1;
    1980:   }
    1981:   for (ez = startz; ez < startz + nz; ++ez) {
    1982:     for (ey = starty; ey < starty + ny; ++ey) {
    1983:       for (ex = startx; ex < startx + nx; ++ex) {
    1984:         DMStagStencil from, to;
    1985:         PetscScalar   val;

    1987:         /* component 0 on element is viscosity */
    1988:         from.i   = ex;
    1989:         from.j   = ey;
    1990:         from.k   = ez;
    1991:         from.c   = 0;
    1992:         from.loc = DMSTAG_ELEMENT;
    1993:         PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 1, &from, &val));
    1994:         val = scale / val; /* inverse viscosity, scaled */
    1995:         to  = from;
    1996:         PetscCall(DMStagMatSetValuesStencil(dm, mat, 1, &to, 1, &to, &val, INSERT_VALUES));
    1997:       }
    1998:     }
    1999:   }
    2000:   PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));
    2001:   /* Note that this function does not call MatAssembly{Begin,End} */
    2002:   PetscFunctionReturn(PETSC_SUCCESS);
    2003: }

    2005: /* Create a pressure-only DMStag and use it to generate a nullspace vector
    2006:    - Create a compatible DMStag with one dof per element (and nothing else).
    2007:    - Create a constant vector and normalize it
    2008:    - Migrate it to a vector on the original dmSol (making use of the fact
    2009:    that this will fill in zeros for "extra" dof)
    2010:    - Set the nullspace for the operator
    2011:    - Destroy everything (the operator keeps the references it needs) */
    2012: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
    2013: {
    2014:   DM           dmPressure;
    2015:   Vec          constantPressure, basis;
    2016:   PetscReal    nrm;
    2017:   MatNullSpace matNullSpace;

    2019:   PetscFunctionBeginUser;
    2020:   PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 0, 1, 0, &dmPressure));
    2021:   PetscCall(DMGetGlobalVector(dmPressure, &constantPressure));
    2022:   PetscCall(VecSet(constantPressure, 1.0));
    2023:   PetscCall(VecNorm(constantPressure, NORM_2, &nrm));
    2024:   PetscCall(VecScale(constantPressure, 1.0 / nrm));
    2025:   PetscCall(DMCreateGlobalVector(dmSol, &basis));
    2026:   PetscCall(DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis));
    2027:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace));
    2028:   PetscCall(VecDestroy(&basis));
    2029:   PetscCall(MatSetNullSpace(A, matNullSpace));
    2030:   PetscCall(MatNullSpaceDestroy(&matNullSpace));
    2031:   PetscCall(DMRestoreGlobalVector(dmPressure, &constantPressure));
    2032:   PetscCall(DMDestroy(&dmPressure));
    2033:   PetscFunctionReturn(PETSC_SUCCESS);
    2034: }

    2036: static PetscErrorCode DumpSolution(Ctx ctx, PetscInt level, Vec x)
    2037: {
    2038:   DM       dm_stokes, dm_coefficients;
    2039:   Vec      coeff;
    2040:   DM       dm_vel_avg;
    2041:   Vec      vel_avg;
    2042:   DM       da_vel_avg, da_p, da_eta_element;
    2043:   Vec      vec_vel_avg, vec_p, vec_eta_element;
    2044:   DM       da_eta_down_left, da_rho_down_left, da_eta_back_left, da_rho_back_left, da_eta_back_down, da_rho_back_down;
    2045:   Vec      vec_eta_down_left, vec_rho_down_left, vec_eta_back_left, vec_rho_back_left, vec_eta_back_down, vec_rho_back_down;
    2046:   PetscInt ex, ey, ez, startx, starty, startz, nx, ny, nz;
    2047:   Vec      stokesLocal;

    2049:   PetscFunctionBeginUser;
    2050:   dm_stokes       = ctx->levels[level]->dm_stokes;
    2051:   dm_coefficients = ctx->levels[level]->dm_coefficients;
    2052:   coeff           = ctx->levels[level]->coeff;

    2054:   /* For convenience, create a new DM and Vec which will hold averaged velocities
    2055:      Note that this could also be accomplished with direct array access, using
    2056:      DMStagVecGetArray() and related functions */
    2057:   if (ctx->dim == 2) {
    2058:     PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 2, 0, &dm_vel_avg)); /* 2 dof per element */
    2059:   } else if (ctx->dim == 3) {
    2060:     PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 0, 3, &dm_vel_avg)); /* 3 dof per element */
    2061:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not Implemented for dimension %" PetscInt_FMT, ctx->dim);
    2062:   PetscCall(DMSetUp(dm_vel_avg));
    2063:   PetscCall(DMStagSetUniformCoordinatesProduct(dm_vel_avg, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax, ctx->zmin, ctx->zmax));
    2064:   PetscCall(DMCreateGlobalVector(dm_vel_avg, &vel_avg));
    2065:   PetscCall(DMGetLocalVector(dm_stokes, &stokesLocal));
    2066:   PetscCall(DMGlobalToLocal(dm_stokes, x, INSERT_VALUES, stokesLocal));
    2067:   PetscCall(DMStagGetCorners(dm_vel_avg, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
    2068:   if (ctx->dim == 2) {
    2069:     for (ey = starty; ey < starty + ny; ++ey) {
    2070:       for (ex = startx; ex < startx + nx; ++ex) {
    2071:         DMStagStencil from[4], to[2];
    2072:         PetscScalar   valFrom[4], valTo[2];

    2074:         from[0].i   = ex;
    2075:         from[0].j   = ey;
    2076:         from[0].loc = DMSTAG_UP;
    2077:         from[0].c   = 0;
    2078:         from[1].i   = ex;
    2079:         from[1].j   = ey;
    2080:         from[1].loc = DMSTAG_DOWN;
    2081:         from[1].c   = 0;
    2082:         from[2].i   = ex;
    2083:         from[2].j   = ey;
    2084:         from[2].loc = DMSTAG_LEFT;
    2085:         from[2].c   = 0;
    2086:         from[3].i   = ex;
    2087:         from[3].j   = ey;
    2088:         from[3].loc = DMSTAG_RIGHT;
    2089:         from[3].c   = 0;
    2090:         PetscCall(DMStagVecGetValuesStencil(dm_stokes, stokesLocal, 4, from, valFrom));
    2091:         to[0].i   = ex;
    2092:         to[0].j   = ey;
    2093:         to[0].loc = DMSTAG_ELEMENT;
    2094:         to[0].c   = 0;
    2095:         valTo[0]  = 0.5 * (valFrom[2] + valFrom[3]);
    2096:         to[1].i   = ex;
    2097:         to[1].j   = ey;
    2098:         to[1].loc = DMSTAG_ELEMENT;
    2099:         to[1].c   = 1;
    2100:         valTo[1]  = 0.5 * (valFrom[0] + valFrom[1]);
    2101:         PetscCall(DMStagVecSetValuesStencil(dm_vel_avg, vel_avg, 2, to, valTo, INSERT_VALUES));
    2102:       }
    2103:     }
    2104:   } else if (ctx->dim == 3) {
    2105:     for (ez = startz; ez < startz + nz; ++ez) {
    2106:       for (ey = starty; ey < starty + ny; ++ey) {
    2107:         for (ex = startx; ex < startx + nx; ++ex) {
    2108:           DMStagStencil from[6], to[3];
    2109:           PetscScalar   valFrom[6], valTo[3];

    2111:           from[0].i   = ex;
    2112:           from[0].j   = ey;
    2113:           from[0].k   = ez;
    2114:           from[0].loc = DMSTAG_UP;
    2115:           from[0].c   = 0;
    2116:           from[1].i   = ex;
    2117:           from[1].j   = ey;
    2118:           from[1].k   = ez;
    2119:           from[1].loc = DMSTAG_DOWN;
    2120:           from[1].c   = 0;
    2121:           from[2].i   = ex;
    2122:           from[2].j   = ey;
    2123:           from[2].k   = ez;
    2124:           from[2].loc = DMSTAG_LEFT;
    2125:           from[2].c   = 0;
    2126:           from[3].i   = ex;
    2127:           from[3].j   = ey;
    2128:           from[3].k   = ez;
    2129:           from[3].loc = DMSTAG_RIGHT;
    2130:           from[3].c   = 0;
    2131:           from[4].i   = ex;
    2132:           from[4].j   = ey;
    2133:           from[4].k   = ez;
    2134:           from[4].loc = DMSTAG_BACK;
    2135:           from[4].c   = 0;
    2136:           from[5].i   = ex;
    2137:           from[5].j   = ey;
    2138:           from[5].k   = ez;
    2139:           from[5].loc = DMSTAG_FRONT;
    2140:           from[5].c   = 0;
    2141:           PetscCall(DMStagVecGetValuesStencil(dm_stokes, stokesLocal, 6, from, valFrom));
    2142:           to[0].i   = ex;
    2143:           to[0].j   = ey;
    2144:           to[0].k   = ez;
    2145:           to[0].loc = DMSTAG_ELEMENT;
    2146:           to[0].c   = 0;
    2147:           valTo[0]  = 0.5 * (valFrom[2] + valFrom[3]);
    2148:           to[1].i   = ex;
    2149:           to[1].j   = ey;
    2150:           to[1].k   = ez;
    2151:           to[1].loc = DMSTAG_ELEMENT;
    2152:           to[1].c   = 1;
    2153:           valTo[1]  = 0.5 * (valFrom[0] + valFrom[1]);
    2154:           to[2].i   = ex;
    2155:           to[2].j   = ey;
    2156:           to[2].k   = ez;
    2157:           to[2].loc = DMSTAG_ELEMENT;
    2158:           to[2].c   = 2;
    2159:           valTo[2]  = 0.5 * (valFrom[4] + valFrom[5]);
    2160:           PetscCall(DMStagVecSetValuesStencil(dm_vel_avg, vel_avg, 3, to, valTo, INSERT_VALUES));
    2161:         }
    2162:       }
    2163:     }
    2164:   }
    2165:   PetscCall(VecAssemblyBegin(vel_avg));
    2166:   PetscCall(VecAssemblyEnd(vel_avg));
    2167:   PetscCall(DMRestoreLocalVector(dm_stokes, &stokesLocal));

    2169:   /* Create individual DMDAs for sub-grids of our DMStag objects. This is
    2170:      somewhat inefficient, but allows use of the DMDA API without re-implementing
    2171:      all utilities for DMStag */

    2173:   PetscCall(DMStagVecSplitToDMDA(dm_stokes, x, DMSTAG_ELEMENT, 0, &da_p, &vec_p));
    2174:   PetscCall(PetscObjectSetName((PetscObject)vec_p, "p (scaled)"));

    2176:   PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_ELEMENT, 0, &da_eta_element, &vec_eta_element));
    2177:   PetscCall(PetscObjectSetName((PetscObject)vec_eta_element, "eta"));

    2179:   PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_DOWN_LEFT, 0, &da_eta_down_left, &vec_eta_down_left));
    2180:   PetscCall(PetscObjectSetName((PetscObject)vec_eta_down_left, "eta"));

    2182:   PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_DOWN_LEFT, 1, &da_rho_down_left, &vec_rho_down_left));
    2183:   PetscCall(PetscObjectSetName((PetscObject)vec_rho_down_left, "density"));

    2185:   if (ctx->dim == 3) {
    2186:     PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_LEFT, 0, &da_eta_back_left, &vec_eta_back_left));
    2187:     PetscCall(PetscObjectSetName((PetscObject)vec_eta_back_left, "eta"));

    2189:     PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_LEFT, 1, &da_rho_back_left, &vec_rho_back_left));
    2190:     PetscCall(PetscObjectSetName((PetscObject)vec_rho_back_left, "rho"));

    2192:     PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_DOWN, 0, &da_eta_back_down, &vec_eta_back_down));
    2193:     PetscCall(PetscObjectSetName((PetscObject)vec_eta_back_down, "eta"));

    2195:     PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_DOWN, 1, &da_rho_back_down, &vec_rho_back_down));
    2196:     PetscCall(PetscObjectSetName((PetscObject)vec_rho_back_down, "rho"));
    2197:   }

    2199:   PetscCall(DMStagVecSplitToDMDA(dm_vel_avg, vel_avg, DMSTAG_ELEMENT, -3, &da_vel_avg, &vec_vel_avg)); /* note -3 : pad with zero */
    2200:   PetscCall(PetscObjectSetName((PetscObject)vec_vel_avg, "Velocity (Averaged)"));

    2202:   /* Dump element-based fields to a .vtr file */
    2203:   {
    2204:     PetscViewer viewer;

    2206:     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_vel_avg), "ex4_element.vtr", FILE_MODE_WRITE, &viewer));
    2207:     PetscCall(VecView(vec_vel_avg, viewer));
    2208:     PetscCall(VecView(vec_p, viewer));
    2209:     PetscCall(VecView(vec_eta_element, viewer));
    2210:     PetscCall(PetscViewerDestroy(&viewer));
    2211:   }

    2213:   /* Dump vertex- or edge-based fields to a second .vtr file */
    2214:   {
    2215:     PetscViewer viewer;

    2217:     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_down_left), "ex4_down_left.vtr", FILE_MODE_WRITE, &viewer));
    2218:     PetscCall(VecView(vec_eta_down_left, viewer));
    2219:     PetscCall(VecView(vec_rho_down_left, viewer));
    2220:     PetscCall(PetscViewerDestroy(&viewer));
    2221:   }
    2222:   if (ctx->dim == 3) {
    2223:     PetscViewer viewer;

    2225:     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_back_left), "ex4_back_left.vtr", FILE_MODE_WRITE, &viewer));
    2226:     PetscCall(VecView(vec_eta_back_left, viewer));
    2227:     PetscCall(VecView(vec_rho_back_left, viewer));
    2228:     PetscCall(PetscViewerDestroy(&viewer));
    2229:   }
    2230:   if (ctx->dim == 3) {
    2231:     PetscViewer viewer;

    2233:     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_back_down), "ex4_back_down.vtr", FILE_MODE_WRITE, &viewer));
    2234:     PetscCall(VecView(vec_eta_back_down, viewer));
    2235:     PetscCall(VecView(vec_rho_back_down, viewer));
    2236:     PetscCall(PetscViewerDestroy(&viewer));
    2237:   }

    2239:   /* Destroy DMDAs and Vecs */
    2240:   PetscCall(VecDestroy(&vec_vel_avg));
    2241:   PetscCall(VecDestroy(&vec_p));
    2242:   PetscCall(VecDestroy(&vec_eta_element));
    2243:   PetscCall(VecDestroy(&vec_eta_down_left));
    2244:   if (ctx->dim == 3) {
    2245:     PetscCall(VecDestroy(&vec_eta_back_left));
    2246:     PetscCall(VecDestroy(&vec_eta_back_down));
    2247:   }
    2248:   PetscCall(VecDestroy(&vec_rho_down_left));
    2249:   if (ctx->dim == 3) {
    2250:     PetscCall(VecDestroy(&vec_rho_back_left));
    2251:     PetscCall(VecDestroy(&vec_rho_back_down));
    2252:   }
    2253:   PetscCall(DMDestroy(&da_vel_avg));
    2254:   PetscCall(DMDestroy(&da_p));
    2255:   PetscCall(DMDestroy(&da_eta_element));
    2256:   PetscCall(DMDestroy(&da_eta_down_left));
    2257:   if (ctx->dim == 3) {
    2258:     PetscCall(DMDestroy(&da_eta_back_left));
    2259:     PetscCall(DMDestroy(&da_eta_back_down));
    2260:   }
    2261:   PetscCall(DMDestroy(&da_rho_down_left));
    2262:   if (ctx->dim == 3) {
    2263:     PetscCall(DMDestroy(&da_rho_back_left));
    2264:     PetscCall(DMDestroy(&da_rho_back_down));
    2265:   }
    2266:   PetscCall(VecDestroy(&vel_avg));
    2267:   PetscCall(DMDestroy(&dm_vel_avg));
    2268:   PetscFunctionReturn(PETSC_SUCCESS);
    2269: }

    2271: /*TEST

    2273:    test:
    2274:       suffix: direct_umfpack
    2275:       requires: suitesparse !complex
    2276:       nsize: 1
    2277:       args: -dim 2 -coefficients layers -nondimensional 0 -stag_grid_x 12 -stag_grid_y 7 -pc_type lu -pc_factor_mat_solver_type umfpack -ksp_converged_reason

    2279:    test:
    2280:       suffix: direct_mumps
    2281:       requires: mumps !complex !single
    2282:       nsize: 9
    2283:       args: -dim 2 -coefficients layers -nondimensional 0 -stag_grid_x 13 -stag_grid_y 8 -pc_type lu -pc_factor_mat_solver_type mumps -ksp_converged_reason

    2285:    test:
    2286:       suffix: isovisc_nondim_abf_mg
    2287:       nsize: 1
    2288:       args: -dim 2 -coefficients layers -nondimensional 1 -pc_type fieldsplit -pc_fieldsplit_type schur -ksp_converged_reason -fieldsplit_element_ksp_type preonly -pc_fieldsplit_detect_saddle_point false -fieldsplit_face_pc_type mg -fieldsplit_face_pc_mg_levels 3 -stag_grid_x 24 -stag_grid_y 24 -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_ksp_converged_reason -ksp_type fgmres -fieldsplit_element_pc_type none -fieldsplit_face_mg_levels_ksp_max_it 6 -pc_fieldsplit_schur_fact_type upper -isoviscous

    2290:    test:
    2291:       suffix: isovisc_nondim_abf_mg_2
    2292:       nsize: 1
    2293:       args: -dim 2 -coefficients layers -nondimensional -isoviscous -eta1 1.0 -stag_grid_x 32 -stag_grid_y 32 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -build_auxiliary_operator -fieldsplit_element_ksp_type preonly -fieldsplit_element_pc_type jacobi -fieldsplit_face_pc_type mg -fieldsplit_face_pc_mg_levels 3 -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_mg_levels_pc_type jacobi -fieldsplit_face_mg_levels_ksp_type chebyshev -ksp_converged_reason

    2295:    test:
    2296:       suffix: nondim_abf_mg
    2297:       requires: suitesparse !complex
    2298:       nsize: 4
    2299:       args: -dim 2 -coefficients layers -pc_type fieldsplit -pc_fieldsplit_type schur -ksp_converged_reason -fieldsplit_element_ksp_type preonly -pc_fieldsplit_detect_saddle_point false -fieldsplit_face_pc_type mg -fieldsplit_face_pc_mg_levels 3 -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_mg_coarse_pc_type redundant -fieldsplit_face_mg_coarse_redundant_pc_type lu -fieldsplit_face_mg_coarse_redundant_pc_factor_mat_solver_type umfpack -ksp_type fgmres -fieldsplit_element_pc_type none -fieldsplit_face_mg_levels_ksp_max_it 6 -pc_fieldsplit_schur_fact_type upper -nondimensional -eta1 1e-2 -eta2 1.0 -ksp_monitor -fieldsplit_face_mg_levels_pc_type jacobi -fieldsplit_face_mg_levels_ksp_type gmres -fieldsplit_element_pc_type jacobi -pc_fieldsplit_schur_precondition selfp -stag_grid_x 32 -stag_grid_y 32 -fieldsplit_face_ksp_monitor

    2301:    test:
    2302:       suffix: nondim_abf_lu
    2303:       requires: suitesparse !complex
    2304:       nsize: 1
    2305:       args: -dim 2 -coefficients layers -pc_type fieldsplit -pc_fieldsplit_type schur -ksp_converged_reason -fieldsplit_element_ksp_type preonly -pc_fieldsplit_detect_saddle_point false -ksp_type fgmres -fieldsplit_element_pc_type none -pc_fieldsplit_schur_fact_type upper -nondimensional -eta1 1e-2 -eta2 1.0 -isoviscous 0 -ksp_monitor -fieldsplit_element_pc_type jacobi -build_auxiliary_operator -fieldsplit_face_pc_type lu -fieldsplit_face_pc_factor_mat_solver_type umfpack -stag_grid_x 32 -stag_grid_y 32

    2307:    test:
    2308:       suffix: 3d_nondim_isovisc_abf_mg
    2309:       requires: !single
    2310:       nsize: 1
    2311:       args: -dim 3 -coefficients layers -isoviscous -nondimensional -build_auxiliary_operator -pc_type fieldsplit -pc_fieldsplit_type schur -ksp_converged_reason -fieldsplit_element_ksp_type preonly -pc_fieldsplit_detect_saddle_point false -fieldsplit_face_pc_type mg -fieldsplit_face_pc_mg_levels 3 -s 16 -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_ksp_converged_reason -ksp_type fgmres -fieldsplit_element_pc_type none -fieldsplit_face_mg_levels_ksp_max_it 6 -pc_fieldsplit_schur_fact_type upper

    2313:    test:
    2314:       TODO: unstable across systems
    2315:       suffix: monolithic
    2316:       nsize: 1
    2317:       requires: double !complex
    2318:       args: -dim {{2 3}separate output} -s 16 -custom_pc_mat -pc_type mg -pc_mg_levels 3 -pc_mg_galerkin -mg_levels_ksp_type gmres -mg_levels_ksp_norm_type unpreconditioned -mg_levels_ksp_max_it 10 -mg_levels_pc_type jacobi -ksp_converged_reason

    2320:    test:
    2321:       suffix: 3d_nondim_isovisc_sinker_abf_mg
    2322:       requires: !complex !single
    2323:       nsize: 1
    2324:       args: -dim 3 -coefficients sinker -isoviscous -nondimensional -pc_type fieldsplit -pc_fieldsplit_type schur -ksp_converged_reason -fieldsplit_element_ksp_type preonly -pc_fieldsplit_detect_saddle_point false -fieldsplit_face_pc_type mg -fieldsplit_face_pc_mg_levels 3 -s 16 -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_ksp_converged_reason -ksp_type fgmres -fieldsplit_element_pc_type none -fieldsplit_face_mg_levels_ksp_max_it 6 -pc_fieldsplit_schur_fact_type upper

    2326:    test:
    2327:       TODO: unstable across systems
    2328:       suffix: 3d_nondim_mono_mg_lamemstyle
    2329:       nsize: 1
    2330:       requires: suitesparse
    2331:       args: -dim 3 -coefficients layers -nondimensional -s 16 -custom_pc_mat -pc_type mg -pc_mg_galerkin -pc_mg_levels 2 -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -mg_levels_ksp_richardson_scale 0.5 -mg_levels_ksp_max_it 20 -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason

    2333:    test:
    2334:       suffix: 3d_isovisc_blob_cuda
    2335:       requires: cuda defined(PETSC_USE_LOG)
    2336:       nsize: 2
    2337:       args: -dim 3 -coefficients blob -isoviscous -s 8 -build_auxiliary_operator -fieldsplit_element_ksp_type preonly -fieldsplit_element_pc_type jacobi -fieldsplit_face_ksp_type fgmres -fieldsplit_face_mg_levels_ksp_max_it 6 -fieldsplit_face_mg_levels_ksp_norm_type none -fieldsplit_face_mg_levels_ksp_type chebyshev -fieldsplit_face_mg_levels_pc_type jacobi -fieldsplit_face_pc_mg_galerkin -fieldsplit_face_pc_mg_levels 3 -fieldsplit_face_pc_type mg -pc_fieldsplit_schur_fact_type upper -pc_fieldsplit_schur_precondition user -pc_fieldsplit_type schur -pc_type fieldsplit -dm_mat_type aijcusparse -dm_vec_type cuda -log_view -fieldsplit_face_pc_mg_log
    2338:       filter: awk "/MGInterp/ {print \$NF}"

    2340: TEST*/