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, KSP_DMACTIVE_ALL, 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) PetscCall(MatDestroy(&A_faces[level]));
302: PetscCall(PetscFree(A_faces));
303: }
304: PetscCall(MatDestroy(&P));
305: PetscCall(VecDestroy(&x));
306: PetscCall(VecDestroy(&b));
307: PetscCall(MatDestroy(&S_hat));
308: PetscCall(KSPDestroy(&ksp));
309: PetscCall(CtxDestroy(&ctx));
310: PetscCall(PetscFinalize());
311: return 0;
312: }
314: static PetscScalar GetEta_constant(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
315: {
316: (void)ctx;
317: (void)x;
318: (void)y;
319: (void)z;
320: return 1.0;
321: }
323: static PetscScalar GetRho_layers(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
324: {
325: (void)y;
326: (void)z;
327: return PetscRealPart(x) < (ctx->xmax - ctx->xmin) / 2.0 ? ctx->rho1 : ctx->rho2;
328: }
330: static PetscScalar GetEta_layers(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
331: {
332: (void)y;
333: (void)z;
334: return PetscRealPart(x) < (ctx->xmax - ctx->xmin) / 2.0 ? ctx->eta1 : ctx->eta2;
335: }
337: static PetscScalar GetRho_sinker_box2(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
338: {
339: (void)z;
340: const PetscReal d = ctx->xmax - ctx->xmin;
341: const PetscReal xx = PetscRealPart(x) / d - 0.5;
342: const PetscReal yy = PetscRealPart(y) / d - 0.5;
343: return (xx * xx > 0.15 * 0.15 || yy * yy > 0.15 * 0.15) ? ctx->rho1 : ctx->rho2;
344: }
346: static PetscScalar GetEta_sinker_box2(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
347: {
348: (void)z;
349: const PetscReal d = ctx->xmax - ctx->xmin;
350: const PetscReal xx = PetscRealPart(x) / d - 0.5;
351: const PetscReal yy = PetscRealPart(y) / d - 0.5;
352: return (xx * xx > 0.15 * 0.15 || yy * yy > 0.15 * 0.15) ? ctx->eta1 : ctx->eta2;
353: }
355: static PetscScalar GetRho_sinker_box3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
356: {
357: const PetscReal d = ctx->xmax - ctx->xmin;
358: const PetscReal xx = PetscRealPart(x) / d - 0.5;
359: const PetscReal yy = PetscRealPart(y) / d - 0.5;
360: const PetscReal zz = PetscRealPart(z) / d - 0.5;
361: const PetscReal half_width = 0.15;
362: return (PetscAbsReal(xx) > half_width || PetscAbsReal(yy) > half_width || PetscAbsReal(zz) > half_width) ? ctx->rho1 : ctx->rho2;
363: }
365: static PetscScalar GetEta_sinker_box3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
366: {
367: const PetscReal d = ctx->xmax - ctx->xmin;
368: const PetscReal xx = PetscRealPart(x) / d - 0.5;
369: const PetscReal yy = PetscRealPart(y) / d - 0.5;
370: const PetscReal zz = PetscRealPart(z) / d - 0.5;
371: const PetscReal half_width = 0.15;
372: return (PetscAbsReal(xx) > half_width || PetscAbsReal(yy) > half_width || PetscAbsReal(zz) > half_width) ? ctx->eta1 : ctx->eta2;
373: }
375: static PetscScalar GetRho_sinker_sphere3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
376: {
377: const PetscReal d = ctx->xmax - ctx->xmin;
378: const PetscReal xx = PetscRealPart(x) / d - 0.5;
379: const PetscReal yy = PetscRealPart(y) / d - 0.5;
380: const PetscReal zz = PetscRealPart(z) / d - 0.5;
381: const PetscReal half_width = 0.3;
382: return (xx * xx + yy * yy + zz * zz > half_width * half_width) ? ctx->rho1 : ctx->rho2;
383: }
385: static PetscScalar GetEta_sinker_sphere3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
386: {
387: const PetscReal d = ctx->xmax - ctx->xmin;
388: const PetscReal xx = PetscRealPart(x) / d - 0.5;
389: const PetscReal yy = PetscRealPart(y) / d - 0.5;
390: const PetscReal zz = PetscRealPart(z) / d - 0.5;
391: const PetscReal half_width = 0.3;
392: return (xx * xx + yy * yy + zz * zz > half_width * half_width) ? ctx->eta1 : ctx->eta2;
393: }
395: static PetscScalar GetEta_blob3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
396: {
397: const PetscReal d = ctx->xmax - ctx->xmin;
398: const PetscReal xx = PetscRealPart(x) / d - 0.5;
399: const PetscReal yy = PetscRealPart(y) / d - 0.5;
400: const PetscReal zz = PetscRealPart(z) / d - 0.5;
401: return ctx->eta1 + ctx->eta2 * PetscExpScalar(-20.0 * (xx * xx + yy * yy + zz * zz));
402: }
404: static PetscScalar GetRho_blob3(Ctx ctx, PetscScalar x, PetscScalar y, PetscScalar z)
405: {
406: const PetscReal d = ctx->xmax - ctx->xmin;
407: const PetscReal xx = PetscRealPart(x) / d - 0.5;
408: const PetscReal yy = PetscRealPart(y) / d - 0.5;
409: const PetscReal zz = PetscRealPart(z) / d - 0.5;
410: return ctx->rho1 + ctx->rho2 * PetscExpScalar(-20.0 * (xx * xx + yy * yy + zz * zz));
411: }
413: static PetscErrorCode LevelCtxCreate(LevelCtx *p_level_ctx)
414: {
415: LevelCtx level_ctx;
417: PetscFunctionBeginUser;
418: PetscCall(PetscMalloc1(1, p_level_ctx));
419: level_ctx = *p_level_ctx;
420: level_ctx->dm_stokes = NULL;
421: level_ctx->dm_coefficients = NULL;
422: level_ctx->dm_faces = NULL;
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode LevelCtxDestroy(LevelCtx *p_level_ctx)
427: {
428: LevelCtx level_ctx;
430: PetscFunctionBeginUser;
431: level_ctx = *p_level_ctx;
432: PetscCall(DMDestroy(&level_ctx->dm_stokes));
433: PetscCall(DMDestroy(&level_ctx->dm_coefficients));
434: PetscCall(DMDestroy(&level_ctx->dm_faces));
435: PetscCall(VecDestroy(&level_ctx->coeff));
436: PetscCall(PetscFree(*p_level_ctx));
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: static PetscErrorCode CtxCreateAndSetFromOptions(Ctx *p_ctx)
441: {
442: Ctx ctx;
444: PetscFunctionBeginUser;
445: PetscCall(PetscMalloc1(1, p_ctx));
446: ctx = *p_ctx;
448: ctx->comm = PETSC_COMM_WORLD;
449: ctx->pin_pressure = PETSC_FALSE;
450: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pin_pressure", &ctx->pin_pressure, NULL));
451: ctx->dim = 3;
452: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &ctx->dim, NULL));
453: if (ctx->dim <= 2) {
454: ctx->cells_x = 32;
455: } else {
456: ctx->cells_x = 16;
457: }
458: PetscCall(PetscOptionsGetInt(NULL, NULL, "-s", &ctx->cells_x, NULL)); /* shortcut. Usually, use -stag_grid_x etc. */
459: ctx->cells_z = ctx->cells_y = ctx->cells_x;
460: ctx->xmin = ctx->ymin = ctx->zmin = 0.0;
461: {
462: PetscBool nondimensional = PETSC_TRUE;
464: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nondimensional", &nondimensional, NULL));
465: if (nondimensional) {
466: ctx->xmax = ctx->ymax = ctx->zmax = 1.0;
467: ctx->rho1 = 0.0;
468: ctx->rho2 = 1.0;
469: ctx->eta1 = 1.0;
470: ctx->eta2 = 1e2;
471: ctx->gy = -1.0; /* downwards */
472: } else {
473: ctx->xmax = 1e6;
474: ctx->ymax = 1.5e6;
475: ctx->zmax = 1e6;
476: ctx->rho1 = 3200;
477: ctx->rho2 = 3300;
478: ctx->eta1 = 1e20;
479: ctx->eta2 = 1e22;
480: ctx->gy = -10.0; /* downwards */
481: }
482: }
483: {
484: PetscBool isoviscous;
486: isoviscous = PETSC_FALSE;
487: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-eta1", &ctx->eta1, NULL));
488: PetscCall(PetscOptionsGetBool(NULL, NULL, "-isoviscous", &isoviscous, NULL));
489: if (isoviscous) {
490: ctx->eta2 = ctx->eta1;
491: ctx->GetEta = GetEta_constant; /* override */
492: } else {
493: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-eta2", &ctx->eta2, NULL));
494: }
495: }
496: {
497: char mode[1024] = "sinker";
498: PetscBool is_layers, is_blob, is_sinker_box, is_sinker_sphere;
500: PetscCall(PetscOptionsGetString(NULL, NULL, "-coefficients", mode, sizeof(mode), NULL));
501: PetscCall(PetscStrncmp(mode, "layers", sizeof(mode), &is_layers));
502: PetscCall(PetscStrncmp(mode, "sinker", sizeof(mode), &is_sinker_box));
503: if (!is_sinker_box) PetscCall(PetscStrncmp(mode, "sinker_box", sizeof(mode), &is_sinker_box));
504: PetscCall(PetscStrncmp(mode, "sinker_sphere", sizeof(mode), &is_sinker_sphere));
505: PetscCall(PetscStrncmp(mode, "blob", sizeof(mode), &is_blob));
507: if (is_layers) {
508: ctx->GetRho = GetRho_layers;
509: ctx->GetEta = GetEta_layers;
510: }
511: if (is_blob) {
512: PetscCheck(ctx->dim == 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
513: ctx->GetRho = GetRho_blob3;
514: ctx->GetEta = GetEta_blob3;
515: }
516: if (is_sinker_box) {
517: if (ctx->dim == 2) {
518: ctx->GetRho = GetRho_sinker_box2;
519: ctx->GetEta = GetEta_sinker_box2;
520: } else if (ctx->dim == 3) {
521: ctx->GetRho = GetRho_sinker_box3;
522: ctx->GetEta = GetEta_sinker_box3;
523: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
524: }
525: if (is_sinker_sphere) {
526: PetscCheck(ctx->dim == 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
527: ctx->GetRho = GetRho_sinker_sphere3;
528: ctx->GetEta = GetEta_sinker_sphere3;
529: }
530: }
532: /* Per-level data */
533: ctx->n_levels = 1;
534: PetscCall(PetscOptionsGetInt(NULL, NULL, "-levels", &ctx->n_levels, NULL));
535: PetscCall(PetscMalloc1(ctx->n_levels, &ctx->levels));
536: for (PetscInt i = 0; i < ctx->n_levels; ++i) PetscCall(LevelCtxCreate(&ctx->levels[i]));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode CtxDestroy(Ctx *p_ctx)
541: {
542: Ctx ctx;
544: PetscFunctionBeginUser;
545: ctx = *p_ctx;
546: for (PetscInt i = 0; i < ctx->n_levels; ++i) PetscCall(LevelCtxDestroy(&ctx->levels[i]));
547: PetscCall(PetscFree(ctx->levels));
548: PetscCall(PetscFree(*p_ctx));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode SystemParametersCreate(SystemParameters *parameters, Ctx ctx)
553: {
554: PetscFunctionBeginUser;
555: PetscCall(PetscMalloc1(1, parameters));
556: (*parameters)->ctx = ctx;
557: (*parameters)->level = ctx->n_levels - 1;
558: (*parameters)->include_inverse_visc = PETSC_FALSE;
559: (*parameters)->faces_only = PETSC_FALSE;
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: static PetscErrorCode SystemParametersDestroy(SystemParameters *parameters)
564: {
565: PetscFunctionBeginUser;
566: PetscCall(PetscFree(*parameters));
567: *parameters = NULL;
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode CreateSystem2d(SystemParameters parameters, Mat *pA, Vec *pRhs)
572: {
573: PetscInt N[2];
574: PetscInt ex, ey, startx, starty, nx, ny;
575: Mat A;
576: Vec rhs;
577: PetscReal hx, hy, dv;
578: Vec coefficients_local;
579: PetscBool build_rhs;
580: DM dm_main, dm_coefficients;
581: PetscScalar K_cont, K_bound;
582: Ctx ctx = parameters->ctx;
583: PetscInt level = parameters->level;
585: PetscFunctionBeginUser;
586: if (parameters->faces_only) {
587: dm_main = ctx->levels[level]->dm_faces;
588: } else {
589: dm_main = ctx->levels[level]->dm_stokes;
590: }
591: dm_coefficients = ctx->levels[level]->dm_coefficients;
592: K_cont = ctx->levels[level]->K_cont;
593: K_bound = ctx->levels[level]->K_bound;
594: PetscCall(DMCreateMatrix(dm_main, pA));
595: A = *pA;
596: build_rhs = (PetscBool)(pRhs != NULL);
597: PetscCheck(!(parameters->faces_only && build_rhs), PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "RHS for faces-only not supported");
598: if (build_rhs) {
599: PetscCall(DMCreateGlobalVector(dm_main, pRhs));
600: rhs = *pRhs;
601: } else {
602: rhs = NULL;
603: }
604: PetscCall(DMStagGetCorners(dm_main, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
605: PetscCall(DMStagGetGlobalSizes(dm_main, &N[0], &N[1], NULL));
606: hx = ctx->levels[level]->hx_characteristic;
607: hy = ctx->levels[level]->hy_characteristic;
608: dv = hx * hy;
609: PetscCall(DMGetLocalVector(dm_coefficients, &coefficients_local));
610: PetscCall(DMGlobalToLocal(dm_coefficients, ctx->levels[level]->coeff, INSERT_VALUES, coefficients_local));
612: /* Loop over all local elements */
613: for (ey = starty; ey < starty + ny; ++ey) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
614: for (ex = startx; ex < startx + nx; ++ex) {
615: const PetscBool left_boundary = (PetscBool)(ex == 0);
616: const PetscBool right_boundary = (PetscBool)(ex == N[0] - 1);
617: const PetscBool bottom_boundary = (PetscBool)(ey == 0);
618: const PetscBool top_boundary = (PetscBool)(ey == N[1] - 1);
620: if (ey == N[1] - 1) {
621: /* Top boundary velocity Dirichlet */
622: DMStagStencil row;
623: const PetscScalar val_rhs = 0.0;
624: const PetscScalar val_A = K_bound;
626: row.i = ex;
627: row.j = ey;
628: row.loc = DMSTAG_UP;
629: row.c = 0;
630: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
631: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
632: }
634: if (ey == 0) {
635: /* Bottom boundary velocity Dirichlet */
636: DMStagStencil row;
637: const PetscScalar val_rhs = 0.0;
638: const PetscScalar val_A = K_bound;
640: row.i = ex;
641: row.j = ey;
642: row.loc = DMSTAG_DOWN;
643: row.c = 0;
644: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
645: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
646: } else {
647: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y
648: includes non-zero forcing and free-slip boundary conditions */
649: PetscInt count;
650: DMStagStencil row, col[11];
651: PetscScalar val_A[11];
652: DMStagStencil rhoPoint[2];
653: PetscScalar rho[2], val_rhs;
654: DMStagStencil etaPoint[4];
655: PetscScalar eta[4], eta_left, eta_right, eta_up, eta_down;
657: row.i = ex;
658: row.j = ey;
659: row.loc = DMSTAG_DOWN;
660: row.c = 0;
662: /* get rho values and compute rhs value*/
663: rhoPoint[0].i = ex;
664: rhoPoint[0].j = ey;
665: rhoPoint[0].loc = DMSTAG_DOWN_LEFT;
666: rhoPoint[0].c = 1;
667: rhoPoint[1].i = ex;
668: rhoPoint[1].j = ey;
669: rhoPoint[1].loc = DMSTAG_DOWN_RIGHT;
670: rhoPoint[1].c = 1;
671: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 2, rhoPoint, rho));
672: val_rhs = -ctx->gy * dv * 0.5 * (rho[0] + rho[1]);
674: /* Get eta values */
675: etaPoint[0].i = ex;
676: etaPoint[0].j = ey;
677: etaPoint[0].loc = DMSTAG_DOWN_LEFT;
678: etaPoint[0].c = 0; /* Left */
679: etaPoint[1].i = ex;
680: etaPoint[1].j = ey;
681: etaPoint[1].loc = DMSTAG_DOWN_RIGHT;
682: etaPoint[1].c = 0; /* Right */
683: etaPoint[2].i = ex;
684: etaPoint[2].j = ey;
685: etaPoint[2].loc = DMSTAG_ELEMENT;
686: etaPoint[2].c = 0; /* Up */
687: etaPoint[3].i = ex;
688: etaPoint[3].j = ey - 1;
689: etaPoint[3].loc = DMSTAG_ELEMENT;
690: etaPoint[3].c = 0; /* Down */
691: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 4, etaPoint, eta));
692: eta_left = eta[0];
693: eta_right = eta[1];
694: eta_up = eta[2];
695: eta_down = eta[3];
697: count = 0;
699: col[count] = row;
700: val_A[count] = -2.0 * dv * (eta_down + eta_up) / (hy * hy);
701: if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
702: if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
703: ++count;
705: col[count].i = ex;
706: col[count].j = ey - 1;
707: col[count].loc = DMSTAG_DOWN;
708: col[count].c = 0;
709: val_A[count] = 2.0 * dv * eta_down / (hy * hy);
710: ++count;
711: col[count].i = ex;
712: col[count].j = ey + 1;
713: col[count].loc = DMSTAG_DOWN;
714: col[count].c = 0;
715: val_A[count] = 2.0 * dv * eta_up / (hy * hy);
716: ++count;
717: if (!left_boundary) {
718: col[count].i = ex - 1;
719: col[count].j = ey;
720: col[count].loc = DMSTAG_DOWN;
721: col[count].c = 0;
722: val_A[count] = dv * eta_left / (hx * hx);
723: ++count;
724: }
725: if (!right_boundary) {
726: col[count].i = ex + 1;
727: col[count].j = ey;
728: col[count].loc = DMSTAG_DOWN;
729: col[count].c = 0;
730: val_A[count] = dv * eta_right / (hx * hx);
731: ++count;
732: }
733: col[count].i = ex;
734: col[count].j = ey - 1;
735: col[count].loc = DMSTAG_LEFT;
736: col[count].c = 0;
737: val_A[count] = dv * eta_left / (hx * hy);
738: ++count; /* down left x edge */
739: col[count].i = ex;
740: col[count].j = ey - 1;
741: col[count].loc = DMSTAG_RIGHT;
742: col[count].c = 0;
743: val_A[count] = -1.0 * dv * eta_right / (hx * hy);
744: ++count; /* down right x edge */
745: col[count].i = ex;
746: col[count].j = ey;
747: col[count].loc = DMSTAG_LEFT;
748: col[count].c = 0;
749: val_A[count] = -1.0 * dv * eta_left / (hx * hy);
750: ++count; /* up left x edge */
751: col[count].i = ex;
752: col[count].j = ey;
753: col[count].loc = DMSTAG_RIGHT;
754: col[count].c = 0;
755: val_A[count] = dv * eta_right / (hx * hy);
756: ++count; /* up right x edge */
757: if (!parameters->faces_only) {
758: col[count].i = ex;
759: col[count].j = ey - 1;
760: col[count].loc = DMSTAG_ELEMENT;
761: col[count].c = 0;
762: val_A[count] = K_cont * dv / hy;
763: ++count;
764: col[count].i = ex;
765: col[count].j = ey;
766: col[count].loc = DMSTAG_ELEMENT;
767: col[count].c = 0;
768: val_A[count] = -1.0 * K_cont * dv / hy;
769: ++count;
770: }
772: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
773: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
774: }
776: if (ex == N[0] - 1) {
777: /* Right Boundary velocity Dirichlet */
778: /* Redundant in the corner */
779: DMStagStencil row;
780: const PetscScalar val_rhs = 0.0;
781: const PetscScalar val_A = K_bound;
783: row.i = ex;
784: row.j = ey;
785: row.loc = DMSTAG_RIGHT;
786: row.c = 0;
787: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
788: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
789: }
790: if (ex == 0) {
791: /* Left velocity Dirichlet */
792: DMStagStencil row;
793: const PetscScalar val_rhs = 0.0;
794: const PetscScalar val_A = K_bound;
796: row.i = ex;
797: row.j = ey;
798: row.loc = DMSTAG_LEFT;
799: row.c = 0;
800: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
801: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
802: } else {
803: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x
804: zero RHS, including free-slip boundary conditions */
805: PetscInt count;
806: DMStagStencil row, col[11];
807: PetscScalar val_A[11];
808: DMStagStencil etaPoint[4];
809: PetscScalar eta[4], eta_left, eta_right, eta_up, eta_down;
810: const PetscScalar val_rhs = 0.0;
812: row.i = ex;
813: row.j = ey;
814: row.loc = DMSTAG_LEFT;
815: row.c = 0;
817: /* Get eta values */
818: etaPoint[0].i = ex - 1;
819: etaPoint[0].j = ey;
820: etaPoint[0].loc = DMSTAG_ELEMENT;
821: etaPoint[0].c = 0; /* Left */
822: etaPoint[1].i = ex;
823: etaPoint[1].j = ey;
824: etaPoint[1].loc = DMSTAG_ELEMENT;
825: etaPoint[1].c = 0; /* Right */
826: etaPoint[2].i = ex;
827: etaPoint[2].j = ey;
828: etaPoint[2].loc = DMSTAG_UP_LEFT;
829: etaPoint[2].c = 0; /* Up */
830: etaPoint[3].i = ex;
831: etaPoint[3].j = ey;
832: etaPoint[3].loc = DMSTAG_DOWN_LEFT;
833: etaPoint[3].c = 0; /* Down */
834: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coefficients_local, 4, etaPoint, eta));
835: eta_left = eta[0];
836: eta_right = eta[1];
837: eta_up = eta[2];
838: eta_down = eta[3];
840: count = 0;
841: col[count] = row;
842: val_A[count] = -2.0 * dv * (eta_left + eta_right) / (hx * hx);
843: if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
844: if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
845: ++count;
847: if (!bottom_boundary) {
848: col[count].i = ex;
849: col[count].j = ey - 1;
850: col[count].loc = DMSTAG_LEFT;
851: col[count].c = 0;
852: val_A[count] = dv * eta_down / (hy * hy);
853: ++count;
854: }
855: if (!top_boundary) {
856: col[count].i = ex;
857: col[count].j = ey + 1;
858: col[count].loc = DMSTAG_LEFT;
859: col[count].c = 0;
860: val_A[count] = dv * eta_up / (hy * hy);
861: ++count;
862: }
863: col[count].i = ex - 1;
864: col[count].j = ey;
865: col[count].loc = DMSTAG_LEFT;
866: col[count].c = 0;
867: val_A[count] = 2.0 * dv * eta_left / (hx * hx);
868: ++count;
869: col[count].i = ex + 1;
870: col[count].j = ey;
871: col[count].loc = DMSTAG_LEFT;
872: col[count].c = 0;
873: val_A[count] = 2.0 * dv * eta_right / (hx * hx);
874: ++count;
875: col[count].i = ex - 1;
876: col[count].j = ey;
877: col[count].loc = DMSTAG_DOWN;
878: col[count].c = 0;
879: val_A[count] = dv * eta_down / (hx * hy);
880: ++count; /* down left */
881: col[count].i = ex;
882: col[count].j = ey;
883: col[count].loc = DMSTAG_DOWN;
884: col[count].c = 0;
885: val_A[count] = -1.0 * dv * eta_down / (hx * hy);
886: ++count; /* down right */
887: col[count].i = ex - 1;
888: col[count].j = ey;
889: col[count].loc = DMSTAG_UP;
890: col[count].c = 0;
891: val_A[count] = -1.0 * dv * eta_up / (hx * hy);
892: ++count; /* up left */
893: col[count].i = ex;
894: col[count].j = ey;
895: col[count].loc = DMSTAG_UP;
896: col[count].c = 0;
897: val_A[count] = dv * eta_up / (hx * hy);
898: ++count; /* up right */
899: if (!parameters->faces_only) {
900: col[count].i = ex - 1;
901: col[count].j = ey;
902: col[count].loc = DMSTAG_ELEMENT;
903: col[count].c = 0;
904: val_A[count] = K_cont * dv / hx;
905: ++count;
906: col[count].i = ex;
907: col[count].j = ey;
908: col[count].loc = DMSTAG_ELEMENT;
909: col[count].c = 0;
910: val_A[count] = -1.0 * K_cont * dv / hx;
911: ++count;
912: }
914: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
915: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
916: }
918: /* P equation : u_x + v_y = 0
920: Note that this includes an explicit zero on the diagonal. This is only needed for
921: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace)
923: Note: the scaling by dv is not chosen in a principled way and is likely sub-optimal
924: */
925: if (!parameters->faces_only) {
926: if (ctx->pin_pressure && ex == 0 && ey == 0) { /* Pin the first pressure node to zero, if requested */
927: DMStagStencil row;
928: const PetscScalar val_A = K_bound;
929: const PetscScalar val_rhs = 0.0;
931: row.i = ex;
932: row.j = ey;
933: row.loc = DMSTAG_ELEMENT;
934: row.c = 0;
935: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
936: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
937: } else {
938: DMStagStencil row, col[5];
939: PetscScalar val_A[5];
940: const PetscScalar val_rhs = 0.0;
942: row.i = ex;
943: row.j = ey;
944: row.loc = DMSTAG_ELEMENT;
945: row.c = 0;
946: col[0].i = ex;
947: col[0].j = ey;
948: col[0].loc = DMSTAG_LEFT;
949: col[0].c = 0;
950: val_A[0] = -1.0 * K_cont * dv / hx;
951: col[1].i = ex;
952: col[1].j = ey;
953: col[1].loc = DMSTAG_RIGHT;
954: col[1].c = 0;
955: val_A[1] = K_cont * dv / hx;
956: col[2].i = ex;
957: col[2].j = ey;
958: col[2].loc = DMSTAG_DOWN;
959: col[2].c = 0;
960: val_A[2] = -1.0 * K_cont * dv / hy;
961: col[3].i = ex;
962: col[3].j = ey;
963: col[3].loc = DMSTAG_UP;
964: col[3].c = 0;
965: val_A[3] = K_cont * dv / hy;
966: col[4] = row;
967: val_A[4] = 0.0;
968: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 5, col, val_A, INSERT_VALUES));
969: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
970: }
971: }
972: }
973: }
974: PetscCall(DMRestoreLocalVector(dm_coefficients, &coefficients_local));
976: /* Add additional inverse viscosity terms (for use in building a matrix used to construct the preconditioner) */
977: if (parameters->include_inverse_visc) {
978: PetscCheck(!parameters->faces_only, PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "Does not make sense with faces only");
979: PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_main, dm_coefficients, ctx->levels[level]->coeff, 1.0, A));
980: }
982: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
983: if (build_rhs) PetscCall(VecAssemblyBegin(rhs));
984: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
985: if (build_rhs) PetscCall(VecAssemblyEnd(rhs));
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: static PetscErrorCode CreateSystem3d(SystemParameters parameters, Mat *pA, Vec *pRhs)
990: {
991: PetscInt N[3];
992: PetscInt ex, ey, ez, startx, starty, startz, nx, ny, nz;
993: Mat A;
994: PetscReal hx, hy, hz, dv;
995: PetscInt pinx, piny, pinz;
996: Vec coeff_local, rhs;
997: PetscBool build_rhs;
998: DM dm_main, dm_coefficients;
999: PetscScalar K_cont, K_bound;
1000: Ctx ctx = parameters->ctx;
1001: PetscInt level = parameters->level;
1003: PetscFunctionBeginUser;
1004: if (parameters->faces_only) {
1005: dm_main = ctx->levels[level]->dm_faces;
1006: } else {
1007: dm_main = ctx->levels[level]->dm_stokes;
1008: }
1009: dm_coefficients = ctx->levels[level]->dm_coefficients;
1010: K_cont = ctx->levels[level]->K_cont;
1011: K_bound = ctx->levels[level]->K_bound;
1012: PetscCall(DMCreateMatrix(dm_main, pA));
1013: A = *pA;
1014: build_rhs = (PetscBool)(pRhs != NULL);
1015: if (build_rhs) {
1016: PetscCall(DMCreateGlobalVector(dm_main, pRhs));
1017: rhs = *pRhs;
1018: } else {
1019: rhs = NULL;
1020: }
1021: PetscCall(DMStagGetCorners(dm_main, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
1022: PetscCall(DMStagGetGlobalSizes(dm_main, &N[0], &N[1], &N[2]));
1023: hx = ctx->levels[level]->hx_characteristic;
1024: hy = ctx->levels[level]->hy_characteristic;
1025: hz = ctx->levels[level]->hz_characteristic;
1026: dv = hx * hy * hz;
1027: PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
1028: PetscCall(DMGlobalToLocal(dm_coefficients, ctx->levels[level]->coeff, INSERT_VALUES, coeff_local));
1030: 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");
1031: pinx = 1;
1032: piny = 0;
1033: pinz = 0; /* Depends on assertion above that there are at least two element in the x direction */
1035: /* Loop over all local elements.
1037: For each element, fill 4-7 rows of the matrix, corresponding to
1038: - the pressure degree of freedom (dof), centered on the element
1039: - the 3 velocity dofs on left, bottom, and back faces of the element
1040: - velocity dof on the right, top, and front faces of the element (only on domain boundaries)
1042: */
1043: for (ez = startz; ez < startz + nz; ++ez) {
1044: for (ey = starty; ey < starty + ny; ++ey) {
1045: for (ex = startx; ex < startx + nx; ++ex) {
1046: const PetscBool left_boundary = (PetscBool)(ex == 0);
1047: const PetscBool right_boundary = (PetscBool)(ex == N[0] - 1);
1048: const PetscBool bottom_boundary = (PetscBool)(ey == 0);
1049: const PetscBool top_boundary = (PetscBool)(ey == N[1] - 1);
1050: const PetscBool back_boundary = (PetscBool)(ez == 0);
1051: const PetscBool front_boundary = (PetscBool)(ez == N[2] - 1);
1053: /* Note that below, we depend on the check above that there is never one
1054: element (globally) in a given direction. Thus, for example, an
1055: element is never both on the left and right boundary */
1057: /* X-faces - right boundary */
1058: if (right_boundary) {
1059: /* Right x-velocity Dirichlet */
1060: DMStagStencil row;
1061: const PetscScalar val_rhs = 0.0;
1062: const PetscScalar val_A = K_bound;
1064: row.i = ex;
1065: row.j = ey;
1066: row.k = ez;
1067: row.loc = DMSTAG_RIGHT;
1068: row.c = 0;
1069: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1070: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1071: }
1073: /* X faces - left*/
1074: {
1075: DMStagStencil row;
1077: row.i = ex;
1078: row.j = ey;
1079: row.k = ez;
1080: row.loc = DMSTAG_LEFT;
1081: row.c = 0;
1083: if (left_boundary) {
1084: /* Left x-velocity Dirichlet */
1085: const PetscScalar val_rhs = 0.0;
1086: const PetscScalar val_A = K_bound;
1088: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1089: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1090: } else {
1091: /* X-momentum equation */
1092: PetscInt count;
1093: DMStagStencil col[17];
1094: PetscScalar val_A[17];
1095: DMStagStencil eta_point[6];
1096: PetscScalar eta[6], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the left face */
1097: const PetscScalar val_rhs = 0.0;
1099: /* Get eta values */
1100: eta_point[0].i = ex - 1;
1101: eta_point[0].j = ey;
1102: eta_point[0].k = ez;
1103: eta_point[0].loc = DMSTAG_ELEMENT;
1104: eta_point[0].c = 0; /* Left */
1105: eta_point[1].i = ex;
1106: eta_point[1].j = ey;
1107: eta_point[1].k = ez;
1108: eta_point[1].loc = DMSTAG_ELEMENT;
1109: eta_point[1].c = 0; /* Right */
1110: eta_point[2].i = ex;
1111: eta_point[2].j = ey;
1112: eta_point[2].k = ez;
1113: eta_point[2].loc = DMSTAG_DOWN_LEFT;
1114: eta_point[2].c = 0; /* Down */
1115: eta_point[3].i = ex;
1116: eta_point[3].j = ey;
1117: eta_point[3].k = ez;
1118: eta_point[3].loc = DMSTAG_UP_LEFT;
1119: eta_point[3].c = 0; /* Up */
1120: eta_point[4].i = ex;
1121: eta_point[4].j = ey;
1122: eta_point[4].k = ez;
1123: eta_point[4].loc = DMSTAG_BACK_LEFT;
1124: eta_point[4].c = 0; /* Back */
1125: eta_point[5].i = ex;
1126: eta_point[5].j = ey;
1127: eta_point[5].k = ez;
1128: eta_point[5].loc = DMSTAG_FRONT_LEFT;
1129: eta_point[5].c = 0; /* Front */
1130: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
1131: eta_left = eta[0];
1132: eta_right = eta[1];
1133: eta_down = eta[2];
1134: eta_up = eta[3];
1135: eta_back = eta[4];
1136: eta_front = eta[5];
1138: count = 0;
1140: col[count] = row;
1141: val_A[count] = -2.0 * dv * (eta_left + eta_right) / (hx * hx);
1142: if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
1143: if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
1144: if (!back_boundary) val_A[count] += -1.0 * dv * eta_back / (hz * hz);
1145: if (!front_boundary) val_A[count] += -1.0 * dv * eta_front / (hz * hz);
1146: ++count;
1148: col[count].i = ex - 1;
1149: col[count].j = ey;
1150: col[count].k = ez;
1151: col[count].loc = DMSTAG_LEFT;
1152: col[count].c = 0;
1153: val_A[count] = 2.0 * dv * eta_left / (hx * hx);
1154: ++count;
1155: col[count].i = ex + 1;
1156: col[count].j = ey;
1157: col[count].k = ez;
1158: col[count].loc = DMSTAG_LEFT;
1159: col[count].c = 0;
1160: val_A[count] = 2.0 * dv * eta_right / (hx * hx);
1161: ++count;
1162: if (!bottom_boundary) {
1163: col[count].i = ex;
1164: col[count].j = ey - 1;
1165: col[count].k = ez;
1166: col[count].loc = DMSTAG_LEFT;
1167: col[count].c = 0;
1168: val_A[count] = dv * eta_down / (hy * hy);
1169: ++count;
1170: }
1171: if (!top_boundary) {
1172: col[count].i = ex;
1173: col[count].j = ey + 1;
1174: col[count].k = ez;
1175: col[count].loc = DMSTAG_LEFT;
1176: col[count].c = 0;
1177: val_A[count] = dv * eta_up / (hy * hy);
1178: ++count;
1179: }
1180: if (!back_boundary) {
1181: col[count].i = ex;
1182: col[count].j = ey;
1183: col[count].k = ez - 1;
1184: col[count].loc = DMSTAG_LEFT;
1185: col[count].c = 0;
1186: val_A[count] = dv * eta_back / (hz * hz);
1187: ++count;
1188: }
1189: if (!front_boundary) {
1190: col[count].i = ex;
1191: col[count].j = ey;
1192: col[count].k = ez + 1;
1193: col[count].loc = DMSTAG_LEFT;
1194: col[count].c = 0;
1195: val_A[count] = dv * eta_front / (hz * hz);
1196: ++count;
1197: }
1199: col[count].i = ex - 1;
1200: col[count].j = ey;
1201: col[count].k = ez;
1202: col[count].loc = DMSTAG_DOWN;
1203: col[count].c = 0;
1204: val_A[count] = dv * eta_down / (hx * hy);
1205: ++count; /* down left */
1206: col[count].i = ex;
1207: col[count].j = ey;
1208: col[count].k = ez;
1209: col[count].loc = DMSTAG_DOWN;
1210: col[count].c = 0;
1211: val_A[count] = -1.0 * dv * eta_down / (hx * hy);
1212: ++count; /* down right */
1214: col[count].i = ex - 1;
1215: col[count].j = ey;
1216: col[count].k = ez;
1217: col[count].loc = DMSTAG_UP;
1218: col[count].c = 0;
1219: val_A[count] = -1.0 * dv * eta_up / (hx * hy);
1220: ++count; /* up left */
1221: col[count].i = ex;
1222: col[count].j = ey;
1223: col[count].k = ez;
1224: col[count].loc = DMSTAG_UP;
1225: col[count].c = 0;
1226: val_A[count] = dv * eta_up / (hx * hy);
1227: ++count; /* up right */
1229: col[count].i = ex - 1;
1230: col[count].j = ey;
1231: col[count].k = ez;
1232: col[count].loc = DMSTAG_BACK;
1233: col[count].c = 0;
1234: val_A[count] = dv * eta_back / (hx * hz);
1235: ++count; /* back left */
1236: col[count].i = ex;
1237: col[count].j = ey;
1238: col[count].k = ez;
1239: col[count].loc = DMSTAG_BACK;
1240: col[count].c = 0;
1241: val_A[count] = -1.0 * dv * eta_back / (hx * hz);
1242: ++count; /* back right */
1244: col[count].i = ex - 1;
1245: col[count].j = ey;
1246: col[count].k = ez;
1247: col[count].loc = DMSTAG_FRONT;
1248: col[count].c = 0;
1249: val_A[count] = -1.0 * dv * eta_front / (hx * hz);
1250: ++count; /* front left */
1251: col[count].i = ex;
1252: col[count].j = ey;
1253: col[count].k = ez;
1254: col[count].loc = DMSTAG_FRONT;
1255: col[count].c = 0;
1256: val_A[count] = dv * eta_front / (hx * hz);
1257: ++count; /* front right */
1259: if (!parameters->faces_only) {
1260: col[count].i = ex - 1;
1261: col[count].j = ey;
1262: col[count].k = ez;
1263: col[count].loc = DMSTAG_ELEMENT;
1264: col[count].c = 0;
1265: val_A[count] = K_cont * dv / hx;
1266: ++count;
1267: col[count].i = ex;
1268: col[count].j = ey;
1269: col[count].k = ez;
1270: col[count].loc = DMSTAG_ELEMENT;
1271: col[count].c = 0;
1272: val_A[count] = -1.0 * K_cont * dv / hx;
1273: ++count;
1274: }
1276: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
1277: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1278: }
1279: }
1281: /* Y faces - top boundary */
1282: if (top_boundary) {
1283: /* Top y-velocity Dirichlet */
1284: DMStagStencil row;
1285: const PetscScalar val_rhs = 0.0;
1286: const PetscScalar val_A = K_bound;
1288: row.i = ex;
1289: row.j = ey;
1290: row.k = ez;
1291: row.loc = DMSTAG_UP;
1292: row.c = 0;
1293: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1294: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1295: }
1297: /* Y faces - down */
1298: {
1299: DMStagStencil row;
1301: row.i = ex;
1302: row.j = ey;
1303: row.k = ez;
1304: row.loc = DMSTAG_DOWN;
1305: row.c = 0;
1307: if (bottom_boundary) {
1308: /* Bottom y-velocity Dirichlet */
1309: const PetscScalar val_rhs = 0.0;
1310: const PetscScalar val_A = K_bound;
1312: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1313: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1314: } else {
1315: /* Y-momentum equation (including non-zero forcing) */
1316: PetscInt count;
1317: DMStagStencil col[17];
1318: PetscScalar val_rhs, val_A[17];
1319: DMStagStencil eta_point[6], rho_point[4];
1320: PetscScalar eta[6], rho[4], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the bottom face */
1322: if (build_rhs) {
1323: /* get rho values (note .c = 1) */
1324: /* Note that we have rho at perhaps strange points (edges not corners) */
1325: rho_point[0].i = ex;
1326: rho_point[0].j = ey;
1327: rho_point[0].k = ez;
1328: rho_point[0].loc = DMSTAG_DOWN_LEFT;
1329: rho_point[0].c = 1;
1330: rho_point[1].i = ex;
1331: rho_point[1].j = ey;
1332: rho_point[1].k = ez;
1333: rho_point[1].loc = DMSTAG_DOWN_RIGHT;
1334: rho_point[1].c = 1;
1335: rho_point[2].i = ex;
1336: rho_point[2].j = ey;
1337: rho_point[2].k = ez;
1338: rho_point[2].loc = DMSTAG_BACK_DOWN;
1339: rho_point[2].c = 1;
1340: rho_point[3].i = ex;
1341: rho_point[3].j = ey;
1342: rho_point[3].k = ez;
1343: rho_point[3].loc = DMSTAG_FRONT_DOWN;
1344: rho_point[3].c = 1;
1345: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 4, rho_point, rho));
1347: /* Compute forcing */
1348: val_rhs = ctx->gy * dv * (0.25 * (rho[0] + rho[1] + rho[2] + rho[3]));
1349: }
1351: /* Get eta values */
1352: eta_point[0].i = ex;
1353: eta_point[0].j = ey;
1354: eta_point[0].k = ez;
1355: eta_point[0].loc = DMSTAG_DOWN_LEFT;
1356: eta_point[0].c = 0; /* Left */
1357: eta_point[1].i = ex;
1358: eta_point[1].j = ey;
1359: eta_point[1].k = ez;
1360: eta_point[1].loc = DMSTAG_DOWN_RIGHT;
1361: eta_point[1].c = 0; /* Right */
1362: eta_point[2].i = ex;
1363: eta_point[2].j = ey - 1;
1364: eta_point[2].k = ez;
1365: eta_point[2].loc = DMSTAG_ELEMENT;
1366: eta_point[2].c = 0; /* Down */
1367: eta_point[3].i = ex;
1368: eta_point[3].j = ey;
1369: eta_point[3].k = ez;
1370: eta_point[3].loc = DMSTAG_ELEMENT;
1371: eta_point[3].c = 0; /* Up */
1372: eta_point[4].i = ex;
1373: eta_point[4].j = ey;
1374: eta_point[4].k = ez;
1375: eta_point[4].loc = DMSTAG_BACK_DOWN;
1376: eta_point[4].c = 0; /* Back */
1377: eta_point[5].i = ex;
1378: eta_point[5].j = ey;
1379: eta_point[5].k = ez;
1380: eta_point[5].loc = DMSTAG_FRONT_DOWN;
1381: eta_point[5].c = 0; /* Front */
1382: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
1383: eta_left = eta[0];
1384: eta_right = eta[1];
1385: eta_down = eta[2];
1386: eta_up = eta[3];
1387: eta_back = eta[4];
1388: eta_front = eta[5];
1390: count = 0;
1392: col[count] = row;
1393: val_A[count] = -2.0 * dv * (eta_up + eta_down) / (hy * hy);
1394: if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
1395: if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
1396: if (!back_boundary) val_A[count] += -1.0 * dv * eta_back / (hz * hz);
1397: if (!front_boundary) val_A[count] += -1.0 * dv * eta_front / (hz * hz);
1398: ++count;
1400: col[count].i = ex;
1401: col[count].j = ey - 1;
1402: col[count].k = ez;
1403: col[count].loc = DMSTAG_DOWN;
1404: col[count].c = 0;
1405: val_A[count] = 2.0 * dv * eta_down / (hy * hy);
1406: ++count;
1407: col[count].i = ex;
1408: col[count].j = ey + 1;
1409: col[count].k = ez;
1410: col[count].loc = DMSTAG_DOWN;
1411: col[count].c = 0;
1412: val_A[count] = 2.0 * dv * eta_up / (hy * hy);
1413: ++count;
1415: if (!left_boundary) {
1416: col[count].i = ex - 1;
1417: col[count].j = ey;
1418: col[count].k = ez;
1419: col[count].loc = DMSTAG_DOWN;
1420: col[count].c = 0;
1421: val_A[count] = dv * eta_left / (hx * hx);
1422: ++count;
1423: }
1424: if (!right_boundary) {
1425: col[count].i = ex + 1;
1426: col[count].j = ey;
1427: col[count].k = ez;
1428: col[count].loc = DMSTAG_DOWN;
1429: col[count].c = 0;
1430: val_A[count] = dv * eta_right / (hx * hx);
1431: ++count;
1432: }
1433: if (!back_boundary) {
1434: col[count].i = ex;
1435: col[count].j = ey;
1436: col[count].k = ez - 1;
1437: col[count].loc = DMSTAG_DOWN;
1438: col[count].c = 0;
1439: val_A[count] = dv * eta_back / (hz * hz);
1440: ++count;
1441: }
1442: if (!front_boundary) {
1443: col[count].i = ex;
1444: col[count].j = ey;
1445: col[count].k = ez + 1;
1446: col[count].loc = DMSTAG_DOWN;
1447: col[count].c = 0;
1448: val_A[count] = dv * eta_front / (hz * hz);
1449: ++count;
1450: }
1452: col[count].i = ex;
1453: col[count].j = ey - 1;
1454: col[count].k = ez;
1455: col[count].loc = DMSTAG_LEFT;
1456: col[count].c = 0;
1457: val_A[count] = dv * eta_left / (hx * hy);
1458: ++count; /* down left*/
1459: col[count].i = ex;
1460: col[count].j = ey;
1461: col[count].k = ez;
1462: col[count].loc = DMSTAG_LEFT;
1463: col[count].c = 0;
1464: val_A[count] = -1.0 * dv * eta_left / (hx * hy);
1465: ++count; /* up left*/
1467: col[count].i = ex;
1468: col[count].j = ey - 1;
1469: col[count].k = ez;
1470: col[count].loc = DMSTAG_RIGHT;
1471: col[count].c = 0;
1472: val_A[count] = -1.0 * dv * eta_right / (hx * hy);
1473: ++count; /* down right*/
1474: col[count].i = ex;
1475: col[count].j = ey;
1476: col[count].k = ez;
1477: col[count].loc = DMSTAG_RIGHT;
1478: col[count].c = 0;
1479: val_A[count] = dv * eta_right / (hx * hy);
1480: ++count; /* up right*/
1482: col[count].i = ex;
1483: col[count].j = ey - 1;
1484: col[count].k = ez;
1485: col[count].loc = DMSTAG_BACK;
1486: col[count].c = 0;
1487: val_A[count] = dv * eta_back / (hy * hz);
1488: ++count; /* back down */
1489: col[count].i = ex;
1490: col[count].j = ey;
1491: col[count].k = ez;
1492: col[count].loc = DMSTAG_BACK;
1493: col[count].c = 0;
1494: val_A[count] = -1.0 * dv * eta_back / (hy * hz);
1495: ++count; /* back up */
1497: col[count].i = ex;
1498: col[count].j = ey - 1;
1499: col[count].k = ez;
1500: col[count].loc = DMSTAG_FRONT;
1501: col[count].c = 0;
1502: val_A[count] = -1.0 * dv * eta_front / (hy * hz);
1503: ++count; /* front down */
1504: col[count].i = ex;
1505: col[count].j = ey;
1506: col[count].k = ez;
1507: col[count].loc = DMSTAG_FRONT;
1508: col[count].c = 0;
1509: val_A[count] = dv * eta_front / (hy * hz);
1510: ++count; /* front up */
1512: if (!parameters->faces_only) {
1513: col[count].i = ex;
1514: col[count].j = ey - 1;
1515: col[count].k = ez;
1516: col[count].loc = DMSTAG_ELEMENT;
1517: col[count].c = 0;
1518: val_A[count] = K_cont * dv / hy;
1519: ++count;
1520: col[count].i = ex;
1521: col[count].j = ey;
1522: col[count].k = ez;
1523: col[count].loc = DMSTAG_ELEMENT;
1524: col[count].c = 0;
1525: val_A[count] = -1.0 * K_cont * dv / hy;
1526: ++count;
1527: }
1529: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
1530: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1531: }
1532: }
1534: if (front_boundary) {
1535: /* Front z-velocity Dirichlet */
1536: DMStagStencil row;
1537: const PetscScalar val_rhs = 0.0;
1538: const PetscScalar val_A = K_bound;
1540: row.i = ex;
1541: row.j = ey;
1542: row.k = ez;
1543: row.loc = DMSTAG_FRONT;
1544: row.c = 0;
1545: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1546: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1547: }
1549: /* Z faces - back */
1550: {
1551: DMStagStencil row;
1553: row.i = ex;
1554: row.j = ey;
1555: row.k = ez;
1556: row.loc = DMSTAG_BACK;
1557: row.c = 0;
1559: if (back_boundary) {
1560: /* Back z-velocity Dirichlet */
1561: const PetscScalar val_rhs = 0.0;
1562: const PetscScalar val_A = K_bound;
1564: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1565: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1566: } else {
1567: /* Z-momentum equation */
1568: PetscInt count;
1569: DMStagStencil col[17];
1570: PetscScalar val_A[17];
1571: DMStagStencil eta_point[6];
1572: PetscScalar eta[6], eta_left, eta_right, eta_up, eta_down, eta_back, eta_front; /* relative to the back face */
1573: const PetscScalar val_rhs = 0.0;
1575: /* Get eta values */
1576: eta_point[0].i = ex;
1577: eta_point[0].j = ey;
1578: eta_point[0].k = ez;
1579: eta_point[0].loc = DMSTAG_BACK_LEFT;
1580: eta_point[0].c = 0; /* Left */
1581: eta_point[1].i = ex;
1582: eta_point[1].j = ey;
1583: eta_point[1].k = ez;
1584: eta_point[1].loc = DMSTAG_BACK_RIGHT;
1585: eta_point[1].c = 0; /* Right */
1586: eta_point[2].i = ex;
1587: eta_point[2].j = ey;
1588: eta_point[2].k = ez;
1589: eta_point[2].loc = DMSTAG_BACK_DOWN;
1590: eta_point[2].c = 0; /* Down */
1591: eta_point[3].i = ex;
1592: eta_point[3].j = ey;
1593: eta_point[3].k = ez;
1594: eta_point[3].loc = DMSTAG_BACK_UP;
1595: eta_point[3].c = 0; /* Up */
1596: eta_point[4].i = ex;
1597: eta_point[4].j = ey;
1598: eta_point[4].k = ez - 1;
1599: eta_point[4].loc = DMSTAG_ELEMENT;
1600: eta_point[4].c = 0; /* Back */
1601: eta_point[5].i = ex;
1602: eta_point[5].j = ey;
1603: eta_point[5].k = ez;
1604: eta_point[5].loc = DMSTAG_ELEMENT;
1605: eta_point[5].c = 0; /* Front */
1606: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 6, eta_point, eta));
1607: eta_left = eta[0];
1608: eta_right = eta[1];
1609: eta_down = eta[2];
1610: eta_up = eta[3];
1611: eta_back = eta[4];
1612: eta_front = eta[5];
1614: count = 0;
1616: col[count] = row;
1617: val_A[count] = -2.0 * dv * (eta_back + eta_front) / (hz * hz);
1618: if (!left_boundary) val_A[count] += -1.0 * dv * eta_left / (hx * hx);
1619: if (!right_boundary) val_A[count] += -1.0 * dv * eta_right / (hx * hx);
1620: if (!top_boundary) val_A[count] += -1.0 * dv * eta_up / (hy * hy);
1621: if (!bottom_boundary) val_A[count] += -1.0 * dv * eta_down / (hy * hy);
1622: ++count;
1624: col[count].i = ex;
1625: col[count].j = ey;
1626: col[count].k = ez - 1;
1627: col[count].loc = DMSTAG_BACK;
1628: col[count].c = 0;
1629: val_A[count] = 2.0 * dv * eta_back / (hz * hz);
1630: ++count;
1631: col[count].i = ex;
1632: col[count].j = ey;
1633: col[count].k = ez + 1;
1634: col[count].loc = DMSTAG_BACK;
1635: col[count].c = 0;
1636: val_A[count] = 2.0 * dv * eta_front / (hz * hz);
1637: ++count;
1639: if (!left_boundary) {
1640: col[count].i = ex - 1;
1641: col[count].j = ey;
1642: col[count].k = ez;
1643: col[count].loc = DMSTAG_BACK;
1644: col[count].c = 0;
1645: val_A[count] = dv * eta_left / (hx * hx);
1646: ++count;
1647: }
1648: if (!right_boundary) {
1649: col[count].i = ex + 1;
1650: col[count].j = ey;
1651: col[count].k = ez;
1652: col[count].loc = DMSTAG_BACK;
1653: col[count].c = 0;
1654: val_A[count] = dv * eta_right / (hx * hx);
1655: ++count;
1656: }
1657: if (!bottom_boundary) {
1658: col[count].i = ex;
1659: col[count].j = ey - 1;
1660: col[count].k = ez;
1661: col[count].loc = DMSTAG_BACK;
1662: col[count].c = 0;
1663: val_A[count] = dv * eta_down / (hy * hy);
1664: ++count;
1665: }
1666: if (!top_boundary) {
1667: col[count].i = ex;
1668: col[count].j = ey + 1;
1669: col[count].k = ez;
1670: col[count].loc = DMSTAG_BACK;
1671: col[count].c = 0;
1672: val_A[count] = dv * eta_up / (hy * hy);
1673: ++count;
1674: }
1676: col[count].i = ex;
1677: col[count].j = ey;
1678: col[count].k = ez - 1;
1679: col[count].loc = DMSTAG_LEFT;
1680: col[count].c = 0;
1681: val_A[count] = dv * eta_left / (hx * hz);
1682: ++count; /* back left*/
1683: col[count].i = ex;
1684: col[count].j = ey;
1685: col[count].k = ez;
1686: col[count].loc = DMSTAG_LEFT;
1687: col[count].c = 0;
1688: val_A[count] = -1.0 * dv * eta_left / (hx * hz);
1689: ++count; /* front left*/
1691: col[count].i = ex;
1692: col[count].j = ey;
1693: col[count].k = ez - 1;
1694: col[count].loc = DMSTAG_RIGHT;
1695: col[count].c = 0;
1696: val_A[count] = -1.0 * dv * eta_right / (hx * hz);
1697: ++count; /* back right */
1698: col[count].i = ex;
1699: col[count].j = ey;
1700: col[count].k = ez;
1701: col[count].loc = DMSTAG_RIGHT;
1702: col[count].c = 0;
1703: val_A[count] = dv * eta_right / (hx * hz);
1704: ++count; /* front right*/
1706: col[count].i = ex;
1707: col[count].j = ey;
1708: col[count].k = ez - 1;
1709: col[count].loc = DMSTAG_DOWN;
1710: col[count].c = 0;
1711: val_A[count] = dv * eta_down / (hy * hz);
1712: ++count; /* back down */
1713: col[count].i = ex;
1714: col[count].j = ey;
1715: col[count].k = ez;
1716: col[count].loc = DMSTAG_DOWN;
1717: col[count].c = 0;
1718: val_A[count] = -1.0 * dv * eta_down / (hy * hz);
1719: ++count; /* back down */
1721: col[count].i = ex;
1722: col[count].j = ey;
1723: col[count].k = ez - 1;
1724: col[count].loc = DMSTAG_UP;
1725: col[count].c = 0;
1726: val_A[count] = -1.0 * dv * eta_up / (hy * hz);
1727: ++count; /* back up */
1728: col[count].i = ex;
1729: col[count].j = ey;
1730: col[count].k = ez;
1731: col[count].loc = DMSTAG_UP;
1732: col[count].c = 0;
1733: val_A[count] = dv * eta_up / (hy * hz);
1734: ++count; /* back up */
1736: if (!parameters->faces_only) {
1737: col[count].i = ex;
1738: col[count].j = ey;
1739: col[count].k = ez - 1;
1740: col[count].loc = DMSTAG_ELEMENT;
1741: col[count].c = 0;
1742: val_A[count] = K_cont * dv / hz;
1743: ++count;
1744: col[count].i = ex;
1745: col[count].j = ey;
1746: col[count].k = ez;
1747: col[count].loc = DMSTAG_ELEMENT;
1748: col[count].c = 0;
1749: val_A[count] = -1.0 * K_cont * dv / hz;
1750: ++count;
1751: }
1753: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, count, col, val_A, INSERT_VALUES));
1754: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1755: }
1756: }
1758: /* Elements */
1759: if (!parameters->faces_only) {
1760: DMStagStencil row;
1762: row.i = ex;
1763: row.j = ey;
1764: row.k = ez;
1765: row.loc = DMSTAG_ELEMENT;
1766: row.c = 0;
1768: if (ctx->pin_pressure && ex == pinx && ey == piny && ez == pinz) {
1769: /* Pin a pressure node to zero, if requested */
1770: const PetscScalar val_A = K_bound;
1771: const PetscScalar val_rhs = 0.0;
1773: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 1, &row, &val_A, INSERT_VALUES));
1774: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1775: } else {
1776: /* Continuity equation */
1777: /* Note that this includes an explicit zero on the diagonal. This is only needed for
1778: some direct solvers (not required if using an iterative solver and setting a constant-pressure nullspace) */
1779: /* Note: the scaling by dv is not chosen in a principled way and is likely sub-optimal */
1780: DMStagStencil col[7];
1781: PetscScalar val_A[7];
1782: const PetscScalar val_rhs = 0.0;
1784: col[0].i = ex;
1785: col[0].j = ey;
1786: col[0].k = ez;
1787: col[0].loc = DMSTAG_LEFT;
1788: col[0].c = 0;
1789: val_A[0] = -1.0 * K_cont * dv / hx;
1790: col[1].i = ex;
1791: col[1].j = ey;
1792: col[1].k = ez;
1793: col[1].loc = DMSTAG_RIGHT;
1794: col[1].c = 0;
1795: val_A[1] = K_cont * dv / hx;
1796: col[2].i = ex;
1797: col[2].j = ey;
1798: col[2].k = ez;
1799: col[2].loc = DMSTAG_DOWN;
1800: col[2].c = 0;
1801: val_A[2] = -1.0 * K_cont * dv / hy;
1802: col[3].i = ex;
1803: col[3].j = ey;
1804: col[3].k = ez;
1805: col[3].loc = DMSTAG_UP;
1806: col[3].c = 0;
1807: val_A[3] = K_cont * dv / hy;
1808: col[4].i = ex;
1809: col[4].j = ey;
1810: col[4].k = ez;
1811: col[4].loc = DMSTAG_BACK;
1812: col[4].c = 0;
1813: val_A[4] = -1.0 * K_cont * dv / hz;
1814: col[5].i = ex;
1815: col[5].j = ey;
1816: col[5].k = ez;
1817: col[5].loc = DMSTAG_FRONT;
1818: col[5].c = 0;
1819: val_A[5] = K_cont * dv / hz;
1820: col[6] = row;
1821: val_A[6] = 0.0;
1822: PetscCall(DMStagMatSetValuesStencil(dm_main, A, 1, &row, 7, col, val_A, INSERT_VALUES));
1823: if (build_rhs) PetscCall(DMStagVecSetValuesStencil(dm_main, rhs, 1, &row, &val_rhs, INSERT_VALUES));
1824: }
1825: }
1826: }
1827: }
1828: }
1829: PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));
1831: /* Add additional inverse viscosity terms (for use in building a matrix used to construct the preconditioner) */
1832: if (parameters->include_inverse_visc) {
1833: PetscCheck(!parameters->faces_only, PetscObjectComm((PetscObject)dm_main), PETSC_ERR_SUP, "Does not make sense with faces only");
1834: PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_main, dm_coefficients, ctx->levels[level]->coeff, 1.0, A));
1835: }
1837: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1838: if (build_rhs) PetscCall(VecAssemblyBegin(rhs));
1839: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1840: if (build_rhs) PetscCall(VecAssemblyEnd(rhs));
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: static PetscErrorCode CreateSystem(SystemParameters parameters, Mat *pA, Vec *pRhs)
1845: {
1846: PetscFunctionBeginUser;
1847: if (parameters->ctx->dim == 2) {
1848: PetscCall(CreateSystem2d(parameters, pA, pRhs));
1849: } else if (parameters->ctx->dim == 3) {
1850: PetscCall(CreateSystem3d(parameters, pA, pRhs));
1851: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, parameters->ctx->dim);
1852: PetscFunctionReturn(PETSC_SUCCESS);
1853: }
1855: PetscErrorCode PopulateCoefficientData(Ctx ctx, PetscInt level)
1856: {
1857: PetscInt dim;
1858: PetscInt N[3];
1859: PetscInt ex, ey, ez, startx, starty, startz, nx, ny, nz;
1860: PetscInt slot_prev, slot_center;
1861: PetscInt slot_rho_downleft, slot_rho_backleft, slot_rho_backdown, slot_eta_element, slot_eta_downleft, slot_eta_backleft, slot_eta_backdown;
1862: Vec coeff_local;
1863: PetscReal **arr_coordinates_x, **arr_coordinates_y, **arr_coordinates_z;
1864: DM dm_coefficients;
1865: Vec coeff;
1867: PetscFunctionBeginUser;
1868: dm_coefficients = ctx->levels[level]->dm_coefficients;
1869: PetscCall(DMGetDimension(dm_coefficients, &dim));
1871: /* Create global coefficient vector */
1872: PetscCall(DMCreateGlobalVector(dm_coefficients, &ctx->levels[level]->coeff));
1873: coeff = ctx->levels[level]->coeff;
1875: /* Get temporary access to a local representation of the coefficient data */
1876: PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
1877: PetscCall(DMGlobalToLocal(dm_coefficients, coeff, INSERT_VALUES, coeff_local));
1879: /* Use direct array access to coefficient and coordinate arrays, to popoulate coefficient data */
1880: PetscCall(DMStagGetGhostCorners(dm_coefficients, &startx, &starty, &startz, &nx, &ny, &nz));
1881: PetscCall(DMStagGetGlobalSizes(dm_coefficients, &N[0], &N[1], &N[2]));
1882: PetscCall(DMStagGetProductCoordinateArraysRead(dm_coefficients, &arr_coordinates_x, &arr_coordinates_y, &arr_coordinates_z));
1883: PetscCall(DMStagGetProductCoordinateLocationSlot(dm_coefficients, DMSTAG_ELEMENT, &slot_center));
1884: PetscCall(DMStagGetProductCoordinateLocationSlot(dm_coefficients, DMSTAG_LEFT, &slot_prev));
1885: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_ELEMENT, 0, &slot_eta_element));
1886: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_DOWN_LEFT, 0, &slot_eta_downleft));
1887: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_DOWN_LEFT, 1, &slot_rho_downleft));
1888: if (dim == 2) {
1889: PetscScalar ***arr_coefficients;
1891: PetscCall(DMStagVecGetArray(dm_coefficients, coeff_local, &arr_coefficients));
1892: /* Note that these ranges are with respect to the local representation */
1893: for (ey = starty; ey < starty + ny; ++ey) {
1894: for (ex = startx; ex < startx + nx; ++ex) {
1895: arr_coefficients[ey][ex][slot_eta_element] = ctx->GetEta(ctx, arr_coordinates_x[ex][slot_center], arr_coordinates_y[ey][slot_center], 0.0);
1896: arr_coefficients[ey][ex][slot_eta_downleft] = ctx->GetEta(ctx, arr_coordinates_x[ex][slot_prev], arr_coordinates_y[ey][slot_prev], 0.0);
1897: arr_coefficients[ey][ex][slot_rho_downleft] = ctx->GetRho(ctx, arr_coordinates_x[ex][slot_prev], arr_coordinates_y[ey][slot_prev], 0.0);
1898: }
1899: }
1900: PetscCall(DMStagVecRestoreArray(dm_coefficients, coeff_local, &arr_coefficients));
1901: } else if (dim == 3) {
1902: PetscScalar ****arr_coefficients;
1904: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_LEFT, 0, &slot_eta_backleft));
1905: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_LEFT, 1, &slot_rho_backleft));
1906: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_DOWN, 0, &slot_eta_backdown));
1907: PetscCall(DMStagGetLocationSlot(dm_coefficients, DMSTAG_BACK_DOWN, 1, &slot_rho_backdown));
1908: PetscCall(DMStagVecGetArray(dm_coefficients, coeff_local, &arr_coefficients));
1909: /* Note that these are with respect to the entire local representation, including ghosts */
1910: for (ez = startz; ez < startz + nz; ++ez) {
1911: for (ey = starty; ey < starty + ny; ++ey) {
1912: for (ex = startx; ex < startx + nx; ++ex) {
1913: const PetscScalar x_prev = arr_coordinates_x[ex][slot_prev];
1914: const PetscScalar y_prev = arr_coordinates_y[ey][slot_prev];
1915: const PetscScalar z_prev = arr_coordinates_z[ez][slot_prev];
1916: const PetscScalar x_center = arr_coordinates_x[ex][slot_center];
1917: const PetscScalar y_center = arr_coordinates_y[ey][slot_center];
1918: const PetscScalar z_center = arr_coordinates_z[ez][slot_center];
1920: arr_coefficients[ez][ey][ex][slot_eta_element] = ctx->GetEta(ctx, x_center, y_center, z_center);
1921: arr_coefficients[ez][ey][ex][slot_eta_downleft] = ctx->GetEta(ctx, x_prev, y_prev, z_center);
1922: arr_coefficients[ez][ey][ex][slot_rho_downleft] = ctx->GetRho(ctx, x_prev, y_prev, z_center);
1923: arr_coefficients[ez][ey][ex][slot_eta_backleft] = ctx->GetEta(ctx, x_prev, y_center, z_prev);
1924: arr_coefficients[ez][ey][ex][slot_rho_backleft] = ctx->GetRho(ctx, x_prev, y_center, z_prev);
1925: arr_coefficients[ez][ey][ex][slot_eta_backdown] = ctx->GetEta(ctx, x_center, y_prev, z_prev);
1926: arr_coefficients[ez][ey][ex][slot_rho_backdown] = ctx->GetRho(ctx, x_center, y_prev, z_prev);
1927: }
1928: }
1929: }
1930: PetscCall(DMStagVecRestoreArray(dm_coefficients, coeff_local, &arr_coefficients));
1931: } else SETERRQ(PetscObjectComm((PetscObject)dm_coefficients), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
1932: PetscCall(DMStagRestoreProductCoordinateArraysRead(dm_coefficients, &arr_coordinates_x, &arr_coordinates_y, &arr_coordinates_z));
1933: PetscCall(DMLocalToGlobal(dm_coefficients, coeff_local, INSERT_VALUES, coeff));
1934: PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: static PetscErrorCode CreateAuxiliaryOperator(Ctx ctx, PetscInt level, Mat *p_S_hat)
1939: {
1940: DM dm_element;
1941: Mat S_hat;
1942: DM dm_stokes, dm_coefficients;
1943: Vec coeff;
1945: PetscFunctionBeginUser;
1946: dm_stokes = ctx->levels[level]->dm_stokes;
1947: dm_coefficients = ctx->levels[level]->dm_coefficients;
1948: coeff = ctx->levels[level]->coeff;
1949: if (ctx->dim == 2) {
1950: PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 1, 0, &dm_element));
1951: } else if (ctx->dim == 3) {
1952: PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 0, 1, &dm_element));
1953: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented for dimension %" PetscInt_FMT, ctx->dim);
1954: PetscCall(DMCreateMatrix(dm_element, p_S_hat));
1955: S_hat = *p_S_hat;
1956: PetscCall(OperatorInsertInverseViscosityPressureTerms(dm_element, dm_coefficients, coeff, 1.0, S_hat));
1957: PetscCall(MatAssemblyBegin(S_hat, MAT_FINAL_ASSEMBLY));
1958: PetscCall(MatAssemblyEnd(S_hat, MAT_FINAL_ASSEMBLY));
1959: PetscCall(DMDestroy(&dm_element));
1960: PetscFunctionReturn(PETSC_SUCCESS);
1961: }
1963: static PetscErrorCode OperatorInsertInverseViscosityPressureTerms(DM dm, DM dm_coefficients, Vec coefficients, PetscScalar scale, Mat mat)
1964: {
1965: PetscInt dim, ex, ey, ez, startx, starty, startz, nx, ny, nz;
1966: Vec coeff_local;
1968: PetscFunctionBeginUser;
1969: PetscCall(DMGetDimension(dm, &dim));
1970: PetscCall(DMGetLocalVector(dm_coefficients, &coeff_local));
1971: PetscCall(DMGlobalToLocal(dm_coefficients, coefficients, INSERT_VALUES, coeff_local));
1972: PetscCall(DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
1973: if (dim == 2) { /* Trick to have one loop nest */
1974: startz = 0;
1975: nz = 1;
1976: }
1977: for (ez = startz; ez < startz + nz; ++ez) {
1978: for (ey = starty; ey < starty + ny; ++ey) {
1979: for (ex = startx; ex < startx + nx; ++ex) {
1980: DMStagStencil from, to;
1981: PetscScalar val;
1983: /* component 0 on element is viscosity */
1984: from.i = ex;
1985: from.j = ey;
1986: from.k = ez;
1987: from.c = 0;
1988: from.loc = DMSTAG_ELEMENT;
1989: PetscCall(DMStagVecGetValuesStencil(dm_coefficients, coeff_local, 1, &from, &val));
1990: val = scale / val; /* inverse viscosity, scaled */
1991: to = from;
1992: PetscCall(DMStagMatSetValuesStencil(dm, mat, 1, &to, 1, &to, &val, INSERT_VALUES));
1993: }
1994: }
1995: }
1996: PetscCall(DMRestoreLocalVector(dm_coefficients, &coeff_local));
1997: /* Note that this function does not call MatAssembly{Begin,End} */
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: /* Create a pressure-only DMStag and use it to generate a nullspace vector
2002: - Create a compatible DMStag with one dof per element (and nothing else).
2003: - Create a constant vector and normalize it
2004: - Migrate it to a vector on the original dmSol (making use of the fact
2005: that this will fill in zeros for "extra" dof)
2006: - Set the nullspace for the operator
2007: - Destroy everything (the operator keeps the references it needs) */
2008: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
2009: {
2010: DM dmPressure;
2011: Vec constantPressure, basis;
2012: PetscReal nrm;
2013: MatNullSpace matNullSpace;
2015: PetscFunctionBeginUser;
2016: PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 0, 1, 0, &dmPressure));
2017: PetscCall(DMGetGlobalVector(dmPressure, &constantPressure));
2018: PetscCall(VecSet(constantPressure, 1.0));
2019: PetscCall(VecNorm(constantPressure, NORM_2, &nrm));
2020: PetscCall(VecScale(constantPressure, 1.0 / nrm));
2021: PetscCall(DMCreateGlobalVector(dmSol, &basis));
2022: PetscCall(DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis));
2023: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace));
2024: PetscCall(VecDestroy(&basis));
2025: PetscCall(MatSetNullSpace(A, matNullSpace));
2026: PetscCall(MatNullSpaceDestroy(&matNullSpace));
2027: PetscCall(DMRestoreGlobalVector(dmPressure, &constantPressure));
2028: PetscCall(DMDestroy(&dmPressure));
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2032: static PetscErrorCode DumpSolution(Ctx ctx, PetscInt level, Vec x)
2033: {
2034: DM dm_stokes, dm_coefficients;
2035: Vec coeff;
2036: DM dm_vel_avg;
2037: Vec vel_avg;
2038: DM da_vel_avg, da_p, da_eta_element;
2039: Vec vec_vel_avg, vec_p, vec_eta_element;
2040: 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;
2041: 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;
2042: PetscInt ex, ey, ez, startx, starty, startz, nx, ny, nz;
2043: Vec stokesLocal;
2045: PetscFunctionBeginUser;
2046: dm_stokes = ctx->levels[level]->dm_stokes;
2047: dm_coefficients = ctx->levels[level]->dm_coefficients;
2048: coeff = ctx->levels[level]->coeff;
2050: /* For convenience, create a new DM and Vec which will hold averaged velocities
2051: Note that this could also be accomplished with direct array access, using
2052: DMStagVecGetArray() and related functions */
2053: if (ctx->dim == 2) {
2054: PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 2, 0, &dm_vel_avg)); /* 2 dof per element */
2055: } else if (ctx->dim == 3) {
2056: PetscCall(DMStagCreateCompatibleDMStag(dm_stokes, 0, 0, 0, 3, &dm_vel_avg)); /* 3 dof per element */
2057: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not Implemented for dimension %" PetscInt_FMT, ctx->dim);
2058: PetscCall(DMSetUp(dm_vel_avg));
2059: PetscCall(DMStagSetUniformCoordinatesProduct(dm_vel_avg, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax, ctx->zmin, ctx->zmax));
2060: PetscCall(DMCreateGlobalVector(dm_vel_avg, &vel_avg));
2061: PetscCall(DMGetLocalVector(dm_stokes, &stokesLocal));
2062: PetscCall(DMGlobalToLocal(dm_stokes, x, INSERT_VALUES, stokesLocal));
2063: PetscCall(DMStagGetCorners(dm_vel_avg, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
2064: if (ctx->dim == 2) {
2065: for (ey = starty; ey < starty + ny; ++ey) {
2066: for (ex = startx; ex < startx + nx; ++ex) {
2067: DMStagStencil from[4], to[2];
2068: PetscScalar valFrom[4], valTo[2];
2070: from[0].i = ex;
2071: from[0].j = ey;
2072: from[0].loc = DMSTAG_UP;
2073: from[0].c = 0;
2074: from[1].i = ex;
2075: from[1].j = ey;
2076: from[1].loc = DMSTAG_DOWN;
2077: from[1].c = 0;
2078: from[2].i = ex;
2079: from[2].j = ey;
2080: from[2].loc = DMSTAG_LEFT;
2081: from[2].c = 0;
2082: from[3].i = ex;
2083: from[3].j = ey;
2084: from[3].loc = DMSTAG_RIGHT;
2085: from[3].c = 0;
2086: PetscCall(DMStagVecGetValuesStencil(dm_stokes, stokesLocal, 4, from, valFrom));
2087: to[0].i = ex;
2088: to[0].j = ey;
2089: to[0].loc = DMSTAG_ELEMENT;
2090: to[0].c = 0;
2091: valTo[0] = 0.5 * (valFrom[2] + valFrom[3]);
2092: to[1].i = ex;
2093: to[1].j = ey;
2094: to[1].loc = DMSTAG_ELEMENT;
2095: to[1].c = 1;
2096: valTo[1] = 0.5 * (valFrom[0] + valFrom[1]);
2097: PetscCall(DMStagVecSetValuesStencil(dm_vel_avg, vel_avg, 2, to, valTo, INSERT_VALUES));
2098: }
2099: }
2100: } else if (ctx->dim == 3) {
2101: for (ez = startz; ez < startz + nz; ++ez) {
2102: for (ey = starty; ey < starty + ny; ++ey) {
2103: for (ex = startx; ex < startx + nx; ++ex) {
2104: DMStagStencil from[6], to[3];
2105: PetscScalar valFrom[6], valTo[3];
2107: from[0].i = ex;
2108: from[0].j = ey;
2109: from[0].k = ez;
2110: from[0].loc = DMSTAG_UP;
2111: from[0].c = 0;
2112: from[1].i = ex;
2113: from[1].j = ey;
2114: from[1].k = ez;
2115: from[1].loc = DMSTAG_DOWN;
2116: from[1].c = 0;
2117: from[2].i = ex;
2118: from[2].j = ey;
2119: from[2].k = ez;
2120: from[2].loc = DMSTAG_LEFT;
2121: from[2].c = 0;
2122: from[3].i = ex;
2123: from[3].j = ey;
2124: from[3].k = ez;
2125: from[3].loc = DMSTAG_RIGHT;
2126: from[3].c = 0;
2127: from[4].i = ex;
2128: from[4].j = ey;
2129: from[4].k = ez;
2130: from[4].loc = DMSTAG_BACK;
2131: from[4].c = 0;
2132: from[5].i = ex;
2133: from[5].j = ey;
2134: from[5].k = ez;
2135: from[5].loc = DMSTAG_FRONT;
2136: from[5].c = 0;
2137: PetscCall(DMStagVecGetValuesStencil(dm_stokes, stokesLocal, 6, from, valFrom));
2138: to[0].i = ex;
2139: to[0].j = ey;
2140: to[0].k = ez;
2141: to[0].loc = DMSTAG_ELEMENT;
2142: to[0].c = 0;
2143: valTo[0] = 0.5 * (valFrom[2] + valFrom[3]);
2144: to[1].i = ex;
2145: to[1].j = ey;
2146: to[1].k = ez;
2147: to[1].loc = DMSTAG_ELEMENT;
2148: to[1].c = 1;
2149: valTo[1] = 0.5 * (valFrom[0] + valFrom[1]);
2150: to[2].i = ex;
2151: to[2].j = ey;
2152: to[2].k = ez;
2153: to[2].loc = DMSTAG_ELEMENT;
2154: to[2].c = 2;
2155: valTo[2] = 0.5 * (valFrom[4] + valFrom[5]);
2156: PetscCall(DMStagVecSetValuesStencil(dm_vel_avg, vel_avg, 3, to, valTo, INSERT_VALUES));
2157: }
2158: }
2159: }
2160: }
2161: PetscCall(VecAssemblyBegin(vel_avg));
2162: PetscCall(VecAssemblyEnd(vel_avg));
2163: PetscCall(DMRestoreLocalVector(dm_stokes, &stokesLocal));
2165: /* Create individual DMDAs for sub-grids of our DMStag objects. This is
2166: somewhat inefficient, but allows use of the DMDA API without re-implementing
2167: all utilities for DMStag */
2169: PetscCall(DMStagVecSplitToDMDA(dm_stokes, x, DMSTAG_ELEMENT, 0, &da_p, &vec_p));
2170: PetscCall(PetscObjectSetName((PetscObject)vec_p, "p (scaled)"));
2172: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_ELEMENT, 0, &da_eta_element, &vec_eta_element));
2173: PetscCall(PetscObjectSetName((PetscObject)vec_eta_element, "eta"));
2175: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_DOWN_LEFT, 0, &da_eta_down_left, &vec_eta_down_left));
2176: PetscCall(PetscObjectSetName((PetscObject)vec_eta_down_left, "eta"));
2178: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_DOWN_LEFT, 1, &da_rho_down_left, &vec_rho_down_left));
2179: PetscCall(PetscObjectSetName((PetscObject)vec_rho_down_left, "density"));
2181: if (ctx->dim == 3) {
2182: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_LEFT, 0, &da_eta_back_left, &vec_eta_back_left));
2183: PetscCall(PetscObjectSetName((PetscObject)vec_eta_back_left, "eta"));
2185: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_LEFT, 1, &da_rho_back_left, &vec_rho_back_left));
2186: PetscCall(PetscObjectSetName((PetscObject)vec_rho_back_left, "rho"));
2188: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_DOWN, 0, &da_eta_back_down, &vec_eta_back_down));
2189: PetscCall(PetscObjectSetName((PetscObject)vec_eta_back_down, "eta"));
2191: PetscCall(DMStagVecSplitToDMDA(dm_coefficients, coeff, DMSTAG_BACK_DOWN, 1, &da_rho_back_down, &vec_rho_back_down));
2192: PetscCall(PetscObjectSetName((PetscObject)vec_rho_back_down, "rho"));
2193: }
2195: PetscCall(DMStagVecSplitToDMDA(dm_vel_avg, vel_avg, DMSTAG_ELEMENT, -3, &da_vel_avg, &vec_vel_avg)); /* note -3 : pad with zero */
2196: PetscCall(PetscObjectSetName((PetscObject)vec_vel_avg, "Velocity (Averaged)"));
2198: /* Dump element-based fields to a .vtr file */
2199: {
2200: PetscViewer viewer;
2202: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_vel_avg), "ex4_element.vtr", FILE_MODE_WRITE, &viewer));
2203: PetscCall(VecView(vec_vel_avg, viewer));
2204: PetscCall(VecView(vec_p, viewer));
2205: PetscCall(VecView(vec_eta_element, viewer));
2206: PetscCall(PetscViewerDestroy(&viewer));
2207: }
2209: /* Dump vertex- or edge-based fields to a second .vtr file */
2210: {
2211: PetscViewer viewer;
2213: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_down_left), "ex4_down_left.vtr", FILE_MODE_WRITE, &viewer));
2214: PetscCall(VecView(vec_eta_down_left, viewer));
2215: PetscCall(VecView(vec_rho_down_left, viewer));
2216: PetscCall(PetscViewerDestroy(&viewer));
2217: }
2218: if (ctx->dim == 3) {
2219: PetscViewer viewer;
2221: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_back_left), "ex4_back_left.vtr", FILE_MODE_WRITE, &viewer));
2222: PetscCall(VecView(vec_eta_back_left, viewer));
2223: PetscCall(VecView(vec_rho_back_left, viewer));
2224: PetscCall(PetscViewerDestroy(&viewer));
2225: }
2226: if (ctx->dim == 3) {
2227: PetscViewer viewer;
2229: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)da_eta_back_down), "ex4_back_down.vtr", FILE_MODE_WRITE, &viewer));
2230: PetscCall(VecView(vec_eta_back_down, viewer));
2231: PetscCall(VecView(vec_rho_back_down, viewer));
2232: PetscCall(PetscViewerDestroy(&viewer));
2233: }
2235: /* Destroy DMDAs and Vecs */
2236: PetscCall(VecDestroy(&vec_vel_avg));
2237: PetscCall(VecDestroy(&vec_p));
2238: PetscCall(VecDestroy(&vec_eta_element));
2239: PetscCall(VecDestroy(&vec_eta_down_left));
2240: if (ctx->dim == 3) {
2241: PetscCall(VecDestroy(&vec_eta_back_left));
2242: PetscCall(VecDestroy(&vec_eta_back_down));
2243: }
2244: PetscCall(VecDestroy(&vec_rho_down_left));
2245: if (ctx->dim == 3) {
2246: PetscCall(VecDestroy(&vec_rho_back_left));
2247: PetscCall(VecDestroy(&vec_rho_back_down));
2248: }
2249: PetscCall(DMDestroy(&da_vel_avg));
2250: PetscCall(DMDestroy(&da_p));
2251: PetscCall(DMDestroy(&da_eta_element));
2252: PetscCall(DMDestroy(&da_eta_down_left));
2253: if (ctx->dim == 3) {
2254: PetscCall(DMDestroy(&da_eta_back_left));
2255: PetscCall(DMDestroy(&da_eta_back_down));
2256: }
2257: PetscCall(DMDestroy(&da_rho_down_left));
2258: if (ctx->dim == 3) {
2259: PetscCall(DMDestroy(&da_rho_back_left));
2260: PetscCall(DMDestroy(&da_rho_back_down));
2261: }
2262: PetscCall(VecDestroy(&vel_avg));
2263: PetscCall(DMDestroy(&dm_vel_avg));
2264: PetscFunctionReturn(PETSC_SUCCESS);
2265: }
2267: /*TEST
2269: test:
2270: suffix: direct_umfpack
2271: requires: suitesparse !complex
2272: nsize: 1
2273: 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
2275: test:
2276: suffix: direct_mumps
2277: requires: mumps !complex !single
2278: nsize: 9
2279: 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
2281: test:
2282: suffix: isovisc_nondim_abf_mg
2283: nsize: 1
2284: 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
2286: test:
2287: suffix: isovisc_nondim_abf_mg_2
2288: nsize: 1
2289: 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
2291: test:
2292: suffix: nondim_abf_mg
2293: requires: suitesparse !complex
2294: nsize: 4
2295: 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
2297: test:
2298: suffix: nondim_abf_lu
2299: requires: suitesparse !complex
2300: nsize: 1
2301: 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
2303: test:
2304: suffix: 3d_nondim_isovisc_abf_mg
2305: requires: !single
2306: nsize: 1
2307: 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
2309: test:
2310: TODO: unstable across systems
2311: suffix: monolithic
2312: nsize: 1
2313: requires: double !complex
2314: 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
2316: test:
2317: suffix: 3d_nondim_isovisc_sinker_abf_mg
2318: requires: !complex !single
2319: nsize: 1
2320: 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
2322: test:
2323: TODO: unstable across systems
2324: suffix: 3d_nondim_mono_mg_lamemstyle
2325: nsize: 1
2326: requires: suitesparse
2327: 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
2329: test:
2330: suffix: 3d_isovisc_blob_cuda
2331: requires: cuda defined(PETSC_USE_LOG)
2332: nsize: 2
2333: 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
2334: filter: awk "/MGInterp/ {print \$NF}"
2336: TEST*/