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 {SOL_QUADRATIC, SOL_TRIG, SOL_UNKNOWN} SolType;
38: const char *SolTypes[] = {"quadratic", "trig", "unknown", "SolType", "SOL_", 0};
40: typedef struct {
41: PetscScalar mu; /* dynamic shear viscosity */
42: } Parameter;
44: typedef struct {
45: /* Domain and mesh definition */
46: char filename[PETSC_MAX_PATH_LEN];
47: /* Problem definition */
48: PetscBag bag; /* Problem parameters */
49: SolType sol; /* MMS solution */
50: } AppCtx;
52: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
53: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
54: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
55: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
56: {
57: const PetscReal mu = PetscRealPart(constants[0]);
58: const PetscInt Nc = uOff[1]-uOff[0];
59: PetscInt c, d;
61: for (c = 0; c < Nc; ++c) {
62: for (d = 0; d < dim; ++d) {
63: f1[c*dim+d] = mu * (u_x[c*dim+d] + u_x[d*dim+c]);
64: }
65: f1[c*dim+c] -= u[uOff[1]];
66: }
67: }
69: static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
70: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
71: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
72: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
73: {
74: PetscInt d;
75: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] -= u_x[d*dim+d];
76: }
78: static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
79: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
80: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
81: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
82: {
83: PetscInt d;
84: for (d = 0; d < dim; ++d) g1[d*dim+d] = -1.0; /* < q, -\nabla\cdot u > */
85: }
87: static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux,
88: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
89: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
90: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
91: {
92: PetscInt d;
93: for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* -< \nabla\cdot v, p > */
94: }
96: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
97: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
98: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
99: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
100: {
101: const PetscReal mu = PetscRealPart(constants[0]);
102: const PetscInt Nc = uOff[1]-uOff[0];
103: PetscInt c, d;
105: for (c = 0; c < Nc; ++c) {
106: for (d = 0; d < dim; ++d) {
107: g3[((c*Nc+c)*dim+d)*dim+d] += mu; /* < \nabla v, \nabla u > */
108: g3[((c*Nc+d)*dim+d)*dim+c] += mu; /* < \nabla v, {\nabla u}^T > */
109: }
110: }
111: }
113: static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
114: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
115: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
116: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
117: {
118: const PetscReal mu = PetscRealPart(constants[0]);
120: g0[0] = 1.0/mu;
121: }
123: /* Quadratic MMS Solution
124: 2D:
126: u = x^2 + y^2
127: v = 2 x^2 - 2xy
128: p = x + y - 1
129: f = <1 - 4 mu, 1 - 4 mu>
131: so that
133: e(u) = (grad u + grad u^T) = / 4x 4x \
134: \ 4x -4x /
135: div mu e(u) - \nabla p + f = mu <4, 4> - <1, 1> + <1 - 4 mu, 1 - 4 mu> = 0
136: \nabla \cdot u = 2x - 2x = 0
138: 3D:
140: u = 2 x^2 + y^2 + z^2
141: v = 2 x^2 - 2xy
142: w = 2 x^2 - 2xz
143: p = x + y + z - 3/2
144: f = <1 - 8 mu, 1 - 4 mu, 1 - 4 mu>
146: so that
148: e(u) = (grad u + grad u^T) = / 8x 4x 4x \
149: | 4x -4x 0 |
150: \ 4x 0 -4x /
151: div mu e(u) - \nabla p + f = mu <8, 4, 4> - <1, 1, 1> + <1 - 8 mu, 1 - 4 mu, 1 - 4 mu> = 0
152: \nabla \cdot u = 4x - 2x - 2x = 0
153: */
154: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
155: {
156: PetscInt c;
158: u[0] = (dim-1)*PetscSqr(x[0]);
159: for (c = 1; c < Nc; ++c) {
160: u[0] += PetscSqr(x[c]);
161: u[c] = 2.0*PetscSqr(x[0]) - 2.0*x[0]*x[c];
162: }
163: return 0;
164: }
166: static PetscErrorCode quadratic_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
167: {
168: PetscInt d;
170: u[0] = -0.5*dim;
171: for (d = 0; d < dim; ++d) u[0] += x[d];
172: return 0;
173: }
175: static void f0_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
176: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
177: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
178: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
179: {
180: const PetscReal mu = PetscRealPart(constants[0]);
181: PetscInt d;
183: f0[0] = (dim-1)*4.0*mu - 1.0;
184: for (d = 1; d < dim; ++d) f0[d] = 4.0*mu - 1.0;
185: }
187: /* Trigonometric MMS Solution
188: 2D:
190: u = sin(pi x) + sin(pi y)
191: v = -pi cos(pi x) y
192: p = sin(2 pi x) + sin(2 pi y)
193: 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>
195: so that
197: e(u) = (grad u + grad u^T) = / 2pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y \
198: \ pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) /
199: 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
200: \nabla \cdot u = pi cos(pi x) - pi cos(pi x) = 0
202: 3D:
204: u = 2 sin(pi x) + sin(pi y) + sin(pi z)
205: v = -pi cos(pi x) y
206: w = -pi cos(pi x) z
207: p = sin(2 pi x) + sin(2 pi y) + sin(2 pi z)
208: 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>
210: so that
212: 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 \
213: | pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) 0 |
214: \ pi cos(pi z) + pi^2 sin(pi x) z 0 -2pi cos(pi x) /
215: 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
216: \nabla \cdot u = 2 pi cos(pi x) - pi cos(pi x) - pi cos(pi x) = 0
217: */
218: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
219: {
220: PetscInt c;
222: u[0] = (dim-1)*PetscSinReal(PETSC_PI*x[0]);
223: for (c = 1; c < Nc; ++c) {
224: u[0] += PetscSinReal(PETSC_PI*x[c]);
225: u[c] = -PETSC_PI*PetscCosReal(PETSC_PI*x[0]) * x[c];
226: }
227: return 0;
228: }
230: static PetscErrorCode trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
231: {
232: PetscInt d;
234: for (d = 0, u[0] = 0.0; d < dim; ++d) u[0] += PetscSinReal(2.0*PETSC_PI*x[d]);
235: return 0;
236: }
238: static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
239: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
240: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
241: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
242: {
243: const PetscReal mu = PetscRealPart(constants[0]);
244: PetscInt d;
246: f0[0] = -2.0*PETSC_PI*PetscCosReal(2.0*PETSC_PI*x[0]) - (dim-1)*mu*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x[0]);
247: for (d = 1; d < dim; ++d) {
248: f0[0] -= mu*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x[d]);
249: 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];
250: }
251: }
253: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
254: {
255: PetscInt sol;
259: options->filename[0] = '\0';
260: options->sol = SOL_QUADRATIC;
262: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
263: sol = options->sol;
264: PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, (sizeof(SolTypes)/sizeof(SolTypes[0]))-3, SolTypes[options->sol], &sol, NULL);
265: options->sol = (SolType) sol;
266: PetscOptionsEnd();
267: return(0);
268: }
270: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
271: {
272: size_t len;
276: PetscStrlen(user->filename, &len);
277: if (len) {DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);}
278: else {DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);}
279: DMSetFromOptions(*dm);
280: DMViewFromOptions(*dm, NULL, "-dm_view");
281: return(0);
282: }
284: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
285: {
286: Parameter *p;
290: /* setup PETSc parameter bag */
291: PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx->bag);
292: PetscBagGetData(ctx->bag, (void **) &p);
293: PetscBagSetName(ctx->bag, "par", "Stokes Parameters");
294: PetscBagRegisterScalar(ctx->bag, &p->mu, 1.0, "mu", "Dynamic Shear Viscosity, Pa s");
295: PetscBagSetFromOptions(ctx->bag);
296: {
297: PetscViewer viewer;
298: PetscViewerFormat format;
299: PetscBool flg;
301: PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg);
302: if (flg) {
303: PetscViewerPushFormat(viewer, format);
304: PetscBagView(ctx->bag, viewer);
305: PetscViewerFlush(viewer);
306: PetscViewerPopFormat(viewer);
307: PetscViewerDestroy(&viewer);
308: }
309: }
310: return(0);
311: }
313: static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
314: {
315: PetscErrorCode (*exactFuncs[2])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
316: PetscDS ds;
317: const PetscInt id = 1;
318: PetscErrorCode ierr;
321: DMGetDS(dm, &ds);
322: switch (user->sol) {
323: case SOL_QUADRATIC:
324: PetscDSSetResidual(ds, 0, f0_quadratic_u, f1_u);
325: exactFuncs[0] = quadratic_u;
326: exactFuncs[1] = quadratic_p;
327: break;
328: case SOL_TRIG:
329: PetscDSSetResidual(ds, 0, f0_trig_u, f1_u);
330: exactFuncs[0] = trig_u;
331: exactFuncs[1] = trig_p;
332: break;
333: default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%D)", SolTypes[PetscMin(user->sol, SOL_UNKNOWN)], user->sol);
334: }
335: PetscDSSetResidual(ds, 1, f0_p, NULL);
336: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
337: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
338: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
339: PetscDSSetJacobianPreconditioner(ds, 0, 0, NULL, NULL, NULL, g3_uu);
340: PetscDSSetJacobianPreconditioner(ds, 1, 1, g0_pp, NULL, NULL, NULL);
342: PetscDSSetExactSolution(ds, 0, exactFuncs[0], user);
343: PetscDSSetExactSolution(ds, 1, exactFuncs[1], user);
345: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, 1, &id, user);
347: /* Make constant values available to pointwise functions */
348: {
349: Parameter *param;
350: PetscScalar constants[1];
352: PetscBagGetData(user->bag, (void **) ¶m);
353: constants[0] = param->mu; /* dynamic shear viscosity, Pa s */
354: PetscDSSetConstants(ds, 1, constants);
355: }
356: return(0);
357: }
359: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
360: {
361: PetscInt c;
362: for (c = 0; c < Nc; ++c) u[c] = 0.0;
363: return 0;
364: }
365: static PetscErrorCode one(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
366: {
367: PetscInt c;
368: for (c = 0; c < Nc; ++c) u[c] = 1.0;
369: return 0;
370: }
372: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
373: {
374: Vec vec;
375: PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx) = {zero, one};
376: PetscErrorCode ierr;
379: if (origField != 1) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Field %D should be 1 for pressure", origField);
380: funcs[field] = one;
381: {
382: PetscDS ds;
383: DMGetDS(dm, &ds);
384: PetscObjectViewFromOptions((PetscObject) ds, NULL, "-ds_view");
385: }
386: DMCreateGlobalVector(dm, &vec);
387: DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);
388: VecNormalize(vec, NULL);
389: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullspace);
390: VecDestroy(&vec);
391: /* New style for field null spaces */
392: {
393: PetscObject pressure;
394: MatNullSpace nullspacePres;
396: DMGetField(dm, field, NULL, &pressure);
397: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres);
398: PetscObjectCompose(pressure, "nullspace", (PetscObject) nullspacePres);
399: MatNullSpaceDestroy(&nullspacePres);
400: }
401: return(0);
402: }
404: static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
405: {
406: DM cdm = dm;
407: PetscQuadrature q = NULL;
408: DMPolytopeType ct;
409: PetscBool simplex;
410: PetscInt dim, Nf = 2, f, Nc[2], cStart;
411: const char *name[2] = {"velocity", "pressure"};
412: const char *prefix[2] = {"vel_", "pres_"};
413: PetscErrorCode ierr;
416: DMGetDimension(dm, &dim);
417: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
418: DMPlexGetCellType(dm, cStart, &ct);
419: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
420: Nc[0] = dim;
421: Nc[1] = 1;
422: for (f = 0; f < Nf; ++f) {
423: PetscFE fe;
425: PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe);
426: PetscObjectSetName((PetscObject) fe, name[f]);
427: if (!q) {PetscFEGetQuadrature(fe, &q);}
428: PetscFESetQuadrature(fe, q);
429: DMSetField(dm, f, NULL, (PetscObject) fe);
430: PetscFEDestroy(&fe);
431: }
432: DMCreateDS(dm);
433: (*setupEqn)(dm, user);
434: while (cdm) {
435: DMCopyDisc(dm, cdm);
436: DMSetNullSpaceConstructor(cdm, 1, CreatePressureNullSpace);
437: DMGetCoarseDM(cdm, &cdm);
438: }
439: return(0);
440: }
442: int main(int argc, char **argv)
443: {
444: SNES snes;
445: DM dm;
446: Vec u;
447: AppCtx user;
450: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
451: ProcessOptions(PETSC_COMM_WORLD, &user);
452: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
453: SNESCreate(PetscObjectComm((PetscObject) dm), &snes);
454: SNESSetDM(snes, dm);
455: DMSetApplicationContext(dm, &user);
457: SetupParameters(PETSC_COMM_WORLD, &user);
458: SetupProblem(dm, SetupEqn, &user);
459: DMPlexCreateClosureIndex(dm, NULL);
461: DMCreateGlobalVector(dm, &u);
462: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
463: SNESSetFromOptions(snes);
464: DMSNESCheckFromOptions(snes, u);
465: PetscObjectSetName((PetscObject) u, "Solution");
466: {
467: Mat J;
468: MatNullSpace sp;
470: SNESSetUp(snes);
471: CreatePressureNullSpace(dm, 1, 1, &sp);
472: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
473: MatSetNullSpace(J, sp);
474: MatNullSpaceDestroy(&sp);
475: PetscObjectSetName((PetscObject) J, "Jacobian");
476: MatViewFromOptions(J, NULL, "-J_view");
477: }
478: SNESSolve(snes, NULL, u);
480: VecDestroy(&u);
481: SNESDestroy(&snes);
482: DMDestroy(&dm);
483: PetscBagDestroy(&user.bag);
484: PetscFinalize();
485: return ierr;
486: }
487: /*TEST
489: test:
490: suffix: 2d_p2_p1_check
491: requires: triangle
492: args: -sol quadratic -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
494: test:
495: suffix: 2d_p2_p1_check_parallel
496: nsize: {{2 3 5}}
497: requires: triangle
498: args: -sol quadratic -dm_refine 2 -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
500: test:
501: suffix: 3d_p2_p1_check
502: requires: ctetgen
503: args: -sol quadratic -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
505: test:
506: suffix: 3d_p2_p1_check_parallel
507: nsize: {{2 3 5}}
508: requires: ctetgen
509: args: -sol quadratic -dm_refine 2 -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
511: test:
512: suffix: 2d_p2_p1_conv
513: requires: triangle
514: # Using -dm_refine 3 gives L_2 convergence rate: [3.0, 2.1]
515: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
516: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
517: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
518: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
520: test:
521: suffix: 2d_p2_p1_conv_gamg
522: requires: triangle
523: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
524: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
525: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_coarse_pc_type svd
527: test:
528: suffix: 3d_p2_p1_conv
529: requires: ctetgen !single
530: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.8]
531: args: -sol trig -dm_plex_box_dim 3 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
532: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
533: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
534: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
536: test:
537: suffix: 2d_q2_q1_check
538: args: -sol quadratic -dm_plex_box_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
540: test:
541: suffix: 3d_q2_q1_check
542: args: -sol quadratic -dm_plex_box_simplex 0 -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
544: test:
545: suffix: 2d_q2_q1_conv
546: # Using -dm_refine 3 -convest_num_refine 1 gives L_2 convergence rate: [3.0, 2.1]
547: args: -sol trig -dm_plex_box_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -ksp_error_if_not_converged \
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: 3d_q2_q1_conv
554: requires: !single
555: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.4]
556: args: -sol trig -dm_plex_box_simplex 0 -dm_plex_box_dim 3 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
557: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
558: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
559: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
561: test:
562: suffix: 2d_p3_p2_check
563: requires: triangle
564: args: -sol quadratic -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
566: test:
567: suffix: 3d_p3_p2_check
568: requires: ctetgen !single
569: args: -sol quadratic -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
571: test:
572: suffix: 2d_p3_p2_conv
573: requires: triangle
574: # Using -dm_refine 2 gives L_2 convergence rate: [3.8, 3.0]
575: args: -sol trig -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
576: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
577: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
578: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
580: test:
581: suffix: 3d_p3_p2_conv
582: requires: ctetgen long_runtime
583: # Using -dm_refine 1 -convest_num_refine 2 gives L_2 convergence rate: [3.6, 3.9]
584: args: -sol trig -dm_plex_box_dim 3 -dm_refine 1 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 \
585: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
586: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
587: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
589: test:
590: suffix: 2d_q1_p0_conv
591: requires: !single
592: # Using -dm_refine 3 gives L_2 convergence rate: [1.9, 1.0]
593: args: -sol quadratic -dm_plex_box_simplex 0 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 2 \
594: -ksp_atol 1e-10 -ksp_error_if_not_converged -petscds_jac_pre 0 \
595: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
596: -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
598: test:
599: suffix: 3d_q1_p0_conv
600: requires: !single
601: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [1.7, 1.0]
602: args: -sol quadratic -dm_plex_box_simplex 0 -dm_plex_box_dim 3 -dm_refine 1 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 1 \
603: -ksp_atol 1e-10 -ksp_error_if_not_converged -petscds_jac_pre 0 \
604: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
605: -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
607: # Stokes preconditioners
608: # Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
609: test:
610: suffix: 2d_p2_p1_block_diagonal
611: requires: triangle
612: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
613: -snes_error_if_not_converged \
614: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-4 -ksp_error_if_not_converged \
615: -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
616: # Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
617: test:
618: suffix: 2d_p2_p1_block_triangular
619: requires: triangle
620: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
621: -snes_error_if_not_converged \
622: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
623: -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
624: # Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
625: test:
626: suffix: 2d_p2_p1_schur_diagonal
627: requires: triangle
628: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
629: -snes_error_if_not_converged \
630: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
631: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -pc_fieldsplit_off_diag_use_amat \
632: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
633: # Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
634: test:
635: suffix: 2d_p2_p1_schur_upper
636: requires: triangle
637: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001 \
638: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
639: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -pc_fieldsplit_off_diag_use_amat \
640: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
641: # Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
642: test:
643: suffix: 2d_p2_p1_schur_lower
644: requires: triangle
645: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
646: -snes_error_if_not_converged \
647: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
648: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -pc_fieldsplit_off_diag_use_amat \
649: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
650: # 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}
651: test:
652: suffix: 2d_p2_p1_schur_full
653: requires: triangle
654: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
655: -snes_error_if_not_converged \
656: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
657: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_off_diag_use_amat \
658: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
659: # Full Schur + Velocity GMG
660: test:
661: suffix: 2d_p2_p1_gmg_vcycle
662: requires: triangle
663: args: -sol quadratic -dm_refine_hierarchy 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
664: -snes_error_if_not_converged \
665: -ksp_type fgmres -ksp_atol 1e-9 -ksp_error_if_not_converged -pc_use_amat \
666: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_off_diag_use_amat \
667: -fieldsplit_velocity_pc_type mg -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_coarse_pc_type svd
668: # 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}
669: test:
670: suffix: 2d_p2_p1_simple
671: requires: triangle
672: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
673: -snes_error_if_not_converged \
674: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
675: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
676: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
677: -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi
678: # FETI-DP solvers (these solvers are quite inefficient, they are here to exercise the code)
679: test:
680: suffix: 2d_p2_p1_fetidp
681: requires: triangle mumps
682: nsize: 5
683: args: -sol quadratic -dm_refine 2 -dm_mat_type is -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
684: -snes_error_if_not_converged \
685: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
686: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
687: -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 \
688: -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
689: test:
690: suffix: 2d_q2_q1_fetidp
691: requires: mumps
692: nsize: 5
693: args: -sol quadratic -dm_plex_box_simplex 0 -dm_refine 2 -dm_mat_type is -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
694: -snes_error_if_not_converged \
695: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
696: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
697: -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 \
698: -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
699: test:
700: suffix: 3d_p2_p1_fetidp
701: requires: ctetgen mumps suitesparse
702: nsize: 5
703: args: -sol quadratic -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
704: -snes_error_if_not_converged \
705: -ksp_type fetidp -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
706: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
707: -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 \
708: -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_benign_trick -fetidp_bddc_pc_bddc_deluxe_singlemat \
709: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
710: -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 \
711: -fetidp_bddelta_pc_factor_mat_ordering_type external \
712: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
713: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
714: test:
715: suffix: 3d_q2_q1_fetidp
716: requires: suitesparse
717: nsize: 5
718: args: -sol quadratic -dm_plex_box_simplex 0 -dm_plex_box_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
719: -snes_error_if_not_converged \
720: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
721: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
722: -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 \
723: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
724: -fetidp_bddc_pc_bddc_symmetric -fetidp_fieldsplit_lag_ksp_type preonly \
725: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
726: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
727: # BDDC solvers (these solvers are quite inefficient, they are here to exercise the code)
728: test:
729: suffix: 2d_p2_p1_bddc
730: nsize: 2
731: requires: triangle !single
732: args: -sol quadratic -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -dm_distribute -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
733: -snes_error_if_not_converged \
734: -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
735: -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
736: # Vanka
737: test:
738: suffix: 2d_q1_p0_vanka
739: requires: double !complex
740: args: -sol quadratic -dm_plex_box_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
741: -snes_rtol 1.0e-4 \
742: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
743: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
744: -sub_ksp_type preonly -sub_pc_type lu
745: test:
746: suffix: 2d_q1_p0_vanka_denseinv
747: requires: double !complex
748: args: -sol quadratic -dm_plex_box_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
749: -snes_rtol 1.0e-4 \
750: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
751: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
752: -pc_patch_dense_inverse -pc_patch_sub_mat_type seqdense
753: # Vanka smoother
754: test:
755: suffix: 2d_q1_p0_gmg_vanka
756: requires: double !complex
757: args: -sol quadratic -dm_plex_box_simplex 0 -dm_refine_hierarchy 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
758: -snes_rtol 1.0e-4 \
759: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
760: -pc_type mg \
761: -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 30 \
762: -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 \
763: -mg_levels_sub_ksp_type preonly -mg_levels_sub_pc_type lu \
764: -mg_coarse_pc_type svd
766: TEST*/