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,
34: SOL_QUADRATIC,
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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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: f0[0] = 0.0;
156: for (PetscInt d = 0; d < dim; ++d) f0[0] += 2.0 * x[d] * (1.0 - x[d]);
157: 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);
158: }
160: /* 2D Dirichlet potential example
162: u = x (1-x) y (1-y)
163: q = <(1-2x) y (1-y), x (1-x) (1-2y)>
164: f = -y (1-y) - x (1-x)
166: du/dn_\Gamma = {(1-2x) y (1-y), x (1-x) (1-2y)} . n produces an essential condition on q
167: */
169: 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[])
170: {
171: PetscInt c;
172: for (c = 0; c < dim; ++c) f0[c] = u[uOff[0] + c];
173: }
175: /* <\nabla\cdot w, u> = <\nabla w, Iu> */
176: 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[])
177: {
178: for (PetscInt c = 0; c < dim; ++c) f1[c * dim + c] = u[uOff[1]];
179: }
181: /* <t, q> */
182: 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[])
183: {
184: for (PetscInt c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
185: }
187: /* <\nabla\cdot t, u> = <\nabla t, Iu> */
188: 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[])
189: {
190: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
191: }
193: /* <v, \nabla\cdot q> */
194: 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[])
195: {
196: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
197: }
199: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
200: {
201: PetscFunctionBeginUser;
202: options->solType = SOL_LINEAR;
203: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
204: PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex24.c", SolTypeNames, (PetscEnum)options->solType, (PetscEnum *)&options->solType, NULL));
205: PetscOptionsEnd();
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
210: {
211: PetscFunctionBeginUser;
212: PetscCall(DMCreate(comm, dm));
213: PetscCall(DMSetType(*dm, DMPLEX));
214: PetscCall(PetscObjectSetName((PetscObject)*dm, "Example Mesh"));
215: PetscCall(DMSetApplicationContext(*dm, user));
216: PetscCall(DMSetFromOptions(*dm));
217: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
222: {
223: PetscDS ds;
224: DMLabel label;
225: PetscWeakForm wf;
226: const PetscInt id = 1;
227: PetscInt bd;
229: PetscFunctionBeginUser;
230: PetscCall(DMGetLabel(dm, "marker", &label));
231: PetscCall(DMGetDS(dm, &ds));
232: PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
233: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
234: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL));
235: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL));
236: switch (user->solType) {
237: case SOL_LINEAR:
238: PetscCall(PetscDSSetResidual(ds, 1, f0_linear_u, NULL));
239: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
240: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
241: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_linear_q, 0, NULL));
242: PetscCall(PetscDSSetExactSolution(ds, 0, linear_q, user));
243: PetscCall(PetscDSSetExactSolution(ds, 1, linear_u, user));
244: break;
245: case SOL_QUADRATIC:
246: PetscCall(PetscDSSetResidual(ds, 1, f0_quadratic_u, NULL));
247: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
248: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
249: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_quadratic_q, 0, NULL));
250: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user));
251: PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user));
252: break;
253: case SOL_QUARTIC:
254: PetscCall(PetscDSSetResidual(ds, 1, f0_quartic_u, NULL));
255: PetscCall(PetscDSSetExactSolution(ds, 0, quartic_q, user));
256: PetscCall(PetscDSSetExactSolution(ds, 1, quartic_u, user));
257: break;
258: case SOL_QUARTIC_NEUMANN:
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: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "Flux condition", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quartic_q, NULL, user, NULL));
263: break;
264: default:
265: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static PetscErrorCode SetupDiscretization(DM dm, PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
271: {
272: DM cdm = dm;
273: PetscFE feq, feu;
274: DMPolytopeType ct;
275: PetscBool simplex;
276: PetscInt dim, cStart;
278: PetscFunctionBeginUser;
279: PetscCall(DMGetDimension(dm, &dim));
280: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
281: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
282: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
283: /* Create finite element */
284: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", -1, &feq));
285: PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
286: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", -1, &feu));
287: PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
288: PetscCall(PetscFECopyQuadrature(feq, feu));
289: /* Set discretization and boundary conditions for each mesh */
290: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
291: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu));
292: PetscCall(DMCreateDS(dm));
293: PetscCall((*setup)(dm, user));
294: while (cdm) {
295: PetscCall(DMCopyDisc(dm, cdm));
296: PetscCall(DMGetCoarseDM(cdm, &cdm));
297: }
298: PetscCall(PetscFEDestroy(&feq));
299: PetscCall(PetscFEDestroy(&feu));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: int main(int argc, char **argv)
304: {
305: DM dm; /* Problem specification */
306: SNES snes; /* Nonlinear solver */
307: Vec u; /* Solutions */
308: AppCtx user; /* User-defined work context */
310: PetscFunctionBeginUser;
311: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
312: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
313: /* Primal system */
314: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
315: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
316: PetscCall(SNESSetDM(snes, dm));
317: PetscCall(SetupDiscretization(dm, SetupPrimalProblem, &user));
318: PetscCall(DMCreateGlobalVector(dm, &u));
319: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
320: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
321: PetscCall(SNESSetFromOptions(snes));
322: PetscCall(DMSNESCheckFromOptions(snes, u));
323: PetscCall(SNESSolve(snes, NULL, u));
324: PetscCall(SNESGetSolution(snes, &u));
325: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
326: /* Cleanup */
327: PetscCall(VecDestroy(&u));
328: PetscCall(SNESDestroy(&snes));
329: PetscCall(DMDestroy(&dm));
330: PetscCall(PetscFinalize());
331: return 0;
332: }
334: /*TEST
336: test:
337: suffix:2d_rt0_tri
338: requires: triangle
339: args: -sol_type linear -dmsnes_check 0.001 \
340: -potential_petscspace_degree 0 \
341: -potential_petscdualspace_lagrange_continuity 0 \
342: -field_petscspace_type ptrimmed \
343: -field_petscspace_components 2 \
344: -field_petscspace_ptrimmed_form_degree -1 \
345: -field_petscdualspace_order 1 \
346: -field_petscdualspace_form_degree -1 \
347: -field_petscdualspace_lagrange_trimmed true \
348: -field_petscfe_default_quadrature_order 2 \
349: -snes_error_if_not_converged \
350: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
351: -pc_type fieldsplit -pc_fieldsplit_type schur \
352: -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
353: -fieldsplit_field_pc_type lu \
354: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
356: test:
357: suffix:2d_rt0_quad
358: requires: triangle
359: args: -dm_plex_simplex 0 -sol_type linear -dmsnes_check 0.001 \
360: -potential_petscspace_degree 0 \
361: -potential_petscdualspace_lagrange_continuity 0 \
362: -field_petscspace_degree 1 \
363: -field_petscspace_type sum \
364: -field_petscspace_variables 2 \
365: -field_petscspace_components 2 \
366: -field_petscspace_sum_spaces 2 \
367: -field_petscspace_sum_concatenate true \
368: -field_sumcomp_0_petscspace_variables 2 \
369: -field_sumcomp_0_petscspace_type tensor \
370: -field_sumcomp_0_petscspace_tensor_spaces 2 \
371: -field_sumcomp_0_petscspace_tensor_uniform false \
372: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
373: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
374: -field_sumcomp_1_petscspace_variables 2 \
375: -field_sumcomp_1_petscspace_type tensor \
376: -field_sumcomp_1_petscspace_tensor_spaces 2 \
377: -field_sumcomp_1_petscspace_tensor_uniform false \
378: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
379: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
380: -field_petscdualspace_form_degree -1 \
381: -field_petscdualspace_order 1 \
382: -field_petscdualspace_lagrange_trimmed true \
383: -field_petscfe_default_quadrature_order 2 \
384: -pc_type fieldsplit -pc_fieldsplit_type schur \
385: -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
386: -fieldsplit_field_pc_type lu \
387: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
389: test:
390: suffix: 2d_bdm1_p0
391: requires: triangle
392: args: -sol_type linear \
393: -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 1 \
394: -dmsnes_check .001 -snes_error_if_not_converged \
395: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
396: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
397: -fieldsplit_field_pc_type lu \
398: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
399: test:
400: nsize: 4
401: suffix: 2d_bdm1_p0_bddc
402: requires: triangle !single
403: args: -sol_type linear \
404: -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 1 \
405: -dmsnes_check .001 -snes_error_if_not_converged \
406: -ksp_error_if_not_converged -ksp_type cg -ksp_norm_type natural -ksp_divtol 1e10 \
407: -petscpartitioner_type simple -dm_mat_type is \
408: -pc_type bddc -pc_bddc_use_local_mat_graph 0 \
409: -pc_bddc_benign_trick -pc_bddc_nonetflux -pc_bddc_detect_disconnected -pc_bddc_use_qr_single \
410: -pc_bddc_coarse_redundant_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_dirichlet_pc_type svd
412: test:
413: nsize: 9
414: requires: !single
415: suffix: 2d_rt1_p0_bddc
416: args: -sol_type quadratic \
417: -potential_petscspace_degree 0 \
418: -potential_petscdualspace_lagrange_use_moments \
419: -potential_petscdualspace_lagrange_moment_order 2 \
420: -field_petscfe_default_quadrature_order 1 \
421: -field_petscspace_degree 1 \
422: -field_petscspace_type sum \
423: -field_petscspace_variables 2 \
424: -field_petscspace_components 2 \
425: -field_petscspace_sum_spaces 2 \
426: -field_petscspace_sum_concatenate true \
427: -field_sumcomp_0_petscspace_variables 2 \
428: -field_sumcomp_0_petscspace_type tensor \
429: -field_sumcomp_0_petscspace_tensor_spaces 2 \
430: -field_sumcomp_0_petscspace_tensor_uniform false \
431: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
432: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
433: -field_sumcomp_1_petscspace_variables 2 \
434: -field_sumcomp_1_petscspace_type tensor \
435: -field_sumcomp_1_petscspace_tensor_spaces 2 \
436: -field_sumcomp_1_petscspace_tensor_uniform false \
437: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
438: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
439: -field_petscdualspace_form_degree -1 \
440: -field_petscdualspace_order 1 \
441: -field_petscdualspace_lagrange_trimmed true \
442: -dm_plex_box_faces 3,3 dm_refine 1 -dm_plex_simplex 0 \
443: -dmsnes_check .001 -snes_error_if_not_converged \
444: -ksp_error_if_not_converged -ksp_type cg \
445: -petscpartitioner_type simple -dm_mat_type is \
446: -pc_type bddc -pc_bddc_use_local_mat_graph 0 \
447: -pc_bddc_benign_trick -pc_bddc_nonetflux -pc_bddc_detect_disconnected -pc_bddc_use_qr_single \
448: -pc_bddc_coarse_redundant_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_dirichlet_pc_type svd
450: test:
451: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.0, 1.0]
452: # Using -sol_type quadratic -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.9, 1.0]
453: suffix: 2d_bdm1_p0_conv
454: requires: triangle
455: args: -sol_type quartic \
456: -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
457: -snes_error_if_not_converged \
458: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
459: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
460: -fieldsplit_field_pc_type lu \
461: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
463: test:
464: # 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
465: # VTK output: -potential_view vtk: -exact_vec_view vtk:
466: # VTU output: -potential_view vtk:multifield.vtu -exact_vec_view vtk:exact.vtu
467: suffix: 2d_q2_p0
468: requires: triangle
469: args: -sol_type linear -dm_plex_simplex 0 \
470: -field_petscspace_degree 2 -dm_refine 0 \
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
476: test:
477: # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [0.0057, 1.0]
478: suffix: 2d_q2_p0_conv
479: requires: triangle
480: args: -sol_type linear -dm_plex_simplex 0 \
481: -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
482: -snes_error_if_not_converged \
483: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
484: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
485: -fieldsplit_field_pc_type lu \
486: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
487: test:
488: # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [-0.014, 0.0066]
489: suffix: 2d_q2_p0_neumann_conv
490: requires: triangle
491: args: -sol_type quartic_neumann -dm_plex_simplex 0 \
492: -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \
493: -snes_error_if_not_converged \
494: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
495: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
496: -fieldsplit_field_pc_type lu \
497: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type svd
499: TEST*/
501: /* These tests will be active once tensor P^- is working
503: test:
504: suffix: 2d_bdmq1_p0_0
505: requires: triangle
506: args: -dm_plex_simplex 0 -sol_type linear \
507: -field_petscspace_poly_type pminus_hdiv -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate \
508: -dmsnes_check .001 -snes_error_if_not_converged \
509: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
510: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
511: -fieldsplit_field_pc_type lu \
512: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
513: test:
514: suffix: 2d_bdmq1_p0_2
515: requires: triangle
516: args: -dm_plex_simplex 0 -sol_type quartic \
517: -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 \
518: -dmsnes_check .001 -snes_error_if_not_converged \
519: -ksp_rtol 1e-10 -ksp_error_if_not_converged \
520: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \
521: -fieldsplit_field_pc_type lu \
522: -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu
524: */