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 matrix from which the preconditioner is constructed */
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:       PetscCheck(ctx->dim == 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
515:       ctx->GetRho = GetRho_blob3;
516:       ctx->GetEta = GetEta_blob3;
517:     }
518:     if (is_sinker_box) {
519:       if (ctx->dim == 2) {
520:         ctx->GetRho = GetRho_sinker_box2;
521:         ctx->GetEta = GetEta_sinker_box2;
522:       } else if (ctx->dim == 3) {
523:         ctx->GetRho = GetRho_sinker_box3;
524:         ctx->GetEta = GetEta_sinker_box3;
525:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
526:     }
527:     if (is_sinker_sphere) {
528:       PetscCheck(ctx->dim == 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
529:       ctx->GetRho = GetRho_sinker_sphere3;
530:       ctx->GetEta = GetEta_sinker_sphere3;
531:     }
532:   }

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

542: static PetscErrorCode CtxDestroy(Ctx *p_ctx)
543: {
544:   Ctx ctx;

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

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

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

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

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

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

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

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

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

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

659:         row.i   = ex;
660:         row.j   = ey;
661:         row.loc = DMSTAG_DOWN;
662:         row.c   = 0;

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

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

699:         count = 0;

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

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

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

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

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

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

814:         row.i   = ex;
815:         row.j   = ey;
816:         row.loc = DMSTAG_LEFT;
817:         row.c   = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

1032:   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");
1033:   pinx = 1;
1034:   piny = 0;
1035:   pinz = 0; /* Depends on assertion above that there are at least two element in the x direction */

1037:   /* Loop over all local elements.

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

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

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

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

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

1075:         /* X faces - left*/
1076:         {
1077:           DMStagStencil row;

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

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

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

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

1140:             count = 0;

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

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

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

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

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

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

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

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

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

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

1299:         /* Y faces - down */
1300:         {
1301:           DMStagStencil row;

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

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

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

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

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

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

1392:             count = 0;

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

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

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

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

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

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

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

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

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

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

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

1551:         /* Z faces - back */
1552:         {
1553:           DMStagStencil row;

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

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

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

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

1616:             count = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2034: static PetscErrorCode DumpSolution(Ctx ctx, PetscInt level, Vec x)
2035: {
2036:   DM       dm_stokes, dm_coefficients;
2037:   Vec      coeff;
2038:   DM       dm_vel_avg;
2039:   Vec      vel_avg;
2040:   DM       da_vel_avg, da_p, da_eta_element;
2041:   Vec      vec_vel_avg, vec_p, vec_eta_element;
2042:   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;
2043:   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;
2044:   PetscInt ex, ey, ez, startx, starty, startz, nx, ny, nz;
2045:   Vec      stokesLocal;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2269: /*TEST

2271:    test:
2272:       suffix: direct_umfpack
2273:       requires: suitesparse !complex
2274:       nsize: 1
2275:       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

2277:    test:
2278:       suffix: direct_mumps
2279:       requires: mumps !complex !single
2280:       nsize: 9
2281:       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

2283:    test:
2284:       suffix: isovisc_nondim_abf_mg
2285:       nsize: 1
2286:       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

2288:    test:
2289:       suffix: isovisc_nondim_abf_mg_2
2290:       nsize: 1
2291:       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

2293:    test:
2294:       suffix: nondim_abf_mg
2295:       requires: suitesparse !complex
2296:       nsize: 4
2297:       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

2299:    test:
2300:       suffix: nondim_abf_lu
2301:       requires: suitesparse !complex
2302:       nsize: 1
2303:       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

2305:    test:
2306:       suffix: 3d_nondim_isovisc_abf_mg
2307:       requires: !single
2308:       nsize: 1
2309:       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

2311:    test:
2312:       TODO: unstable across systems
2313:       suffix: monolithic
2314:       nsize: 1
2315:       requires: double !complex
2316:       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

2318:    test:
2319:       suffix: 3d_nondim_isovisc_sinker_abf_mg
2320:       requires: !complex !single
2321:       nsize: 1
2322:       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

2324:    test:
2325:       TODO: unstable across systems
2326:       suffix: 3d_nondim_mono_mg_lamemstyle
2327:       nsize: 1
2328:       requires: suitesparse
2329:       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

2331:    test:
2332:       suffix: 3d_isovisc_blob_cuda
2333:       requires: cuda defined(PETSC_USE_LOG)
2334:       nsize: 2
2335:       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
2336:       filter: awk "/MGInterp/ {print \$NF}"

2338: TEST*/