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";
8
9 #include <petscdm.h>
10 #include <petscksp.h>
11 #include <petscdmstag.h>
12 #include <petscdmda.h>
13
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;
23
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;
38
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;
46
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 *);
58
main(int argc,char ** argv)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;
67
68 PetscFunctionBeginUser;
69 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
70
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));
80
81 /* Populate application context */
82 PetscCall(CtxCreateAndSetFromOptions(&ctx));
83
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));
104
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));
111
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 }
117
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;
123
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 }
136
137 /* Populate coefficient data */
138 for (PetscInt level = 0; level < ctx->n_levels; ++level) PetscCall(PopulateCoefficientData(ctx, level));
139
140 /* Construct main system */
141 {
142 SystemParameters system_parameters;
143
144 PetscCall(SystemParametersCreate(&system_parameters, ctx));
145 PetscCall(CreateSystem(system_parameters, &A, &b));
146 PetscCall(SystemParametersDestroy(&system_parameters));
147 }
148
149 /* Attach a constant-pressure nullspace to the fine-level operator */
150 if (!ctx->pin_pressure) PetscCall(AttachNullspace(dm_stokes, A));
151
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;
158
159 PetscCallMPI(MPI_Comm_size(ctx->comm, &size));
160 if (PetscDefined(HAVE_SUITESPARSE) && size == 1) {
161 PC pc;
162
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;
169
170 PetscCall(KSPGetPC(ksp, &pc));
171 PetscCall(PCSetType(pc, PCLU));
172 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); /* A default, requires MUMPS */
173 }
174 }
175
176 /* Create and set a custom matrix from which the preconditioner is constructed */
177 if (custom_pc_mat) {
178 SystemParameters system_parameters;
179
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 }
188
189 PetscCall(KSPSetDM(ksp, dm_stokes));
190 PetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
191
192 /* Finish setting up solver (can override options set above) */
193 PetscCall(KSPSetFromOptions(ksp));
194
195 /* Additional solver configuration that can involve setting up and CANNOT be
196 overridden from the command line */
197
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;
202
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 }
208
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;
213
214 if (ctx->n_levels < 2) PetscCall(PetscPrintf(ctx->comm, "Warning: not using multiple levels!\n"));
215
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 }
224
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;
235
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 }
243
244 /* Set up to populate enough to define the sub-solver */
245 PetscCall(KSPSetUp(ksp));
246
247 /* Set multigrid components and settings */
248 {
249 KSP *sub_ksp;
250
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;
262
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));
268
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;
274
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 }
284
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 }
293
294 /* Dump solution by converting to DMDAs and dumping */
295 if (dump_solution) PetscCall(DumpSolution(ctx, ctx->n_levels - 1, x));
296
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 }
315
GetEta_constant(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
324
GetRho_layers(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
331
GetEta_layers(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
338
GetRho_sinker_box2(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
347
GetEta_sinker_box2(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
356
GetRho_sinker_box3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
366
GetEta_sinker_box3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
376
GetRho_sinker_sphere3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
386
GetEta_sinker_sphere3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
396
GetEta_blob3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
405
GetRho_blob3(Ctx ctx,PetscScalar x,PetscScalar y,PetscScalar z)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 }
414
LevelCtxCreate(LevelCtx * p_level_ctx)415 static PetscErrorCode LevelCtxCreate(LevelCtx *p_level_ctx)
416 {
417 LevelCtx level_ctx;
418
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 }
427
LevelCtxDestroy(LevelCtx * p_level_ctx)428 static PetscErrorCode LevelCtxDestroy(LevelCtx *p_level_ctx)
429 {
430 LevelCtx level_ctx;
431
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 }
441
CtxCreateAndSetFromOptions(Ctx * p_ctx)442 static PetscErrorCode CtxCreateAndSetFromOptions(Ctx *p_ctx)
443 {
444 Ctx ctx;
445
446 PetscFunctionBeginUser;
447 PetscCall(PetscMalloc1(1, p_ctx));
448 ctx = *p_ctx;
449
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;
465
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;
487
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;
501
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));
508
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 }
533
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 }
541
CtxDestroy(Ctx * p_ctx)542 static PetscErrorCode CtxDestroy(Ctx *p_ctx)
543 {
544 Ctx ctx;
545
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 }
553
SystemParametersCreate(SystemParameters * parameters,Ctx ctx)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 }
564
SystemParametersDestroy(SystemParameters * parameters)565 static PetscErrorCode SystemParametersDestroy(SystemParameters *parameters)
566 {
567 PetscFunctionBeginUser;
568 PetscCall(PetscFree(*parameters));
569 *parameters = NULL;
570 PetscFunctionReturn(PETSC_SUCCESS);
571 }
572
CreateSystem2d(SystemParameters parameters,Mat * pA,Vec * pRhs)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;
586
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));
613
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);
621
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;
627
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 }
635
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;
641
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;
658
659 row.i = ex;
660 row.j = ey;
661 row.loc = DMSTAG_DOWN;
662 row.c = 0;
663
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]);
675
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];
698
699 count = 0;
700
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;
706
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 }
773
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 }
777
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;
784
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;
797
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;
813
814 row.i = ex;
815 row.j = ey;
816 row.loc = DMSTAG_LEFT;
817 row.c = 0;
818
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];
841
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;
848
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 }
915
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 }
919
920 /* P equation : u_x + v_y = 0
921
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)
924
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;
932
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;
943
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));
977
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 }
983
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 }
990
CreateSystem3d(SystemParameters parameters,Mat * pA,Vec * pRhs)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;
1004
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));
1031
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 */
1036
1037 /* Loop over all local elements.
1038
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)
1043
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);
1054
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 */
1058
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;
1065
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 }
1074
1075 /* X faces - left*/
1076 {
1077 DMStagStencil row;
1078
1079 row.i = ex;
1080 row.j = ey;
1081 row.k = ez;
1082 row.loc = DMSTAG_LEFT;
1083 row.c = 0;
1084
1085 if (left_boundary) {
1086 /* Left x-velocity Dirichlet */
1087 const PetscScalar val_rhs = 0.0;
1088 const PetscScalar val_A = K_bound;
1089
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;
1100
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];
1139
1140 count = 0;
1141
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;
1149
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 }
1200
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 */
1215
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 */
1230
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 */
1245
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 */
1260
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 }
1277
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 }
1282
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;
1289
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 }
1298
1299 /* Y faces - down */
1300 {
1301 DMStagStencil row;
1302
1303 row.i = ex;
1304 row.j = ey;
1305 row.k = ez;
1306 row.loc = DMSTAG_DOWN;
1307 row.c = 0;
1308
1309 if (bottom_boundary) {
1310 /* Bottom y-velocity Dirichlet */
1311 const PetscScalar val_rhs = 0.0;
1312 const PetscScalar val_A = K_bound;
1313
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 */
1323
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));
1348
1349 /* Compute forcing */
1350 val_rhs = ctx->gy * dv * (0.25 * (rho[0] + rho[1] + rho[2] + rho[3]));
1351 }
1352
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];
1391
1392 count = 0;
1393
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;
1401
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;
1416
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 }
1453
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*/
1468
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*/
1483
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 */
1498
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 */
1513
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 }
1530
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 }
1535
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;
1541
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 }
1550
1551 /* Z faces - back */
1552 {
1553 DMStagStencil row;
1554
1555 row.i = ex;
1556 row.j = ey;
1557 row.k = ez;
1558 row.loc = DMSTAG_BACK;
1559 row.c = 0;
1560
1561 if (back_boundary) {
1562 /* Back z-velocity Dirichlet */
1563 const PetscScalar val_rhs = 0.0;
1564 const PetscScalar val_A = K_bound;
1565
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;
1576
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];
1615
1616 count = 0;
1617
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;
1625
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;
1640
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 }
1677
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*/
1692
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*/
1707
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 */
1722
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 */
1737
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 }
1754
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 }
1759
1760 /* Elements */
1761 if (!parameters->faces_only) {
1762 DMStagStencil row;
1763
1764 row.i = ex;
1765 row.j = ey;
1766 row.k = ez;
1767 row.loc = DMSTAG_ELEMENT;
1768 row.c = 0;
1769
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;
1774
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;
1785
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));
1832
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 }
1838
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 }
1845
CreateSystem(SystemParameters parameters,Mat * pA,Vec * pRhs)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 }
1856
PopulateCoefficientData(Ctx ctx,PetscInt level)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;
1868
1869 PetscFunctionBeginUser;
1870 dm_coefficients = ctx->levels[level]->dm_coefficients;
1871 PetscCall(DMGetDimension(dm_coefficients, &dim));
1872
1873 /* Create global coefficient vector */
1874 PetscCall(DMCreateGlobalVector(dm_coefficients, &ctx->levels[level]->coeff));
1875 coeff = ctx->levels[level]->coeff;
1876
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));
1880
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;
1892
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;
1905
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];
1921
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 }
1939
CreateAuxiliaryOperator(Ctx ctx,PetscInt level,Mat * p_S_hat)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;
1946
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 }
1964
OperatorInsertInverseViscosityPressureTerms(DM dm,DM dm_coefficients,Vec coefficients,PetscScalar scale,Mat mat)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;
1969
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;
1984
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 }
2002
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) */
AttachNullspace(DM dmSol,Mat A)2010 static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
2011 {
2012 DM dmPressure;
2013 Vec constantPressure, basis;
2014 PetscReal nrm;
2015 MatNullSpace matNullSpace;
2016
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 }
2033
DumpSolution(Ctx ctx,PetscInt level,Vec x)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;
2046
2047 PetscFunctionBeginUser;
2048 dm_stokes = ctx->levels[level]->dm_stokes;
2049 dm_coefficients = ctx->levels[level]->dm_coefficients;
2050 coeff = ctx->levels[level]->coeff;
2051
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];
2071
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];
2108
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));
2166
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 */
2170
2171 PetscCall(DMStagVecSplitToDMDA(dm_stokes, x, DMSTAG_ELEMENT, 0, &da_p, &vec_p));
2172 PetscCall(PetscObjectSetName((PetscObject)vec_p, "p (scaled)"));
2173
2174 PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_ELEMENT, 0, &da_eta_element, &vec_eta_element));
2175 PetscCall(PetscObjectSetName((PetscObject)vec_eta_element, "eta"));
2176
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"));
2179
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"));
2182
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"));
2186
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"));
2189
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"));
2192
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 }
2196
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)"));
2199
2200 /* Dump element-based fields to a .vtr file */
2201 {
2202 PetscViewer viewer;
2203
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 }
2210
2211 /* Dump vertex- or edge-based fields to a second .vtr file */
2212 {
2213 PetscViewer viewer;
2214
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;
2222
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;
2230
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 }
2236
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 }
2268
2269 /*TEST
2270
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
2276
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
2282
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
2287
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
2292
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
2298
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
2304
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
2310
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
2317
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
2323
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
2330
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}"
2337
2338 TEST*/
2339