Actual source code: ex62.c
1: static char help[] = "Stokes Problem discretized with finite elements,\n\
2: using a parallel unstructured mesh (DMPLEX) to represent the domain.\n\n\n";
4: /*
5: For the isoviscous Stokes problem, which we discretize using the finite
6: element method on an unstructured mesh, the weak form equations are
8: < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > - < v, f > = 0
9: < q, -\nabla\cdot u > = 0
11: Viewing:
13: To produce nice output, use
15: -dm_refine 3 -dm_view hdf5:sol1.h5 -error_vec_view hdf5:sol1.h5::append -snes_view_solution hdf5:sol1.h5::append -exact_vec_view hdf5:sol1.h5::append
17: You can get a LaTeX view of the mesh, with point numbering using
19: -dm_view :mesh.tex:ascii_latex -dm_plex_view_scale 8.0
21: The data layout can be viewed using
23: -dm_petscsection_view
25: Lots of information about the FEM assembly can be printed using
27: -dm_plex_print_fem 3
28: */
30: #include <petscdmplex.h>
31: #include <petscsnes.h>
32: #include <petscds.h>
33: #include <petscbag.h>
35: // TODO: Plot residual by fields after each smoother iterate
37: typedef enum {
38: SOL_QUADRATIC,
39: SOL_TRIG,
40: SOL_UNKNOWN
41: } SolType;
42: const char *SolTypes[] = {"quadratic", "trig", "unknown", "SolType", "SOL_", 0};
44: typedef struct {
45: PetscScalar mu; /* dynamic shear viscosity */
46: } Parameter;
48: typedef struct {
49: PetscBag bag; /* Problem parameters */
50: SolType sol; /* MMS solution */
51: } AppCtx;
53: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
54: {
55: const PetscReal mu = PetscRealPart(constants[0]);
56: const PetscInt Nc = uOff[1] - uOff[0];
57: PetscInt c, d;
59: for (c = 0; c < Nc; ++c) {
60: for (d = 0; d < dim; ++d) f1[c * dim + d] = mu * (u_x[c * dim + d] + u_x[d * dim + c]);
61: f1[c * dim + c] -= u[uOff[1]];
62: }
63: }
65: static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
66: {
67: PetscInt d;
68: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] -= u_x[d * dim + d];
69: }
71: static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
72: {
73: PetscInt d;
74: for (d = 0; d < dim; ++d) g1[d * dim + d] = -1.0; /* < q, -\nabla\cdot u > */
75: }
77: static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
78: {
79: PetscInt d;
80: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* -< \nabla\cdot v, p > */
81: }
83: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
84: {
85: const PetscReal mu = PetscRealPart(constants[0]);
86: const PetscInt Nc = uOff[1] - uOff[0];
87: PetscInt c, d;
89: for (c = 0; c < Nc; ++c) {
90: for (d = 0; d < dim; ++d) {
91: g3[((c * Nc + c) * dim + d) * dim + d] += mu; /* < \nabla v, \nabla u > */
92: g3[((c * Nc + d) * dim + d) * dim + c] += mu; /* < \nabla v, {\nabla u}^T > */
93: }
94: }
95: }
97: static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
98: {
99: const PetscReal mu = PetscRealPart(constants[0]);
101: g0[0] = 1.0 / mu;
102: }
104: /* Quadratic MMS Solution
105: 2D:
107: u = x^2 + y^2
108: v = 2 x^2 - 2xy
109: p = x + y - 1
110: f = <1 - 4 mu, 1 - 4 mu>
112: so that
114: e(u) = (grad u + grad u^T) = / 4x 4x \
115: \ 4x -4x /
116: div mu e(u) - \nabla p + f = mu <4, 4> - <1, 1> + <1 - 4 mu, 1 - 4 mu> = 0
117: \nabla \cdot u = 2x - 2x = 0
119: 3D:
121: u = 2 x^2 + y^2 + z^2
122: v = 2 x^2 - 2xy
123: w = 2 x^2 - 2xz
124: p = x + y + z - 3/2
125: f = <1 - 8 mu, 1 - 4 mu, 1 - 4 mu>
127: so that
129: e(u) = (grad u + grad u^T) = / 8x 4x 4x \
130: | 4x -4x 0 |
131: \ 4x 0 -4x /
132: div mu e(u) - \nabla p + f = mu <8, 4, 4> - <1, 1, 1> + <1 - 8 mu, 1 - 4 mu, 1 - 4 mu> = 0
133: \nabla \cdot u = 4x - 2x - 2x = 0
134: */
135: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
136: {
137: PetscInt c;
139: u[0] = (dim - 1) * PetscSqr(x[0]);
140: for (c = 1; c < Nc; ++c) {
141: u[0] += PetscSqr(x[c]);
142: u[c] = 2.0 * PetscSqr(x[0]) - 2.0 * x[0] * x[c];
143: }
144: return 0;
145: }
147: static PetscErrorCode quadratic_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
148: {
149: PetscInt d;
151: u[0] = -0.5 * dim;
152: for (d = 0; d < dim; ++d) u[0] += x[d];
153: return 0;
154: }
156: static void f0_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
157: {
158: const PetscReal mu = PetscRealPart(constants[0]);
159: PetscInt d;
161: f0[0] = (dim - 1) * 4.0 * mu - 1.0;
162: for (d = 1; d < dim; ++d) f0[d] = 4.0 * mu - 1.0;
163: }
165: /* Trigonometric MMS Solution
166: 2D:
168: u = sin(pi x) + sin(pi y)
169: v = -pi cos(pi x) y
170: p = sin(2 pi x) + sin(2 pi y)
171: f = <2pi cos(2 pi x) + mu pi^2 sin(pi x) + mu pi^2 sin(pi y), 2pi cos(2 pi y) - mu pi^3 cos(pi x) y>
173: so that
175: e(u) = (grad u + grad u^T) = / 2pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y \
176: \ pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) /
177: div mu e(u) - \nabla p + f = mu <-pi^2 sin(pi x) - pi^2 sin(pi y), pi^3 cos(pi x) y> - <2pi cos(2 pi x), 2pi cos(2 pi y)> + <f_x, f_y> = 0
178: \nabla \cdot u = pi cos(pi x) - pi cos(pi x) = 0
180: 3D:
182: u = 2 sin(pi x) + sin(pi y) + sin(pi z)
183: v = -pi cos(pi x) y
184: w = -pi cos(pi x) z
185: p = sin(2 pi x) + sin(2 pi y) + sin(2 pi z)
186: f = <2pi cos(2 pi x) + mu 2pi^2 sin(pi x) + mu pi^2 sin(pi y) + mu pi^2 sin(pi z), 2pi cos(2 pi y) - mu pi^3 cos(pi x) y, 2pi cos(2 pi z) - mu pi^3 cos(pi x) z>
188: so that
190: e(u) = (grad u + grad u^T) = / 4pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y pi cos(pi z) + pi^2 sin(pi x) z \
191: | pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) 0 |
192: \ pi cos(pi z) + pi^2 sin(pi x) z 0 -2pi cos(pi x) /
193: div mu e(u) - \nabla p + f = mu <-2pi^2 sin(pi x) - pi^2 sin(pi y) - pi^2 sin(pi z), pi^3 cos(pi x) y, pi^3 cos(pi x) z> - <2pi cos(2 pi x), 2pi cos(2 pi y), 2pi cos(2 pi z)> + <f_x, f_y, f_z> = 0
194: \nabla \cdot u = 2 pi cos(pi x) - pi cos(pi x) - pi cos(pi x) = 0
195: */
196: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
197: {
198: PetscInt c;
200: u[0] = (dim - 1) * PetscSinReal(PETSC_PI * x[0]);
201: for (c = 1; c < Nc; ++c) {
202: u[0] += PetscSinReal(PETSC_PI * x[c]);
203: u[c] = -PETSC_PI * PetscCosReal(PETSC_PI * x[0]) * x[c];
204: }
205: return 0;
206: }
208: static PetscErrorCode trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
209: {
210: PetscInt d;
212: for (d = 0, u[0] = 0.0; d < dim; ++d) u[0] += PetscSinReal(2.0 * PETSC_PI * x[d]);
213: return 0;
214: }
216: static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
217: {
218: const PetscReal mu = PetscRealPart(constants[0]);
219: PetscInt d;
221: f0[0] = -2.0 * PETSC_PI * PetscCosReal(2.0 * PETSC_PI * x[0]) - (dim - 1) * mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[0]);
222: for (d = 1; d < dim; ++d) {
223: f0[0] -= mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[d]);
224: f0[d] = -2.0 * PETSC_PI * PetscCosReal(2.0 * PETSC_PI * x[d]) + mu * PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * x[0]) * x[d];
225: }
226: }
228: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
229: {
230: PetscInt sol;
233: options->sol = SOL_QUADRATIC;
234: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
235: sol = options->sol;
236: PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, PETSC_STATIC_ARRAY_LENGTH(SolTypes) - 3, SolTypes[options->sol], &sol, NULL);
237: options->sol = (SolType)sol;
238: PetscOptionsEnd();
239: return 0;
240: }
242: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
243: {
245: DMCreate(comm, dm);
246: DMSetType(*dm, DMPLEX);
247: DMSetFromOptions(*dm);
248: DMViewFromOptions(*dm, NULL, "-dm_view");
249: return 0;
250: }
252: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
253: {
254: Parameter *p;
257: /* setup PETSc parameter bag */
258: PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx->bag);
259: PetscBagGetData(ctx->bag, (void **)&p);
260: PetscBagSetName(ctx->bag, "par", "Stokes Parameters");
261: PetscBagRegisterScalar(ctx->bag, &p->mu, 1.0, "mu", "Dynamic Shear Viscosity, Pa s");
262: PetscBagSetFromOptions(ctx->bag);
263: {
264: PetscViewer viewer;
265: PetscViewerFormat format;
266: PetscBool flg;
268: PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg);
269: if (flg) {
270: PetscViewerPushFormat(viewer, format);
271: PetscBagView(ctx->bag, viewer);
272: PetscViewerFlush(viewer);
273: PetscViewerPopFormat(viewer);
274: PetscViewerDestroy(&viewer);
275: }
276: }
277: return 0;
278: }
280: static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
281: {
282: PetscErrorCode (*exactFuncs[2])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
283: PetscDS ds;
284: DMLabel label;
285: const PetscInt id = 1;
288: DMGetDS(dm, &ds);
289: switch (user->sol) {
290: case SOL_QUADRATIC:
291: PetscDSSetResidual(ds, 0, f0_quadratic_u, f1_u);
292: exactFuncs[0] = quadratic_u;
293: exactFuncs[1] = quadratic_p;
294: break;
295: case SOL_TRIG:
296: PetscDSSetResidual(ds, 0, f0_trig_u, f1_u);
297: exactFuncs[0] = trig_u;
298: exactFuncs[1] = trig_p;
299: break;
300: default:
301: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", SolTypes[PetscMin(user->sol, SOL_UNKNOWN)], user->sol);
302: }
303: PetscDSSetResidual(ds, 1, f0_p, NULL);
304: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
305: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
306: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
307: PetscDSSetJacobianPreconditioner(ds, 0, 0, NULL, NULL, NULL, g3_uu);
308: PetscDSSetJacobianPreconditioner(ds, 1, 1, g0_pp, NULL, NULL, NULL);
310: PetscDSSetExactSolution(ds, 0, exactFuncs[0], user);
311: PetscDSSetExactSolution(ds, 1, exactFuncs[1], user);
313: DMGetLabel(dm, "marker", &label);
314: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, user, NULL);
316: /* Make constant values available to pointwise functions */
317: {
318: Parameter *param;
319: PetscScalar constants[1];
321: PetscBagGetData(user->bag, (void **)¶m);
322: constants[0] = param->mu; /* dynamic shear viscosity, Pa s */
323: PetscDSSetConstants(ds, 1, constants);
324: }
325: return 0;
326: }
328: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
329: {
330: PetscInt c;
331: for (c = 0; c < Nc; ++c) u[c] = 0.0;
332: return 0;
333: }
334: static PetscErrorCode one(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
335: {
336: PetscInt c;
337: for (c = 0; c < Nc; ++c) u[c] = 1.0;
338: return 0;
339: }
341: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
342: {
343: Vec vec;
344: PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {zero, one};
348: funcs[field] = one;
349: {
350: PetscDS ds;
351: DMGetDS(dm, &ds);
352: PetscObjectViewFromOptions((PetscObject)ds, NULL, "-ds_view");
353: }
354: DMCreateGlobalVector(dm, &vec);
355: DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);
356: VecNormalize(vec, NULL);
357: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullspace);
358: VecDestroy(&vec);
359: /* New style for field null spaces */
360: {
361: PetscObject pressure;
362: MatNullSpace nullspacePres;
364: DMGetField(dm, field, NULL, &pressure);
365: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres);
366: PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres);
367: MatNullSpaceDestroy(&nullspacePres);
368: }
369: return 0;
370: }
372: static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
373: {
374: DM cdm = dm;
375: PetscQuadrature q = NULL;
376: PetscBool simplex;
377: PetscInt dim, Nf = 2, f, Nc[2];
378: const char *name[2] = {"velocity", "pressure"};
379: const char *prefix[2] = {"vel_", "pres_"};
381: DMGetDimension(dm, &dim);
382: DMPlexIsSimplex(dm, &simplex);
383: Nc[0] = dim;
384: Nc[1] = 1;
385: for (f = 0; f < Nf; ++f) {
386: PetscFE fe;
388: PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe);
389: PetscObjectSetName((PetscObject)fe, name[f]);
390: if (!q) PetscFEGetQuadrature(fe, &q);
391: PetscFESetQuadrature(fe, q);
392: DMSetField(dm, f, NULL, (PetscObject)fe);
393: PetscFEDestroy(&fe);
394: }
395: DMCreateDS(dm);
396: (*setupEqn)(dm, user);
397: while (cdm) {
398: DMCopyDisc(dm, cdm);
399: DMSetNullSpaceConstructor(cdm, 1, CreatePressureNullSpace);
400: DMGetCoarseDM(cdm, &cdm);
401: }
402: return 0;
403: }
405: int main(int argc, char **argv)
406: {
407: SNES snes;
408: DM dm;
409: Vec u;
410: AppCtx user;
413: PetscInitialize(&argc, &argv, NULL, help);
414: ProcessOptions(PETSC_COMM_WORLD, &user);
415: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
416: SNESCreate(PetscObjectComm((PetscObject)dm), &snes);
417: SNESSetDM(snes, dm);
418: DMSetApplicationContext(dm, &user);
420: SetupParameters(PETSC_COMM_WORLD, &user);
421: SetupProblem(dm, SetupEqn, &user);
422: DMPlexCreateClosureIndex(dm, NULL);
424: DMCreateGlobalVector(dm, &u);
425: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
426: SNESSetFromOptions(snes);
427: DMSNESCheckFromOptions(snes, u);
428: PetscObjectSetName((PetscObject)u, "Solution");
429: {
430: Mat J;
431: MatNullSpace sp;
433: SNESSetUp(snes);
434: CreatePressureNullSpace(dm, 1, 1, &sp);
435: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
436: MatSetNullSpace(J, sp);
437: MatNullSpaceDestroy(&sp);
438: PetscObjectSetName((PetscObject)J, "Jacobian");
439: MatViewFromOptions(J, NULL, "-J_view");
440: }
441: SNESSolve(snes, NULL, u);
443: VecDestroy(&u);
444: SNESDestroy(&snes);
445: DMDestroy(&dm);
446: PetscBagDestroy(&user.bag);
447: PetscFinalize();
448: return 0;
449: }
450: /*TEST
452: test:
453: suffix: 2d_p2_p1_check
454: requires: triangle
455: args: -sol quadratic -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
457: test:
458: suffix: 2d_p2_p1_check_parallel
459: nsize: {{2 3 5}}
460: requires: triangle
461: args: -sol quadratic -dm_refine 2 -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
463: test:
464: suffix: 3d_p2_p1_check
465: requires: ctetgen
466: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
468: test:
469: suffix: 3d_p2_p1_check_parallel
470: nsize: {{2 3 5}}
471: requires: ctetgen
472: args: -sol quadratic -dm_refine 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
474: test:
475: suffix: 2d_p2_p1_conv
476: requires: triangle
477: # Using -dm_refine 3 gives L_2 convergence rate: [3.0, 2.1]
478: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
479: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
480: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
481: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
483: test:
484: suffix: 2d_p2_p1_conv_gamg
485: requires: triangle
486: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 \
487: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
488: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_coarse_pc_type svd
490: test:
491: suffix: 3d_p2_p1_conv
492: requires: ctetgen !single
493: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.8]
494: args: -sol trig -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
495: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
496: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
497: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
499: test:
500: suffix: 2d_q2_q1_check
501: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
503: test:
504: suffix: 3d_q2_q1_check
505: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
507: test:
508: suffix: 2d_q2_q1_conv
509: # Using -dm_refine 3 -convest_num_refine 1 gives L_2 convergence rate: [3.0, 2.1]
510: args: -sol trig -dm_plex_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -ksp_error_if_not_converged \
511: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
512: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
513: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
515: test:
516: suffix: 3d_q2_q1_conv
517: requires: !single
518: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.4]
519: args: -sol trig -dm_plex_simplex 0 -dm_plex_dim 3 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
520: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
521: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
522: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
524: test:
525: suffix: 2d_p3_p2_check
526: requires: triangle
527: args: -sol quadratic -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
529: test:
530: suffix: 3d_p3_p2_check
531: requires: ctetgen !single
532: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
534: test:
535: suffix: 2d_p3_p2_conv
536: requires: triangle
537: # Using -dm_refine 2 gives L_2 convergence rate: [3.8, 3.0]
538: args: -sol trig -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
539: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
540: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
541: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
543: test:
544: suffix: 3d_p3_p2_conv
545: requires: ctetgen long_runtime
546: # Using -dm_refine 1 -convest_num_refine 2 gives L_2 convergence rate: [3.6, 3.9]
547: args: -sol trig -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 \
548: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
549: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
550: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
552: test:
553: suffix: 2d_q1_p0_conv
554: requires: !single
555: # Using -dm_refine 3 gives L_2 convergence rate: [1.9, 1.0]
556: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 2 \
557: -ksp_atol 1e-10 -petscds_jac_pre 0 \
558: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
559: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_levels_pc_type jacobi -fieldsplit_pressure_mg_coarse_pc_type svd -fieldsplit_pressure_pc_gamg_aggressive_coarsening 0
561: test:
562: suffix: 3d_q1_p0_conv
563: requires: !single
564: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [1.7, 1.0]
565: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 1 \
566: -ksp_atol 1e-10 -petscds_jac_pre 0 \
567: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
568: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_levels_pc_type jacobi -fieldsplit_pressure_mg_coarse_pc_type svd -fieldsplit_pressure_pc_gamg_aggressive_coarsening 0
570: # Stokes preconditioners
571: # Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
572: test:
573: suffix: 2d_p2_p1_block_diagonal
574: requires: triangle
575: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
576: -snes_error_if_not_converged \
577: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-4 -ksp_error_if_not_converged \
578: -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
579: # Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
580: test:
581: suffix: 2d_p2_p1_block_triangular
582: requires: triangle
583: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
584: -snes_error_if_not_converged \
585: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
586: -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
587: # Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
588: test:
589: suffix: 2d_p2_p1_schur_diagonal
590: requires: triangle
591: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
592: -snes_error_if_not_converged \
593: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
594: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -pc_fieldsplit_off_diag_use_amat \
595: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
596: # Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
597: test:
598: suffix: 2d_p2_p1_schur_upper
599: requires: triangle
600: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001 \
601: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
602: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -pc_fieldsplit_off_diag_use_amat \
603: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
604: # Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
605: test:
606: suffix: 2d_p2_p1_schur_lower
607: requires: triangle
608: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
609: -snes_error_if_not_converged \
610: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
611: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -pc_fieldsplit_off_diag_use_amat \
612: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
613: # Full Schur complement \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix} \begin{pmatrix} I & A^{-1} B \\ 0 & I \end{pmatrix}
614: test:
615: suffix: 2d_p2_p1_schur_full
616: requires: triangle
617: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
618: -snes_error_if_not_converged \
619: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
620: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_off_diag_use_amat \
621: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
622: # Full Schur + Velocity GMG
623: test:
624: suffix: 2d_p2_p1_gmg_vcycle
625: requires: triangle
626: args: -sol quadratic -dm_refine_hierarchy 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
627: -ksp_type fgmres -ksp_atol 1e-9 -snes_error_if_not_converged -pc_use_amat \
628: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_off_diag_use_amat \
629: -fieldsplit_velocity_pc_type mg -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_pc_gamg_esteig_ksp_max_it 10 -fieldsplit_pressure_mg_levels_pc_type sor -fieldsplit_pressure_mg_coarse_pc_type svd
630: # SIMPLE \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & B^T diag(A)^{-1} B \end{pmatrix} \begin{pmatrix} I & diag(A)^{-1} B \\ 0 & I \end{pmatrix}
631: test:
632: suffix: 2d_p2_p1_simple
633: requires: triangle
634: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
635: -snes_error_if_not_converged \
636: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
637: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
638: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
639: -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi
640: # FETI-DP solvers (these solvers are quite inefficient, they are here to exercise the code)
641: test:
642: suffix: 2d_p2_p1_fetidp
643: requires: triangle mumps
644: nsize: 5
645: args: -sol quadratic -dm_refine 2 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
646: -snes_error_if_not_converged \
647: -ksp_type fetidp -ksp_rtol 1.0e-8 \
648: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
649: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 200 -fetidp_fieldsplit_p_pc_type none \
650: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type mumps -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type mumps -fetidp_fieldsplit_lag_ksp_type preonly
651: test:
652: suffix: 2d_q2_q1_fetidp
653: requires: mumps
654: nsize: 5
655: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
656: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
657: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
658: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 200 -fetidp_fieldsplit_p_pc_type none \
659: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type mumps -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type mumps -fetidp_fieldsplit_lag_ksp_type preonly
660: test:
661: suffix: 3d_p2_p1_fetidp
662: requires: ctetgen mumps suitesparse
663: nsize: 5
664: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
665: -snes_error_if_not_converged \
666: -ksp_type fetidp -ksp_rtol 1.0e-9 \
667: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
668: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 1000 -fetidp_fieldsplit_p_pc_type none \
669: -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_benign_trick -fetidp_bddc_pc_bddc_deluxe_singlemat \
670: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
671: -fetidp_bddelta_pc_factor_mat_solver_type umfpack -fetidp_fieldsplit_lag_ksp_type preonly -fetidp_bddc_sub_schurs_mat_solver_type mumps -fetidp_bddc_sub_schurs_mat_mumps_icntl_14 100000 \
672: -fetidp_bddelta_pc_factor_mat_ordering_type external \
673: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
674: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
675: test:
676: suffix: 3d_q2_q1_fetidp
677: requires: suitesparse
678: nsize: 5
679: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
680: -snes_error_if_not_converged \
681: -ksp_type fetidp -ksp_rtol 1.0e-8 \
682: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
683: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 2000 -fetidp_fieldsplit_p_pc_type none \
684: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
685: -fetidp_bddc_pc_bddc_symmetric -fetidp_fieldsplit_lag_ksp_type preonly \
686: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
687: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
688: # BDDC solvers (these solvers are quite inefficient, they are here to exercise the code)
689: test:
690: suffix: 2d_p2_p1_bddc
691: nsize: 2
692: requires: triangle !single
693: args: -sol quadratic -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
694: -snes_error_if_not_converged \
695: -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
696: -pc_type bddc -pc_bddc_corner_selection -pc_bddc_dirichlet_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_coarse_redundant_pc_type svd
697: # Vanka
698: test:
699: suffix: 2d_q1_p0_vanka
700: requires: double !complex
701: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
702: -snes_rtol 1.0e-4 \
703: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
704: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
705: -sub_ksp_type preonly -sub_pc_type lu
706: test:
707: suffix: 2d_q1_p0_vanka_denseinv
708: requires: double !complex
709: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
710: -snes_rtol 1.0e-4 \
711: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
712: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
713: -pc_patch_dense_inverse -pc_patch_sub_mat_type seqdense
714: # Vanka smoother
715: test:
716: suffix: 2d_q1_p0_gmg_vanka
717: requires: double !complex
718: args: -sol quadratic -dm_plex_simplex 0 -dm_refine_hierarchy 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
719: -snes_rtol 1.0e-4 \
720: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
721: -pc_type mg \
722: -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 30 \
723: -mg_levels_pc_type patch -mg_levels_pc_patch_partition_of_unity 0 -mg_levels_pc_patch_construct_codim 0 -mg_levels_pc_patch_construct_type vanka \
724: -mg_levels_sub_ksp_type preonly -mg_levels_sub_pc_type lu \
725: -mg_coarse_pc_type svd
727: TEST*/