Actual source code: ex24.c

  1: static char help[] = "Poisson Problem in mixed form with 2d and 3d with finite elements.\n\
2: We solve the Poisson problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example supports automatic convergence estimation\n\
5: and Hdiv elements.\n\n\n";

7: /*

9: The mixed form of Poisson's equation, e.g. https://firedrakeproject.org/demos/poisson_mixed.py.html, is given
10: in the strong form by
11: \begin{align}
12:   \vb{q} - \nabla u   &= 0 \\
13:   \nabla \cdot \vb{q} &= f
14: \end{align}
15: where $u$ is the potential, as in the original problem, but we also discretize the gradient of potential $\vb{q}$,
16: or flux, directly. The weak form is then
17: \begin{align}
18:   <t, \vb{q}> + <\nabla \cdot t, u> - <t_n, u>_\Gamma &= 0 \\
19:   <v, \nabla \cdot \vb{q}> &= <v, f>
20: \end{align}
21: For the original Poisson problem, the Dirichlet boundary forces an essential boundary condition on the potential space,
22: and the Neumann boundary gives a natural boundary condition in the weak form. In the mixed formulation, the Neumann
23: boundary gives an essential boundary condition on the flux space, $\vb{q} \cdot \vu{n} = h$, and the Dirichlet condition
24: becomes a natural condition in the weak form, <t_n, g>_\Gamma.
25: */

27: #include <petscdmplex.h>
28: #include <petscsnes.h>
29: #include <petscds.h>
30: #include <petscconvest.h>

32: typedef enum {
33:   SOL_LINEAR,
35:   SOL_QUARTIC,
36:   SOL_QUARTIC_NEUMANN,
37:   SOL_UNKNOWN,
38:   NUM_SOLTYPE
39: } SolType;
40: const char *SolTypeNames[NUM_SOLTYPE + 3] = {"linear", "quadratic", "quartic", "quartic_neumann", "unknown", "SolType", "SOL_", NULL};

42: typedef struct {
43:   SolType solType; /* The type of exact solution */
44: } AppCtx;

46: static void f0_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[])
47: {
48:   PetscInt d;
49:   for (d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
50: }

52: /* 2D Dirichlet potential example

54:   u = x
55:   q = <1, 0>
56:   f = 0

58:   We will need a boundary integral of u over \Gamma.
59: */
60: static PetscErrorCode linear_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
61: {
62:   u[0] = x[0];
63:   return PETSC_SUCCESS;
64: }
65: static PetscErrorCode linear_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
66: {
67:   PetscInt c;
68:   for (c = 0; c < Nc; ++c) u[c] = c ? 0.0 : 1.0;
69:   return PETSC_SUCCESS;
70: }

72: static void f0_linear_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[])
73: {
74:   f0[0] = 0.0;
75:   f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0);
76: }
77: static void f0_bd_linear_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
78: {
79:   PetscScalar potential;
80:   PetscInt    d;

82:   PetscCallAbort(PETSC_COMM_SELF, linear_u(dim, t, x, dim, &potential, NULL));
83:   for (d = 0; d < dim; ++d) f0[d] = -potential * n[d];
84: }

86: /* 2D Dirichlet potential example

88:   u = x^2 + y^2
89:   q = <2x, 2y>
90:   f = 4

92:   We will need a boundary integral of u over \Gamma.
93: */
94: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
95: {
96:   PetscInt d;

98:   u[0] = 0.0;
99:   for (d = 0; d < dim; ++d) u[0] += x[d] * x[d];
100:   return PETSC_SUCCESS;
101: }
102: static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
103: {
104:   PetscInt c;
105:   for (c = 0; c < Nc; ++c) u[c] = 2.0 * x[c];
106:   return PETSC_SUCCESS;
107: }

109: 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[])
110: {
111:   f0[0] = -4.0;
112:   f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0);
113: }
114: static void f0_bd_quadratic_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
115: {
116:   PetscScalar potential;
117:   PetscInt    d;

119:   PetscCallAbort(PETSC_COMM_SELF, quadratic_u(dim, t, x, dim, &potential, NULL));
120:   for (d = 0; d < dim; ++d) f0[d] = -potential * n[d];
121: }

123: /* 2D Dirichlet potential example

125:   u = x (1-x) y (1-y)
126:   q = <(1-2x) y (1-y), x (1-x) (1-2y)>
127:   f = -y (1-y) - x (1-x)

129:   u|_\Gamma = 0 so that the boundary integral vanishes
130: */
131: static PetscErrorCode quartic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
132: {
133:   PetscInt d;

135:   u[0] = 1.0;
136:   for (d = 0; d < dim; ++d) u[0] *= x[d] * (1.0 - x[d]);
137:   return PETSC_SUCCESS;
138: }
139: static PetscErrorCode quartic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
140: {
141:   PetscInt c, d;

143:   for (c = 0; c < Nc; ++c) {
144:     u[c] = 1.0;
145:     for (d = 0; d < dim; ++d) {
146:       if (c == d) u[c] *= 1 - 2.0 * x[d];
147:       else u[c] *= x[d] * (1.0 - x[d]);
148:     }
149:   }
150:   return PETSC_SUCCESS;
151: }

153: static void f0_quartic_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[])
154: {
155:   PetscInt d;
156:   f0[0] = 0.0;
157:   for (d = 0; d < dim; ++d) f0[0] += 2.0 * x[d] * (1.0 - x[d]);
158:   f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0);
159: }

161: /* 2D Dirichlet potential example

163:   u = x (1-x) y (1-y)
164:   q = <(1-2x) y (1-y), x (1-x) (1-2y)>
165:   f = -y (1-y) - x (1-x)

167:   du/dn_\Gamma = {(1-2x) y (1-y), x (1-x) (1-2y)} . n produces an essential condition on q
168: */

170: static void f0_q(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[])
171: {
172:   PetscInt c;
173:   for (c = 0; c < dim; ++c) f0[c] = u[uOff[0] + c];
174: }

176: /* <\nabla\cdot w, u> = <\nabla w, Iu> */
177: static void f1_q(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[])
178: {
179:   PetscInt c;
180:   for (c = 0; c < dim; ++c) f1[c * dim + c] = u[uOff[1]];
181: }

183: /* <t, q> */
184: static void g0_qq(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[])
185: {
186:   PetscInt c;
187:   for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
188: }

190: /* <\nabla\cdot t, u> = <\nabla t, Iu> */
191: static void g2_qu(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[])
192: {
193:   PetscInt d;
194:   for (d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
195: }

197: /* <v, \nabla\cdot q> */
198: static void g1_uq(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[])
199: {
200:   PetscInt d;
201:   for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
202: }

204: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
205: {
206:   PetscFunctionBeginUser;
207:   options->solType = SOL_LINEAR;
208:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
209:   PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex24.c", SolTypeNames, (PetscEnum)options->solType, (PetscEnum *)&options->solType, NULL));
210:   PetscOptionsEnd();
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
215: {
216:   PetscFunctionBeginUser;
217:   PetscCall(DMCreate(comm, dm));
218:   PetscCall(DMSetType(*dm, DMPLEX));
219:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Example Mesh"));
220:   PetscCall(DMSetApplicationContext(*dm, user));
221:   PetscCall(DMSetFromOptions(*dm));
222:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
227: {
228:   PetscDS        ds;
229:   DMLabel        label;
230:   PetscWeakForm  wf;
231:   const PetscInt id = 1;
232:   PetscInt       bd;

234:   PetscFunctionBeginUser;
235:   PetscCall(DMGetLabel(dm, "marker", &label));
236:   PetscCall(DMGetDS(dm, &ds));
237:   PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
238:   PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
239:   PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL));
240:   PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL));
241:   switch (user->solType) {
242:   case SOL_LINEAR:
243:     PetscCall(PetscDSSetResidual(ds, 1, f0_linear_u, NULL));
244:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
245:     PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
246:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_linear_q, 0, NULL));
247:     PetscCall(PetscDSSetExactSolution(ds, 0, linear_q, user));
248:     PetscCall(PetscDSSetExactSolution(ds, 1, linear_u, user));
249:     break;
252:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
253:     PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
254:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_quadratic_q, 0, NULL));
257:     break;
258:   case SOL_QUARTIC:
259:     PetscCall(PetscDSSetResidual(ds, 1, f0_quartic_u, NULL));
260:     PetscCall(PetscDSSetExactSolution(ds, 0, quartic_q, user));
261:     PetscCall(PetscDSSetExactSolution(ds, 1, quartic_u, user));
262:     break;
263:   case SOL_QUARTIC_NEUMANN:
264:     PetscCall(PetscDSSetResidual(ds, 1, f0_quartic_u, NULL));
265:     PetscCall(PetscDSSetExactSolution(ds, 0, quartic_q, user));
266:     PetscCall(PetscDSSetExactSolution(ds, 1, quartic_u, user));
267:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "Flux condition", label, 1, &id, 0, 0, NULL, (void (*)(void))quartic_q, NULL, user, NULL));
268:     break;
269:   default:
270:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
271:   }
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: static PetscErrorCode SetupDiscretization(DM dm, PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
276: {
277:   DM             cdm = dm;
278:   PetscFE        feq, feu;
279:   DMPolytopeType ct;
280:   PetscBool      simplex;
281:   PetscInt       dim, cStart;

283:   PetscFunctionBeginUser;
284:   PetscCall(DMGetDimension(dm, &dim));
285:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
286:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
287:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
288:   /* Create finite element */
289:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", -1, &feq));
290:   PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
291:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", -1, &feu));
292:   PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
294:   /* Set discretization and boundary conditions for each mesh */
295:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
296:   PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu));
297:   PetscCall(DMCreateDS(dm));
298:   PetscCall((*setup)(dm, user));
299:   while (cdm) {
300:     PetscCall(DMCopyDisc(dm, cdm));
301:     PetscCall(DMGetCoarseDM(cdm, &cdm));
302:   }
303:   PetscCall(PetscFEDestroy(&feq));
304:   PetscCall(PetscFEDestroy(&feu));
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }

308: int main(int argc, char **argv)
309: {
310:   DM     dm;   /* Problem specification */
311:   SNES   snes; /* Nonlinear solver */
312:   Vec    u;    /* Solutions */
313:   AppCtx user; /* User-defined work context */

315:   PetscFunctionBeginUser;
316:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
317:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
318:   /* Primal system */
319:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
320:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
321:   PetscCall(SNESSetDM(snes, dm));
322:   PetscCall(SetupDiscretization(dm, SetupPrimalProblem, &user));
323:   PetscCall(DMCreateGlobalVector(dm, &u));
324:   PetscCall(VecSet(u, 0.0));
325:   PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
326:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
327:   PetscCall(SNESSetFromOptions(snes));
328:   PetscCall(DMSNESCheckFromOptions(snes, u));
329:   PetscCall(SNESSolve(snes, NULL, u));
330:   PetscCall(SNESGetSolution(snes, &u));
331:   PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
332:   /* Cleanup */
333:   PetscCall(VecDestroy(&u));
334:   PetscCall(SNESDestroy(&snes));
335:   PetscCall(DMDestroy(&dm));
336:   PetscCall(PetscFinalize());
337:   return 0;
338: }

340: /*TEST

342:   test:
343:     suffix: 2d_bdm1_p0
344:     requires: triangle
345:     args: -sol_type linear \
346:           -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 1 \
347:           -dmsnes_check .001 -snes_error_if_not_converged \
348:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
349:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
350:             -fieldsplit_field_pc_type lu \
351:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
352:   test:
353:     nsize: 4
354:     suffix: 2d_bdm1_p0_bddc
355:     requires: triangle !single
356:     args: -sol_type linear \
357:           -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 1 \
358:           -dmsnes_check .001 -snes_error_if_not_converged \
359:           -ksp_error_if_not_converged -ksp_type cg \
360:           -petscpartitioner_type simple -dm_mat_type is \
361:           -pc_type bddc -pc_bddc_use_local_mat_graph 0 \
362:           -pc_bddc_benign_trick -pc_bddc_nonetflux -pc_bddc_detect_disconnected -pc_bddc_use_qr_single \
363:           -pc_bddc_coarse_redundant_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_dirichlet_pc_type svd

365:   test:
366:     nsize: 9
367:     requires: !single
368:     suffix: 2d_rt1_p0_bddc
370:           -potential_petscspace_degree 0 \
371:           -potential_petscdualspace_lagrange_use_moments \
372:           -potential_petscdualspace_lagrange_moment_order 2 \
374:           -field_petscspace_degree 1 \
375:           -field_petscspace_type sum \
376:           -field_petscspace_variables 2 \
377:           -field_petscspace_components 2 \
378:           -field_petscspace_sum_spaces 2 \
379:           -field_petscspace_sum_concatenate true \
380:           -field_sumcomp_0_petscspace_variables 2 \
381:           -field_sumcomp_0_petscspace_type tensor \
382:           -field_sumcomp_0_petscspace_tensor_spaces 2 \
383:           -field_sumcomp_0_petscspace_tensor_uniform false \
384:           -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
385:           -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
386:           -field_sumcomp_1_petscspace_variables 2 \
387:           -field_sumcomp_1_petscspace_type tensor \
388:           -field_sumcomp_1_petscspace_tensor_spaces 2 \
389:           -field_sumcomp_1_petscspace_tensor_uniform false \
390:           -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
391:           -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
392:           -field_petscdualspace_form_degree -1 \
393:           -field_petscdualspace_order 1 \
394:           -field_petscdualspace_lagrange_trimmed true \
395:           -dm_plex_box_faces 3,3 dm_refine 1 -dm_plex_simplex 0 \
396:           -dmsnes_check .001 -snes_error_if_not_converged \
397:           -ksp_error_if_not_converged -ksp_type cg \
398:           -petscpartitioner_type simple -dm_mat_type is \
399:           -pc_type bddc -pc_bddc_use_local_mat_graph 0 \
400:           -pc_bddc_benign_trick -pc_bddc_nonetflux -pc_bddc_detect_disconnected -pc_bddc_use_qr_single \
401:           -pc_bddc_coarse_redundant_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_dirichlet_pc_type svd

403:   test:
404:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.0, 1.0]
405:     # Using -sol_type quadratic -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.9, 1.0]
406:     suffix: 2d_bdm1_p0_conv
407:     requires: triangle
408:     args: -sol_type quartic \
409:           -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
410:           -snes_error_if_not_converged \
411:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
412:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
413:             -fieldsplit_field_pc_type lu \
414:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu

416:   test:
417:     # HDF5 output: -dm_view hdf5:${PETSC_DIR}/sol.h5 -potential_view hdf5:${PETSC_DIR}/sol.h5::append -exact_vec_view hdf5:\${PETSC_DIR}/sol.h5::append
418:     # VTK output: -potential_view vtk: -exact_vec_view vtk:
419:     # VTU output: -potential_view vtk:multifield.vtu -exact_vec_view vtk:exact.vtu
420:     suffix: 2d_q2_p0
421:     requires: triangle
422:     args: -sol_type linear -dm_plex_simplex 0 \
423:           -field_petscspace_degree 2 -dm_refine 0 \
424:           -dmsnes_check .001 -snes_error_if_not_converged \
425:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
426:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
427:             -fieldsplit_field_pc_type lu \
428:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
429:   test:
430:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [0.0057, 1.0]
431:     suffix: 2d_q2_p0_conv
432:     requires: triangle
433:     args: -sol_type linear -dm_plex_simplex 0 \
434:           -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
435:           -snes_error_if_not_converged \
436:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
437:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
438:             -fieldsplit_field_pc_type lu \
439:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
440:   test:
441:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [-0.014, 0.0066]
442:     suffix: 2d_q2_p0_neumann_conv
443:     requires: triangle
444:     args: -sol_type quartic_neumann -dm_plex_simplex 0 \
445:           -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
446:           -snes_error_if_not_converged \
447:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
448:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
449:             -fieldsplit_field_pc_type lu \
450:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type svd

452: TEST*/

454: /* These tests will be active once tensor P^- is working

456:   test:
457:     suffix: 2d_bdmq1_p0_0
458:     requires: triangle
459:     args: -dm_plex_simplex 0 -sol_type linear \
460:           -field_petscspace_poly_type pminus_hdiv -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate \
461:           -dmsnes_check .001 -snes_error_if_not_converged \
462:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
463:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
464:             -fieldsplit_field_pc_type lu \
465:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
466:   test:
467:     suffix: 2d_bdmq1_p0_2
468:     requires: triangle
469:     args: -dm_plex_simplex 0 -sol_type quartic \
470:           -field_petscspace_poly_type_no pminus_hdiv -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate \
471:           -dmsnes_check .001 -snes_error_if_not_converged \
472:           -ksp_rtol 1e-10 -ksp_error_if_not_converged \
473:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
474:             -fieldsplit_field_pc_type lu \
475:             -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu

477: */