Actual source code: ex30.c
1: static char help[] = "Biological network from https://link.springer.com/article/10.1007/s42967-023-00297-3\n\n\n";
3: #include <petscts.h>
4: #include <petscsf.h>
5: #include <petscdmplex.h>
6: #include <petscdmplextransform.h>
7: #include <petscdmforest.h>
8: #include <petscviewerhdf5.h>
9: #include <petscds.h>
11: /*
12: Here we solve the system of PDEs on \Omega \in R^2:
14: * dC/dt - D^2 \Delta C - c^2 \nabla p \cross \nabla p + \alpha sqrt(||C||^2_F + eps)^(\gamma-2) C = 0
15: * - \nabla \cdot ((r + C) \nabla p) = S
17: where:
18: C = symmetric 2x2 conductivity tensor
19: p = potential
20: S = source
22: with natural boundary conditions on \partial\Omega:
23: \nabla C \cdot n = 0
24: \nabla ((r + C)\nabla p) \cdot n = 0
26: Parameters:
27: D = diffusion constant
28: c = activation parameter
29: \alpha = metabolic coefficient
30: \gamma = metabolic exponent
31: r, eps are regularization parameters
33: We use Lagrange elements for C_ij and P.
34: */
36: typedef enum _fieldidx {
37: C_FIELD_ID = 0,
38: P_FIELD_ID,
39: NUM_FIELDS
40: } FieldIdx;
42: typedef enum _constantidx {
43: R_ID = 0,
44: EPS_ID,
45: ALPHA_ID,
46: GAMMA_ID,
47: D_ID,
48: C2_ID,
49: FIXC_ID,
50: NUM_CONSTANTS
51: } ConstantIdx;
53: PetscLogStage SetupStage, SolveStage;
55: #define NORM2C(c00, c01, c11) PetscSqr(c00) + 2 * PetscSqr(c01) + PetscSqr(c11)
57: /* residual for C when tested against basis functions */
58: static void C_0(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[])
59: {
60: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
61: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
62: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
63: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
64: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
65: const PetscScalar crossp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[1] * gradp[1]};
66: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
67: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
68: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
69: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
70: const PetscScalar nexp = (gamma - 2.0) / 2.0;
71: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
73: for (PetscInt k = 0; k < 3; k++) f0[k] = u_t[uOff[C_FIELD_ID] + k] - c2 * crossp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k];
74: }
76: /* Jacobian for C against C basis functions */
77: static void JC_0_c0c0(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 J[])
78: {
79: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
80: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
81: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
82: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
83: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
84: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
85: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
86: const PetscScalar nexp = (gamma - 2.0) / 2.0;
87: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
88: const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0);
89: const PetscScalar dC[] = {2 * C00, 4 * C01, 2 * C11};
91: for (PetscInt k = 0; k < 3; k++) {
92: for (PetscInt j = 0; j < 3; j++) J[k * 3 + j] = alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k];
93: J[k * 3 + k] += alpha * fnorm + u_tShift;
94: }
95: }
97: /* Jacobian for C against C basis functions and gradients of P basis functions */
98: static void JC_0_c0p1(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 J[])
99: {
100: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
101: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
103: J[0] = -c2 * 2 * gradp[0];
104: J[1] = 0.0;
105: J[2] = -c2 * gradp[1];
106: J[3] = -c2 * gradp[0];
107: J[4] = 0.0;
108: J[5] = -c2 * 2 * gradp[1];
109: }
111: /* residual for C when tested against gradients of basis functions */
112: static void C_1(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[])
113: {
114: const PetscReal D = PetscRealPart(constants[D_ID]);
115: for (PetscInt k = 0; k < 3; k++)
116: for (PetscInt d = 0; d < 2; d++) f1[k * 2 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 2 + d];
117: }
119: /* Jacobian for C against gradients of C basis functions */
120: static void JC_1_c1c1(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 J[])
121: {
122: const PetscReal D = PetscRealPart(constants[D_ID]);
123: for (PetscInt k = 0; k < 3; k++)
124: for (PetscInt d = 0; d < 2; d++) J[k * (3 + 1) * 2 * 2 + d * 2 + d] = PetscSqr(D);
125: }
127: /* residual for P when tested against basis functions.
128: The source term always comes from the auxiliary vec because it needs to have zero mean */
129: static void P_0(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[])
130: {
131: PetscScalar S = a[aOff[P_FIELD_ID]];
133: f0[0] = -S;
134: }
136: static inline void QuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal x[2]);
138: /* compute shift to make C positive definite */
139: static inline PetscReal FIX_C(PetscScalar C00, PetscScalar C01, PetscScalar C11)
140: {
141: #if !PetscDefined(USE_COMPLEX)
142: PetscReal eigs[2], s = 0.0;
144: QuadraticRoots(1, -(C00 + C11), C00 * C11 - PetscSqr(C01), eigs);
145: if (eigs[0] < 0 || eigs[1] < 0) s = -PetscMin(eigs[0], eigs[1]) + PETSC_SMALL;
146: return s;
147: #else
148: return 0.0;
149: #endif
150: }
152: /* residual for P when tested against gradients of basis functions */
153: static void P_1(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[])
154: {
155: const PetscReal r = PetscRealPart(constants[R_ID]);
156: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
157: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
158: const PetscScalar C10 = C01;
159: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r;
160: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
161: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
162: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
164: f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1];
165: f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1];
166: }
168: /* Same as above for the P-only subproblem for initial conditions: the conductivity values come from the auxiliary vec */
169: static void P_1_aux(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[])
170: {
171: const PetscReal r = PetscRealPart(constants[R_ID]);
172: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
173: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
174: const PetscScalar C10 = C01;
175: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r;
176: const PetscScalar gradp[] = {u_x[uOff_x[0]], u_x[uOff_x[0] + 1]};
177: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
178: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
180: f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1];
181: f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1];
182: }
184: /* Jacobian for P against gradients of P basis functions */
185: static void JP_1_p1p1(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 J[])
186: {
187: const PetscReal r = PetscRealPart(constants[R_ID]);
188: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
189: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
190: const PetscScalar C10 = C01;
191: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r;
192: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
193: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
195: J[0] = C00 + s;
196: J[1] = C01;
197: J[2] = C10;
198: J[3] = C11 + s;
199: }
201: /* Same as above for the P-only subproblem for initial conditions */
202: static void JP_1_p1p1_aux(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 J[])
203: {
204: const PetscReal r = PetscRealPart(constants[R_ID]);
205: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
206: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
207: const PetscScalar C10 = C01;
208: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r;
209: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
210: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
212: J[0] = C00 + s;
213: J[1] = C01;
214: J[2] = C10;
215: J[3] = C11 + s;
216: }
218: /* Jacobian for P against gradients of P basis functions and C basis functions */
219: static void JP_1_p1c0(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 J[])
220: {
221: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
223: J[0] = gradp[0];
224: J[1] = 0;
225: J[2] = gradp[1];
226: J[3] = gradp[0];
227: J[4] = 0;
228: J[5] = gradp[1];
229: }
231: /* the source term S(x) = exp(-500*||x - x0||^2) */
232: static PetscErrorCode source_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
233: {
234: PetscReal *x0 = (PetscReal *)ctx;
235: PetscReal n = 0;
237: for (PetscInt d = 0; d < dim; ++d) n += (x[d] - x0[d]) * (x[d] - x0[d]);
238: u[0] = PetscExpReal(-500 * n);
239: return PETSC_SUCCESS;
240: }
242: static PetscErrorCode source_1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
243: {
244: PetscScalar ut[1];
245: const PetscReal x0[] = {0.25, 0.25};
246: const PetscReal x1[] = {0.75, 0.75};
248: PetscCall(source_0(dim, time, x, Nf, ut, (void *)x0));
249: PetscCall(source_0(dim, time, x, Nf, u, (void *)x1));
250: u[0] += ut[0];
251: return PETSC_SUCCESS;
252: }
254: /* functionals to be integrated: average -> \int_\Omega u dx */
255: static void average(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 obj[])
256: {
257: obj[0] = u[uOff[P_FIELD_ID]];
258: }
260: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
261: static inline void QuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal x[2])
262: {
263: PetscReal delta = PetscMax(b * b - 4 * a * c, 0); /* eigenvalues symmetric matrix */
264: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(delta));
266: x[0] = temp / a;
267: x[1] = c / temp;
268: }
270: /* functionals to be integrated: energy -> D^2/2 * ||\nabla C||^2 + c^2\nabla p * (r + C) * \nabla p + \alpha/ \gamma * ||C||^\gamma */
271: static void energy(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 obj[])
272: {
273: const PetscReal D = PetscRealPart(constants[D_ID]);
274: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
275: const PetscReal r = PetscRealPart(constants[R_ID]);
276: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
277: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
278: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
279: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
280: const PetscScalar C10 = C01;
281: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
282: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
283: const PetscScalar gradC00[] = {u_x[uOff_x[C_FIELD_ID] + 0], u_x[uOff_x[C_FIELD_ID] + 1]};
284: const PetscScalar gradC01[] = {u_x[uOff_x[C_FIELD_ID] + 2], u_x[uOff_x[C_FIELD_ID] + 3]};
285: const PetscScalar gradC11[] = {u_x[uOff_x[C_FIELD_ID] + 4], u_x[uOff_x[C_FIELD_ID] + 5]};
286: const PetscScalar normC = NORM2C(C00, C01, C11);
287: const PetscScalar normgradC = NORM2C(gradC00[0], gradC01[0], gradC11[0]) + NORM2C(gradC00[1], gradC01[1], gradC11[1]);
288: const PetscScalar nexp = gamma / 2.0;
290: const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC;
291: const PetscScalar t1 = c2 * (gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1]));
292: const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp);
294: obj[0] = t0 + t1 + t2;
295: }
297: /* functionals to be integrated: ellipticity_fail -> 0 means C+r is elliptic at quadrature point, otherwise it returns the absolute value of the most negative eigenvalue */
298: static void ellipticity_fail(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 obj[])
299: {
300: const PetscReal r = PetscRealPart(constants[R_ID]);
301: const PetscReal C00 = PetscRealPart(u[uOff[C_FIELD_ID]] + r);
302: const PetscReal C01 = PetscRealPart(u[uOff[C_FIELD_ID] + 1]);
303: const PetscReal C11 = PetscRealPart(u[uOff[C_FIELD_ID] + 2] + r);
305: PetscReal eigs[2];
306: QuadraticRoots(1, -(C00 + C11), C00 * C11 - PetscSqr(C01), eigs);
307: if (eigs[0] < 0 || eigs[1] < 0) obj[0] = -PetscMin(eigs[0], eigs[1]);
308: else obj[0] = 0.0;
309: }
311: /* initial conditions for C: eq. 16 */
312: static PetscErrorCode initial_conditions_C_0(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
313: {
314: u[0] = 1;
315: u[1] = 0;
316: u[2] = 1;
317: return PETSC_SUCCESS;
318: }
320: /* initial conditions for C: eq. 17 */
321: static PetscErrorCode initial_conditions_C_1(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
322: {
323: const PetscReal x = xx[0];
324: const PetscReal y = xx[1];
326: u[0] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
327: u[1] = 0;
328: u[2] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
329: return PETSC_SUCCESS;
330: }
332: /* initial conditions for C: eq. 18 */
333: static PetscErrorCode initial_conditions_C_2(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
334: {
335: u[0] = 0;
336: u[1] = 0;
337: u[2] = 0;
338: return PETSC_SUCCESS;
339: }
341: /* functionals to be sampled: C * \grad p */
342: static void flux(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 f[])
343: {
344: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
345: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
346: const PetscScalar C10 = C01;
347: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
348: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
350: f[0] = C00 * gradp[0] + C01 * gradp[1];
351: f[1] = C10 * gradp[0] + C11 * gradp[1];
352: }
354: /* functionals to be sampled: ||C|| */
355: static void normc(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 f[])
356: {
357: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
358: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
359: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
361: f[0] = PetscSqrtScalar(NORM2C(C00, C01, C11));
362: }
364: /* functionals to be sampled: zero */
365: static void zero(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 f[])
366: {
367: f[0] = 0.0;
368: }
370: /* functions to be sampled: zero function */
371: static PetscErrorCode zerof(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
372: {
373: for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
374: return PETSC_SUCCESS;
375: }
377: /* functions to be sampled: constant function */
378: static PetscErrorCode constantf(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
379: {
380: PetscInt d;
381: for (d = 0; d < Nc; ++d) u[d] = 1.0;
382: return PETSC_SUCCESS;
383: }
385: /* application context: customizable parameters */
386: typedef struct {
387: PetscReal r;
388: PetscReal eps;
389: PetscReal alpha;
390: PetscReal gamma;
391: PetscReal D;
392: PetscReal c;
393: PetscInt ic_num;
394: PetscInt source_num;
395: PetscReal x0[2];
396: PetscBool lump;
397: PetscBool amr;
398: PetscBool load;
399: char load_filename[PETSC_MAX_PATH_LEN];
400: PetscBool save;
401: char save_filename[PETSC_MAX_PATH_LEN];
402: PetscInt save_every;
403: PetscBool test_restart;
404: PetscBool ellipticity;
405: PetscInt fix_c;
406: } AppCtx;
408: /* process command line options */
409: static PetscErrorCode ProcessOptions(AppCtx *options)
410: {
411: PetscInt dim = PETSC_STATIC_ARRAY_LENGTH(options->x0);
413: PetscFunctionBeginUser;
414: options->r = 1.e-1;
415: options->eps = 1.e-3;
416: options->alpha = 0.75;
417: options->gamma = 0.75;
418: options->c = 5;
419: options->D = 1.e-2;
420: options->ic_num = 0;
421: options->source_num = 0;
422: options->x0[0] = 0.25;
423: options->x0[1] = 0.25;
424: options->lump = PETSC_FALSE;
425: options->amr = PETSC_FALSE;
426: options->load = PETSC_FALSE;
427: options->save = PETSC_FALSE;
428: options->save_every = -1;
429: options->test_restart = PETSC_FALSE;
430: options->ellipticity = PETSC_FALSE;
431: options->fix_c = 1; /* 1 means only Jac, 2 means function and Jac */
433: PetscOptionsBegin(PETSC_COMM_WORLD, "", __FILE__, "DMPLEX");
434: PetscCall(PetscOptionsReal("-alpha", "alpha", __FILE__, options->alpha, &options->alpha, NULL));
435: PetscCall(PetscOptionsReal("-gamma", "gamma", __FILE__, options->gamma, &options->gamma, NULL));
436: PetscCall(PetscOptionsReal("-c", "c", __FILE__, options->c, &options->c, NULL));
437: PetscCall(PetscOptionsReal("-d", "D", __FILE__, options->D, &options->D, NULL));
438: PetscCall(PetscOptionsReal("-eps", "eps", __FILE__, options->eps, &options->eps, NULL));
439: PetscCall(PetscOptionsReal("-r", "r", __FILE__, options->r, &options->r, NULL));
440: PetscCall(PetscOptionsRealArray("-x0", "x0", __FILE__, options->x0, &dim, NULL));
441: PetscCall(PetscOptionsInt("-ic_num", "ic_num", __FILE__, options->ic_num, &options->ic_num, NULL));
442: PetscCall(PetscOptionsInt("-source_num", "source_num", __FILE__, options->source_num, &options->source_num, NULL));
443: PetscCall(PetscOptionsBool("-lump", "use mass lumping", __FILE__, options->lump, &options->lump, NULL));
444: PetscCall(PetscOptionsInt("-fix_c", "shift conductivity to always be positive semi-definite", __FILE__, options->fix_c, &options->fix_c, NULL));
445: PetscCall(PetscOptionsBool("-amr", "use adaptive mesh refinement", __FILE__, options->amr, &options->amr, NULL));
446: PetscCall(PetscOptionsBool("-test_restart", "test restarting files", __FILE__, options->test_restart, &options->test_restart, NULL));
447: if (!options->test_restart) {
448: PetscCall(PetscOptionsString("-load", "filename with data to be loaded for restarting", __FILE__, options->load_filename, options->load_filename, PETSC_MAX_PATH_LEN, &options->load));
449: PetscCall(PetscOptionsString("-save", "filename with data to be saved for restarting", __FILE__, options->save_filename, options->save_filename, PETSC_MAX_PATH_LEN, &options->save));
450: if (options->save) PetscCall(PetscOptionsInt("-save_every", "save every n timestep (-1 saves only the last)", __FILE__, options->save_every, &options->save_every, NULL));
451: }
452: PetscCall(PetscOptionsBool("-monitor_ellipticity", "Dump locations of ellipticity violation", __FILE__, options->ellipticity, &options->ellipticity, NULL));
453: PetscOptionsEnd();
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode SaveToFile(DM dm, Vec u, const char *filename)
458: {
459: #if defined(PETSC_HAVE_HDF5)
460: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
461: PetscViewer viewer;
462: DM cdm = dm;
463: PetscInt numlevels = 0;
465: PetscFunctionBeginUser;
466: while (cdm) {
467: numlevels++;
468: PetscCall(DMGetCoarseDM(cdm, &cdm));
469: }
470: /* Cannot be set programmatically */
471: PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_view_hdf5_storage_version 3.0.0"));
472: PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), filename, FILE_MODE_WRITE, &viewer));
473: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numlevels", PETSC_INT, &numlevels));
474: PetscCall(PetscViewerPushFormat(viewer, format));
475: for (PetscInt level = numlevels - 1; level >= 0; level--) {
476: PetscInt cc, rr;
477: PetscBool isRegular, isUniform;
478: const char *dmname;
479: char groupname[PETSC_MAX_PATH_LEN];
481: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
482: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
483: PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
484: PetscCall(DMGetCoarsenLevel(dm, &cc));
485: PetscCall(DMGetRefineLevel(dm, &rr));
486: PetscCall(DMPlexGetRegularRefinement(dm, &isRegular));
487: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
488: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "meshname", PETSC_STRING, dmname));
489: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refinelevel", PETSC_INT, &rr));
490: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, &cc));
491: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refRegular", PETSC_BOOL, &isRegular));
492: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refUniform", PETSC_BOOL, &isUniform));
493: PetscCall(DMPlexTopologyView(dm, viewer));
494: PetscCall(DMPlexLabelsView(dm, viewer));
495: PetscCall(DMPlexCoordinatesView(dm, viewer));
496: PetscCall(DMPlexSectionView(dm, viewer, NULL));
497: if (level == numlevels - 1) {
498: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
499: PetscCall(DMPlexGlobalVectorView(dm, viewer, NULL, u));
500: }
501: if (level) {
502: PetscInt cStart, cEnd, ccStart, ccEnd, cpStart;
503: DMPolytopeType ct;
504: DMPlexTransform tr;
505: DM sdm;
506: PetscScalar *array;
507: PetscSection section;
508: Vec map;
509: IS gis;
510: const PetscInt *gidx;
512: PetscCall(DMGetCoarseDM(dm, &cdm));
513: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
514: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
515: PetscCall(PetscSectionSetChart(section, cStart, cEnd));
516: for (PetscInt c = cStart; c < cEnd; c++) PetscCall(PetscSectionSetDof(section, c, 1));
517: PetscCall(PetscSectionSetUp(section));
519: PetscCall(DMClone(dm, &sdm));
520: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
521: PetscCall(PetscObjectSetName((PetscObject)section, "pdm_section"));
522: PetscCall(DMSetLocalSection(sdm, section));
523: PetscCall(PetscSectionDestroy(§ion));
525: PetscCall(DMGetLocalVector(sdm, &map));
526: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
527: PetscCall(VecGetArray(map, &array));
528: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
529: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
530: PetscCall(DMPlexTransformSetDM(tr, cdm));
531: PetscCall(DMPlexTransformSetFromOptions(tr));
532: PetscCall(DMPlexTransformSetUp(tr));
533: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
534: PetscCall(DMPlexGetChart(cdm, &cpStart, NULL));
535: PetscCall(DMPlexCreatePointNumbering(cdm, &gis));
536: PetscCall(ISGetIndices(gis, &gidx));
537: for (PetscInt c = ccStart; c < ccEnd; c++) {
538: PetscInt *rsize, *rcone, *rornt, Nt;
539: DMPolytopeType *rct;
540: PetscInt gnum = gidx[c - cpStart] >= 0 ? gidx[c - cpStart] : -(gidx[c - cpStart] + 1);
542: PetscCall(DMPlexGetCellType(cdm, c, &ct));
543: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nt, &rct, &rsize, &rcone, &rornt));
544: for (PetscInt r = 0; r < rsize[Nt - 1]; ++r) {
545: PetscInt pNew;
547: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[Nt - 1], c, r, &pNew));
548: array[pNew - cStart] = gnum;
549: }
550: }
551: PetscCall(ISRestoreIndices(gis, &gidx));
552: PetscCall(ISDestroy(&gis));
553: PetscCall(VecRestoreArray(map, &array));
554: PetscCall(DMPlexTransformDestroy(&tr));
555: PetscCall(DMPlexSectionView(dm, viewer, sdm));
556: PetscCall(DMPlexLocalVectorView(dm, viewer, sdm, map));
557: PetscCall(DMRestoreLocalVector(sdm, &map));
558: PetscCall(DMDestroy(&sdm));
559: }
560: PetscCall(PetscViewerHDF5PopGroup(viewer));
561: PetscCall(DMGetCoarseDM(dm, &dm));
562: }
563: PetscCall(PetscViewerPopFormat(viewer));
564: PetscCall(PetscViewerDestroy(&viewer));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: #else
567: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
568: #endif
569: }
571: static PetscErrorCode LoadFromFile(MPI_Comm comm, const char *filename, DM *odm)
572: {
573: #if defined(PETSC_HAVE_HDF5)
574: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
575: PetscViewer viewer;
576: DM dm, cdm = NULL;
577: PetscSF sfXC = NULL;
578: PetscInt numlevels = -1;
580: PetscFunctionBeginUser;
581: PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer));
582: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numlevels", PETSC_INT, NULL, &numlevels));
583: PetscCall(PetscViewerPushFormat(viewer, format));
584: for (PetscInt level = 0; level < numlevels; level++) {
585: char groupname[PETSC_MAX_PATH_LEN], *dmname;
586: PetscSF sfXB, sfBC, sfG;
587: PetscPartitioner part;
588: PetscInt rr, cc;
589: PetscBool isRegular, isUniform;
591: PetscCall(DMCreate(comm, &dm));
592: PetscCall(DMSetType(dm, DMPLEX));
593: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
594: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
595: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "meshname", PETSC_STRING, NULL, &dmname));
596: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refinelevel", PETSC_INT, NULL, &rr));
597: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, NULL, &cc));
598: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refRegular", PETSC_BOOL, NULL, &isRegular));
599: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refUniform", PETSC_BOOL, NULL, &isUniform));
600: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
601: PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
602: PetscCall(DMPlexLabelsLoad(dm, viewer, sfXB));
603: PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXB));
604: PetscCall(DMPlexGetPartitioner(dm, &part));
605: if (!level) { /* partition the coarse level only */
606: PetscCall(PetscPartitionerSetFromOptions(part));
607: } else { /* propagate partitioning information from coarser to finer level */
608: DM sdm;
609: Vec map;
610: PetscSF sf;
611: PetscLayout clayout;
612: PetscScalar *array;
613: PetscInt *cranks_leaf, *cranks_root, *npoints, *points, *ranks, *starts, *gidxs;
614: PetscInt nparts, cStart, cEnd, nr, ccStart, ccEnd, cpStart, cpEnd;
615: PetscMPIInt size, rank;
617: PetscCall(DMClone(dm, &sdm));
618: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
619: PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXB, NULL, &sf));
620: PetscCall(DMGetLocalVector(sdm, &map));
621: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
622: PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, sf, map));
624: PetscCallMPI(MPI_Comm_size(comm, &size));
625: PetscCallMPI(MPI_Comm_rank(comm, &rank));
626: nparts = size;
627: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
628: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
629: PetscCall(DMPlexGetChart(cdm, &cpStart, &cpEnd));
630: PetscCall(PetscCalloc1(nparts, &npoints));
631: PetscCall(PetscMalloc4(cEnd - cStart, &points, cEnd - cStart, &ranks, nparts + 1, &starts, cEnd - cStart, &gidxs));
632: PetscCall(PetscSFGetGraph(sfXC, &nr, NULL, NULL, NULL));
633: PetscCall(PetscMalloc2(cpEnd - cpStart, &cranks_leaf, nr, &cranks_root));
634: for (PetscInt c = 0; c < cpEnd - cpStart; c++) cranks_leaf[c] = rank;
635: PetscCall(PetscSFReduceBegin(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
636: PetscCall(PetscSFReduceEnd(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
638: PetscCall(VecGetArray(map, &array));
639: for (PetscInt c = 0; c < cEnd - cStart; c++) gidxs[c] = (PetscInt)PetscRealPart(array[c]);
640: PetscCall(VecRestoreArray(map, &array));
642: PetscCall(PetscLayoutCreate(comm, &clayout));
643: PetscCall(PetscLayoutSetLocalSize(clayout, nr));
644: PetscCall(PetscSFSetGraphLayout(sf, clayout, cEnd - cStart, NULL, PETSC_OWN_POINTER, gidxs));
645: PetscCall(PetscLayoutDestroy(&clayout));
647: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
648: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
649: PetscCall(PetscSFDestroy(&sf));
650: PetscCall(PetscFree2(cranks_leaf, cranks_root));
651: for (PetscInt c = 0; c < cEnd - cStart; c++) npoints[ranks[c]]++;
653: starts[0] = 0;
654: for (PetscInt c = 0; c < nparts; c++) starts[c + 1] = starts[c] + npoints[c];
655: for (PetscInt c = 0; c < cEnd - cStart; c++) points[starts[ranks[c]]++] = c;
656: PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
657: PetscCall(PetscPartitionerShellSetPartition(part, nparts, npoints, points));
658: PetscCall(PetscFree(npoints));
659: PetscCall(PetscFree4(points, ranks, starts, gidxs));
660: PetscCall(DMRestoreLocalVector(sdm, &map));
661: PetscCall(DMDestroy(&sdm));
662: }
663: PetscCall(PetscSFDestroy(&sfXC));
664: PetscCall(DMPlexDistribute(dm, 0, &sfBC, odm));
665: if (*odm) {
666: PetscCall(DMDestroy(&dm));
667: dm = *odm;
668: *odm = NULL;
669: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
670: }
671: if (sfBC) PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
672: else {
673: PetscCall(PetscObjectReference((PetscObject)sfXB));
674: sfXC = sfXB;
675: }
676: PetscCall(PetscSFDestroy(&sfXB));
677: PetscCall(PetscSFDestroy(&sfBC));
678: PetscCall(DMSetCoarsenLevel(dm, cc));
679: PetscCall(DMSetRefineLevel(dm, rr));
680: PetscCall(DMPlexSetRegularRefinement(dm, isRegular));
681: PetscCall(DMPlexSetRefinementUniform(dm, isUniform));
682: PetscCall(DMPlexSectionLoad(dm, viewer, NULL, sfXC, &sfG, NULL));
683: if (level == numlevels - 1) {
684: Vec u;
686: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &u));
687: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
688: PetscCall(DMPlexGlobalVectorLoad(dm, viewer, NULL, sfG, u));
689: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &u));
690: }
691: PetscCall(PetscFree(dmname));
692: PetscCall(PetscSFDestroy(&sfG));
693: PetscCall(DMSetCoarseDM(dm, cdm));
694: PetscCall(DMDestroy(&cdm));
695: PetscCall(PetscViewerHDF5PopGroup(viewer));
696: cdm = dm;
697: }
698: *odm = dm;
699: PetscCall(PetscViewerPopFormat(viewer));
700: PetscCall(PetscViewerDestroy(&viewer));
701: PetscCall(PetscSFDestroy(&sfXC));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: #else
704: SETERRQ(comm, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
705: #endif
706: }
708: /* Project source function and make it zero-mean */
709: static PetscErrorCode ProjectSource(DM dm, PetscReal time, AppCtx *ctx)
710: {
711: PetscInt id = C_FIELD_ID;
712: DM dmAux;
713: Vec u, lu;
714: IS is;
715: void *ctxs[NUM_FIELDS];
716: PetscScalar vals[NUM_FIELDS];
717: PetscDS ds;
718: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
720: PetscFunctionBeginUser;
721: switch (ctx->source_num) {
722: case 0:
723: funcs[P_FIELD_ID] = source_0;
724: ctxs[P_FIELD_ID] = ctx->x0;
725: break;
726: case 1:
727: funcs[P_FIELD_ID] = source_1;
728: ctxs[P_FIELD_ID] = NULL;
729: break;
730: default:
731: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown source");
732: }
733: funcs[C_FIELD_ID] = zerof;
734: ctxs[C_FIELD_ID] = NULL;
735: PetscCall(DMGetDS(dm, &ds));
736: PetscCall(DMGetGlobalVector(dm, &u));
737: PetscCall(DMProjectFunction(dm, time, funcs, ctxs, INSERT_ALL_VALUES, u));
738: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
739: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
740: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
741: PetscCall(VecShift(u, -vals[P_FIELD_ID]));
742: PetscCall(DMCreateSubDM(dm, 1, &id, &is, NULL));
743: PetscCall(VecISSet(u, is, 0));
744: PetscCall(ISDestroy(&is));
746: /* Attach source vector as auxiliary vector:
747: Use a different DM to break ref cycles */
748: PetscCall(DMClone(dm, &dmAux));
749: PetscCall(DMCopyDisc(dm, dmAux));
750: PetscCall(DMCreateLocalVector(dmAux, &lu));
751: PetscCall(DMDestroy(&dmAux));
752: PetscCall(DMGlobalToLocal(dm, u, INSERT_VALUES, lu));
753: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, lu));
754: PetscCall(VecViewFromOptions(lu, NULL, "-aux_view"));
755: PetscCall(VecDestroy(&lu));
756: PetscCall(DMRestoreGlobalVector(dm, &u));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /* callback for the creation of the potential null space */
761: static PetscErrorCode CreatePotentialNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
762: {
763: Vec vec;
764: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zerof};
766: PetscFunctionBeginUser;
767: funcs[nfield] = constantf;
768: PetscCall(DMCreateGlobalVector(dm, &vec));
769: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
770: PetscCall(VecNormalize(vec, NULL));
771: PetscCall(PetscObjectSetName((PetscObject)vec, "Potential Null Space"));
772: PetscCall(VecViewFromOptions(vec, NULL, "-potential_nullspace_view"));
773: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
774: /* break ref cycles */
775: PetscCall(VecSetDM(vec, NULL));
776: PetscCall(VecDestroy(&vec));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: static PetscErrorCode DMGetLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
781: {
782: PetscBool has;
784: PetscFunctionBeginUser;
785: if (local) {
786: PetscCall(DMHasNamedLocalVector(dm, "lumped_mass", &has));
787: PetscCall(DMGetNamedLocalVector(dm, "lumped_mass", lumped_mass));
788: } else {
789: PetscCall(DMHasNamedGlobalVector(dm, "lumped_mass", &has));
790: PetscCall(DMGetNamedGlobalVector(dm, "lumped_mass", lumped_mass));
791: }
792: if (!has) {
793: Vec w;
794: IS is;
796: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
797: if (!is) {
798: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
800: PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &is, NULL));
801: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)is));
802: PetscCall(PetscObjectDereference((PetscObject)is));
803: }
804: if (local) {
805: Vec w2, wg;
807: PetscCall(DMCreateMassMatrixLumped(dm, &w, NULL));
808: PetscCall(DMGetGlobalVector(dm, &wg));
809: PetscCall(DMGetLocalVector(dm, &w2));
810: PetscCall(VecSet(w2, 0.0));
811: PetscCall(VecSet(wg, 1.0));
812: PetscCall(VecISSet(wg, is, 0.0));
813: PetscCall(DMGlobalToLocal(dm, wg, INSERT_VALUES, w2));
814: PetscCall(VecPointwiseMult(w, w, w2));
815: PetscCall(DMRestoreGlobalVector(dm, &wg));
816: PetscCall(DMRestoreLocalVector(dm, &w2));
817: } else {
818: PetscCall(DMCreateMassMatrixLumped(dm, NULL, &w));
819: PetscCall(VecISSet(w, is, 0.0));
820: }
821: PetscCall(VecCopy(w, *lumped_mass));
822: PetscCall(VecDestroy(&w));
823: }
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: static PetscErrorCode DMRestoreLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
828: {
829: PetscFunctionBeginUser;
830: if (local) PetscCall(DMRestoreNamedLocalVector(dm, "lumped_mass", lumped_mass));
831: else PetscCall(DMRestoreNamedGlobalVector(dm, "lumped_mass", lumped_mass));
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /* callbacks for lumped mass matrix residual and Jacobian */
836: static PetscErrorCode DMPlexTSComputeIFunctionFEM_Lumped(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
837: {
838: Vec work, local_lumped_mass;
840: PetscFunctionBeginUser;
841: PetscCall(DMGetLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
842: PetscCall(DMGetLocalVector(dm, &work));
843: PetscCall(VecSet(work, 0.0));
844: PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, work, locF, user));
845: PetscCall(VecPointwiseMult(work, locX_t, local_lumped_mass));
846: PetscCall(VecAXPY(locF, 1.0, work));
847: PetscCall(DMRestoreLocalVector(dm, &work));
848: PetscCall(DMRestoreLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: static PetscErrorCode DMPlexTSComputeIJacobianFEM_Lumped(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
853: {
854: Vec lumped_mass, work;
856: PetscFunctionBeginUser;
857: // XXX CHECK DIRK
858: PetscCall(DMGetLumpedMass(dm, PETSC_FALSE, &lumped_mass));
859: PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, 0.0, Jac, JacP, user));
860: PetscCall(DMGetGlobalVector(dm, &work));
861: PetscCall(VecAXPBY(work, X_tShift, 0.0, lumped_mass));
862: PetscCall(MatDiagonalSet(JacP, work, ADD_VALUES));
863: PetscCall(DMRestoreGlobalVector(dm, &work));
864: PetscCall(DMRestoreLumpedMass(dm, PETSC_FALSE, &lumped_mass));
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /* customize residuals and Jacobians */
869: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
870: {
871: PetscDS ds;
872: PetscInt cdim, dim;
873: PetscScalar constants[NUM_CONSTANTS];
875: PetscFunctionBeginUser;
876: constants[R_ID] = ctx->r;
877: constants[EPS_ID] = ctx->eps;
878: constants[ALPHA_ID] = ctx->alpha;
879: constants[GAMMA_ID] = ctx->gamma;
880: constants[D_ID] = ctx->D;
881: constants[C2_ID] = ctx->c * ctx->c;
882: constants[FIXC_ID] = ctx->fix_c;
884: PetscCall(DMGetDimension(dm, &dim));
885: PetscCall(DMGetCoordinateDim(dm, &cdim));
886: PetscCheck(dim == 2 && cdim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2D meshes");
887: PetscCall(DMGetDS(dm, &ds));
888: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
889: PetscCall(PetscDSSetImplicit(ds, C_FIELD_ID, PETSC_TRUE));
890: PetscCall(PetscDSSetImplicit(ds, P_FIELD_ID, PETSC_TRUE));
891: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
892: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
893: PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0, C_1));
894: PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1));
895: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0, NULL, NULL, JC_1_c1c1));
896: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1, NULL, NULL));
897: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0, NULL));
898: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1));
900: /* Attach potential nullspace */
901: PetscCall(DMSetNullSpaceConstructor(dm, P_FIELD_ID, CreatePotentialNullSpace));
903: /* Attach source function as auxiliary vector */
904: PetscCall(ProjectSource(dm, 0, ctx));
906: /* Add callbacks */
907: if (ctx->lump) {
908: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM_Lumped, NULL));
909: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM_Lumped, NULL));
910: } else {
911: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, NULL));
912: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, NULL));
913: }
914: /* This is not really needed because we use Neumann boundaries */
915: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL));
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /* setup discrete spaces and residuals */
920: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
921: {
922: DM plex, cdm = dm;
923: PetscFE feC, feP;
924: PetscBool simplex;
925: PetscInt dim;
926: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
927: MatNullSpace nsp;
929: PetscFunctionBeginUser;
930: PetscCall(DMGetDimension(dm, &dim));
932: PetscCall(DMConvert(dm, DMPLEX, &plex));
933: PetscCall(DMPlexIsSimplex(plex, &simplex));
934: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &simplex, 1, MPIU_BOOL, MPI_LOR, comm));
935: PetscCall(DMDestroy(&plex));
937: /* We model Cij with Cij = Cji -> dim*(dim+1)/2 components */
938: PetscCall(PetscFECreateDefault(comm, dim, (dim * (dim + 1)) / 2, simplex, "c_", -1, &feC));
939: PetscCall(PetscObjectSetName((PetscObject)feC, "conductivity"));
940: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "p_", -1, &feP));
941: PetscCall(PetscObjectSetName((PetscObject)feP, "potential"));
942: PetscCall(MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, &nsp));
943: PetscCall(PetscObjectCompose((PetscObject)feP, "nullspace", (PetscObject)nsp));
944: PetscCall(MatNullSpaceDestroy(&nsp));
945: PetscCall(PetscFECopyQuadrature(feP, feC));
947: PetscCall(DMSetNumFields(dm, 2));
948: PetscCall(DMSetField(dm, C_FIELD_ID, NULL, (PetscObject)feC));
949: PetscCall(DMSetField(dm, P_FIELD_ID, NULL, (PetscObject)feP));
950: PetscCall(PetscFEDestroy(&feC));
951: PetscCall(PetscFEDestroy(&feP));
952: PetscCall(DMCreateDS(dm));
953: while (cdm) {
954: PetscCall(DMCopyDisc(dm, cdm));
955: PetscCall(SetupProblem(cdm, ctx));
956: PetscCall(DMGetCoarseDM(cdm, &cdm));
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: /* Create mesh by command line options */
962: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
963: {
964: PetscFunctionBeginUser;
965: if (ctx->load) {
966: PetscInt refine = 0;
967: PetscBool isHierarchy;
968: DM *dms;
969: char typeName[256];
970: PetscBool flg;
972: PetscCall(LoadFromFile(comm, ctx->load_filename, dm));
973: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
974: PetscCall(PetscOptionsFList("-dm_mat_type", "Matrix type used for created matrices", "DMSetMatType", MatList, MATAIJ, typeName, sizeof(typeName), &flg));
975: if (flg) PetscCall(DMSetMatType(*dm, typeName));
976: PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0));
977: PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0));
978: PetscOptionsEnd();
979: if (refine) {
980: PetscCall(SetupDiscretization(*dm, ctx));
981: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
982: }
983: PetscCall(PetscCalloc1(refine, &dms));
984: if (isHierarchy) PetscCall(DMRefineHierarchy(*dm, refine, dms));
985: for (PetscInt r = 0; r < refine; r++) {
986: Mat M;
987: DM dmr = dms[r];
988: Vec u, ur;
990: if (!isHierarchy) {
991: PetscCall(DMRefine(*dm, PetscObjectComm((PetscObject)*dm), &dmr));
992: PetscCall(DMSetCoarseDM(dmr, *dm));
993: }
994: PetscCall(DMCreateInterpolation(*dm, dmr, &M, NULL));
995: PetscCall(DMGetNamedGlobalVector(*dm, "solution_", &u));
996: PetscCall(DMGetNamedGlobalVector(dmr, "solution_", &ur));
997: PetscCall(MatInterpolate(M, u, ur));
998: PetscCall(DMRestoreNamedGlobalVector(*dm, "solution_", &u));
999: PetscCall(DMRestoreNamedGlobalVector(dmr, "solution_", &ur));
1000: PetscCall(MatDestroy(&M));
1001: if (!isHierarchy) PetscCall(DMSetCoarseDM(dmr, NULL));
1002: PetscCall(DMDestroy(dm));
1003: *dm = dmr;
1004: }
1005: if (refine && !isHierarchy) PetscCall(DMSetRefineLevel(*dm, 0));
1006: PetscCall(PetscFree(dms));
1007: } else {
1008: PetscCall(DMCreate(comm, dm));
1009: PetscCall(DMSetType(*dm, DMPLEX));
1010: PetscCall(DMSetFromOptions(*dm));
1011: PetscCall(DMLocalizeCoordinates(*dm));
1012: {
1013: char convType[256];
1014: PetscBool flg;
1015: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
1016: PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg));
1017: PetscOptionsEnd();
1018: if (flg) {
1019: DM dmConv;
1020: PetscCall(DMConvert(*dm, convType, &dmConv));
1021: if (dmConv) {
1022: PetscCall(DMDestroy(dm));
1023: *dm = dmConv;
1024: PetscCall(DMSetFromOptions(*dm));
1025: PetscCall(DMSetUp(*dm));
1026: }
1027: }
1028: }
1029: }
1030: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }
1034: /* Make potential field zero mean */
1035: static PetscErrorCode ZeroMeanPotential(DM dm, Vec u)
1036: {
1037: PetscScalar vals[NUM_FIELDS];
1038: PetscDS ds;
1039: IS is;
1041: PetscFunctionBeginUser;
1042: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
1043: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
1044: PetscCall(DMGetDS(dm, &ds));
1045: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
1046: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1047: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1048: PetscCall(VecISShift(u, is, -vals[P_FIELD_ID]));
1049: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /* Compute initial conditions and exclude potential from local truncation error
1054: Since we are solving a DAE, once the initial conditions for the differential
1055: variables are set, we need to compute the corresponding value for the
1056: algebraic variables. We do so by creating a subDM for the potential only
1057: and solve a static problem with SNES */
1058: static PetscErrorCode SetInitialConditionsAndTolerances(TS ts, PetscInt nv, Vec vecs[], PetscBool valid)
1059: {
1060: DM dm;
1061: Vec tu, u, p, lsource, subaux, vatol, vrtol;
1062: PetscReal t, atol, rtol;
1063: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
1064: IS isp;
1065: DM dmp;
1066: VecScatter sctp = NULL;
1067: PetscDS ds;
1068: SNES snes;
1069: KSP ksp;
1070: PC pc;
1071: AppCtx *ctx;
1073: PetscFunctionBeginUser;
1074: PetscCall(TSGetDM(ts, &dm));
1075: PetscCall(TSGetApplicationContext(ts, &ctx));
1076: if (valid) {
1077: PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, NULL));
1078: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
1079: PetscCall(DMCreateGlobalVector(dm, &vatol));
1080: PetscCall(DMCreateGlobalVector(dm, &vrtol));
1081: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
1082: PetscCall(VecSet(vatol, atol));
1083: PetscCall(VecISSet(vatol, isp, -1));
1084: PetscCall(VecSet(vrtol, rtol));
1085: PetscCall(VecISSet(vrtol, isp, -1));
1086: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
1087: PetscCall(VecDestroy(&vatol));
1088: PetscCall(VecDestroy(&vrtol));
1089: PetscCall(ISDestroy(&isp));
1090: for (PetscInt i = 0; i < nv; i++) { PetscCall(ZeroMeanPotential(dm, vecs[i])); }
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1093: PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, &dmp));
1094: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
1095: PetscCall(DMSetMatType(dmp, MATAIJ));
1096: PetscCall(DMGetDS(dmp, &ds));
1097: //PetscCall(PetscDSSetResidual(ds, 0, P_0, P_1_aux));
1098: PetscCall(PetscDSSetResidual(ds, 0, 0, P_1_aux));
1099: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux));
1100: PetscCall(DMPlexSetSNESLocalFEM(dmp, PETSC_FALSE, NULL));
1102: PetscCall(DMCreateGlobalVector(dmp, &p));
1104: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dmp), &snes));
1105: PetscCall(SNESSetDM(snes, dmp));
1106: PetscCall(SNESSetOptionsPrefix(snes, "initial_"));
1107: PetscCall(SNESSetErrorIfNotConverged(snes, PETSC_TRUE));
1108: PetscCall(SNESGetKSP(snes, &ksp));
1109: PetscCall(KSPSetType(ksp, KSPFGMRES));
1110: PetscCall(KSPGetPC(ksp, &pc));
1111: PetscCall(PCSetType(pc, PCGAMG));
1112: PetscCall(SNESSetFromOptions(snes));
1113: PetscCall(SNESSetUp(snes));
1115: /* Loop over input vectors and compute corresponding potential */
1116: for (PetscInt i = 0; i < nv; i++) {
1117: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1119: u = vecs[i];
1120: if (!valid) { /* Assumes entries in u are not valid */
1121: PetscCall(TSGetTime(ts, &t));
1122: switch (ctx->ic_num) {
1123: case 0:
1124: funcs[C_FIELD_ID] = initial_conditions_C_0;
1125: break;
1126: case 1:
1127: funcs[C_FIELD_ID] = initial_conditions_C_1;
1128: break;
1129: case 2:
1130: funcs[C_FIELD_ID] = initial_conditions_C_2;
1131: break;
1132: default:
1133: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unknown IC");
1134: }
1135: funcs[P_FIELD_ID] = zerof;
1136: PetscCall(DMProjectFunction(dm, t, funcs, NULL, INSERT_ALL_VALUES, u));
1137: }
1139: /* pass conductivity and source information via auxiliary data */
1140: PetscCall(DMGetGlobalVector(dm, &tu));
1141: PetscCall(VecCopy(u, tu));
1142: PetscCall(VecISSet(tu, isp, 0.0));
1143: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &lsource));
1144: PetscCall(DMCreateLocalVector(dm, &subaux));
1145: PetscCall(DMGlobalToLocal(dm, tu, INSERT_VALUES, subaux));
1146: PetscCall(DMRestoreGlobalVector(dm, &tu));
1147: PetscCall(VecAXPY(subaux, 1.0, lsource));
1148: PetscCall(VecViewFromOptions(subaux, NULL, "-initial_aux_view"));
1149: PetscCall(DMSetAuxiliaryVec(dmp, NULL, 0, 0, subaux));
1150: PetscCall(VecDestroy(&subaux));
1152: /* solve the subproblem */
1153: if (!sctp) PetscCall(VecScatterCreate(u, isp, p, NULL, &sctp));
1154: PetscCall(VecScatterBegin(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1155: PetscCall(VecScatterEnd(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1156: PetscCall(SNESSolve(snes, NULL, p));
1158: /* scatter from potential only to full space */
1159: PetscCall(VecScatterBegin(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1160: PetscCall(VecScatterEnd(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1161: PetscCall(ZeroMeanPotential(dm, u));
1162: }
1163: PetscCall(VecDestroy(&p));
1164: PetscCall(DMDestroy(&dmp));
1165: PetscCall(SNESDestroy(&snes));
1166: PetscCall(VecScatterDestroy(&sctp));
1168: /* exclude potential from computation of the LTE */
1169: PetscCall(DMCreateGlobalVector(dm, &vatol));
1170: PetscCall(DMCreateGlobalVector(dm, &vrtol));
1171: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
1172: PetscCall(VecSet(vatol, atol));
1173: PetscCall(VecISSet(vatol, isp, -1));
1174: PetscCall(VecSet(vrtol, rtol));
1175: PetscCall(VecISSet(vrtol, isp, -1));
1176: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
1177: PetscCall(VecDestroy(&vatol));
1178: PetscCall(VecDestroy(&vrtol));
1179: PetscCall(ISDestroy(&isp));
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1183: /* mesh adaption context */
1184: typedef struct {
1185: VecTagger refineTag;
1186: DMLabel adaptLabel;
1187: PetscInt cnt;
1188: } AdaptCtx;
1190: static PetscErrorCode ResizeSetUp(TS ts, PetscInt nstep, PetscReal time, Vec u, PetscBool *resize, void *vctx)
1191: {
1192: AdaptCtx *ctx = (AdaptCtx *)vctx;
1193: Vec ellVecCells, ellVecCellsF;
1194: DM dm, plex;
1195: PetscDS ds;
1196: PetscReal norm;
1197: PetscInt cStart, cEnd;
1199: PetscFunctionBeginUser;
1200: PetscCall(TSGetDM(ts, &dm));
1201: PetscCall(DMConvert(dm, DMPLEX, &plex));
1202: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
1203: PetscCall(DMDestroy(&plex));
1204: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), NUM_FIELDS * (cEnd - cStart), PETSC_DECIDE, &ellVecCellsF));
1205: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), cEnd - cStart, PETSC_DECIDE, &ellVecCells));
1206: PetscCall(VecSetBlockSize(ellVecCellsF, NUM_FIELDS));
1207: PetscCall(DMGetDS(dm, &ds));
1208: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
1209: PetscCall(DMPlexComputeCellwiseIntegralFEM(dm, u, ellVecCellsF, NULL));
1210: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1211: PetscCall(VecStrideGather(ellVecCellsF, C_FIELD_ID, ellVecCells, INSERT_VALUES));
1212: PetscCall(VecDestroy(&ellVecCellsF));
1213: PetscCall(VecNorm(ellVecCells, NORM_1, &norm));
1214: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "STEP %d norm %g\n", (int)nstep, (double)norm));
1215: if (norm && !ctx->cnt) {
1216: IS refineIS;
1218: *resize = PETSC_TRUE;
1219: if (!ctx->refineTag) {
1220: VecTaggerBox refineBox;
1221: refineBox.min = PETSC_MACHINE_EPSILON;
1222: refineBox.max = PETSC_MAX_REAL;
1224: PetscCall(VecTaggerCreate(PETSC_COMM_SELF, &ctx->refineTag));
1225: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->refineTag, "refine_"));
1226: PetscCall(VecTaggerSetType(ctx->refineTag, VECTAGGERABSOLUTE));
1227: PetscCall(VecTaggerAbsoluteSetBox(ctx->refineTag, &refineBox));
1228: PetscCall(VecTaggerSetFromOptions(ctx->refineTag));
1229: PetscCall(VecTaggerSetUp(ctx->refineTag));
1230: PetscCall(PetscObjectViewFromOptions((PetscObject)ctx->refineTag, NULL, "-tag_view"));
1231: }
1232: PetscCall(DMLabelDestroy(&ctx->adaptLabel));
1233: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)ts), "adapt", &ctx->adaptLabel));
1234: PetscCall(VecTaggerComputeIS(ctx->refineTag, ellVecCells, &refineIS, NULL));
1235: PetscCall(DMLabelSetStratumIS(ctx->adaptLabel, DM_ADAPT_REFINE, refineIS));
1236: PetscCall(ISDestroy(&refineIS));
1237: #if 0
1238: void (*funcs[NUM_FIELDS])(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 f[]);
1239: Vec ellVec;
1241: funcs[P_FIELD_ID] = ellipticity_fail;
1242: funcs[C_FIELD_ID] = NULL;
1244: PetscCall(DMGetGlobalVector(dm, &ellVec));
1245: PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
1246: PetscCall(VecViewFromOptions(ellVec,NULL,"-view_amr_ell"));
1247: PetscCall(DMRestoreGlobalVector(dm, &ellVec));
1248: #endif
1249: ctx->cnt++;
1250: } else {
1251: ctx->cnt = 0;
1252: }
1253: PetscCall(VecDestroy(&ellVecCells));
1254: PetscFunctionReturn(PETSC_SUCCESS);
1255: }
1257: static PetscErrorCode ResizeTransfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *vctx)
1258: {
1259: AdaptCtx *actx = (AdaptCtx *)vctx;
1260: AppCtx *ctx;
1261: DM dm, adm;
1262: PetscReal time;
1264: PetscFunctionBeginUser;
1265: PetscCall(TSGetDM(ts, &dm));
1266: PetscCheck(actx->adaptLabel, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptLabel");
1267: PetscCall(DMAdaptLabel(dm, actx->adaptLabel, &adm));
1268: PetscCheck(adm, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adapted DM");
1269: PetscCall(TSGetTime(ts, &time));
1270: PetscCall(DMLabelDestroy(&actx->adaptLabel));
1271: for (PetscInt i = 0; i < nv; i++) {
1272: PetscCall(DMCreateGlobalVector(adm, &vecsout[i]));
1273: PetscCall(DMForestTransferVec(dm, vecsin[i], adm, vecsout[i], PETSC_TRUE, time));
1274: }
1275: PetscCall(DMForestSetAdaptivityForest(adm, NULL));
1276: PetscCall(DMSetCoarseDM(adm, NULL));
1277: PetscCall(DMSetLocalSection(adm, NULL));
1278: PetscCall(TSSetDM(ts, adm));
1279: PetscCall(TSGetTime(ts, &time));
1280: PetscCall(TSGetApplicationContext(ts, &ctx));
1281: PetscCall(DMSetNullSpaceConstructor(adm, P_FIELD_ID, CreatePotentialNullSpace));
1282: PetscCall(ProjectSource(adm, time, ctx));
1283: PetscCall(SetInitialConditionsAndTolerances(ts, nv, vecsout, PETSC_TRUE));
1284: PetscCall(DMDestroy(&adm));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1289: {
1290: PetscFunctionBeginUser;
1291: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1292: PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
1293: PetscCall(PetscViewerFileSetName(*viewer, filename));
1294: PetscFunctionReturn(PETSC_SUCCESS);
1295: }
1297: /* Monitor relevant functionals */
1298: static PetscErrorCode Monitor(TS ts, PetscInt stepnum, PetscReal time, Vec u, void *vctx)
1299: {
1300: PetscScalar vals[2 * NUM_FIELDS];
1301: DM dm;
1302: PetscDS ds;
1303: AppCtx *ctx;
1305: PetscFunctionBeginUser;
1306: PetscCall(TSGetDM(ts, &dm));
1307: PetscCall(TSGetApplicationContext(ts, &ctx));
1308: PetscCall(DMGetDS(dm, &ds));
1310: /* monitor energy and potential average */
1311: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
1312: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1313: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1315: /* monitor ellipticity_fail */
1316: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
1317: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals + NUM_FIELDS, NULL));
1318: if (ctx->ellipticity) {
1319: void (*funcs[NUM_FIELDS])(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 f[]);
1320: Vec ellVec;
1321: PetscViewer viewer;
1322: char filename[PETSC_MAX_PATH_LEN];
1324: funcs[P_FIELD_ID] = ellipticity_fail;
1325: funcs[C_FIELD_ID] = NULL;
1327: PetscCall(DMGetGlobalVector(dm, &ellVec));
1328: PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
1329: PetscCall(PetscSNPrintf(filename, sizeof filename, "ellipticity_fail-%03" PetscInt_FMT ".vtu", stepnum));
1330: PetscCall(OutputVTK(dm, filename, &viewer));
1331: PetscCall(VecView(ellVec, viewer));
1332: PetscCall(PetscViewerDestroy(&viewer));
1333: PetscCall(DMRestoreGlobalVector(dm, &ellVec));
1334: }
1335: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1337: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%4" PetscInt_FMT " TS: time %g, energy %g, intp %g, ell %g\n", stepnum, (double)time, (double)PetscRealPart(vals[C_FIELD_ID]), (double)PetscRealPart(vals[P_FIELD_ID]), (double)PetscRealPart(vals[NUM_FIELDS + C_FIELD_ID])));
1338: PetscFunctionReturn(PETSC_SUCCESS);
1339: }
1341: /* Save restart information */
1342: static PetscErrorCode MonitorSave(TS ts, PetscInt steps, PetscReal time, Vec u, void *vctx)
1343: {
1344: DM dm;
1345: AppCtx *ctx = (AppCtx *)vctx;
1346: PetscInt save_every = ctx->save_every;
1347: TSConvergedReason reason;
1349: PetscFunctionBeginUser;
1350: if (!ctx->save) PetscFunctionReturn(PETSC_SUCCESS);
1351: PetscCall(TSGetDM(ts, &dm));
1352: PetscCall(TSGetConvergedReason(ts, &reason));
1353: if ((save_every > 0 && steps % save_every == 0) || (save_every == -1 && reason) || save_every < -1) PetscCall(SaveToFile(dm, u, ctx->save_filename));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /* Make potential zero mean after SNES solve */
1358: static PetscErrorCode PostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
1359: {
1360: DM dm;
1361: Vec u = Y[stageindex];
1362: SNES snes;
1363: PetscInt nits, lits, stepnum;
1364: AppCtx *ctx;
1366: PetscFunctionBeginUser;
1367: PetscCall(TSGetDM(ts, &dm));
1368: PetscCall(ZeroMeanPotential(dm, u));
1370: PetscCall(TSGetApplicationContext(ts, &ctx));
1371: if (ctx->test_restart) PetscFunctionReturn(PETSC_SUCCESS);
1373: /* monitor linear and nonlinear iterations */
1374: PetscCall(TSGetStepNumber(ts, &stepnum));
1375: PetscCall(TSGetSNES(ts, &snes));
1376: PetscCall(SNESGetIterationNumber(snes, &nits));
1377: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1379: /* if function evals in TSDIRK are zero in the first stage, it is FSAL */
1380: if (stageindex == 0) {
1381: PetscBool dirk;
1382: PetscInt nf;
1384: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
1385: PetscCall(SNESGetNumberFunctionEvals(snes, &nf));
1386: if (dirk && nf == 0) nits = 0;
1387: }
1388: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), " step %" PetscInt_FMT " stage %" PetscInt_FMT " nonlinear its %" PetscInt_FMT ", linear its %" PetscInt_FMT "\n", stepnum, stageindex, nits, lits));
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: static PetscErrorCode VecViewFlux(Vec u, const char *opts)
1393: {
1394: Vec fluxVec;
1395: DM dmFlux, dm, plex;
1396: PetscInt dim;
1397: PetscFE feC, feFluxC, feNormC;
1398: PetscBool simplex, has;
1400: void (*funcs[])(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 f[]) = {normc, flux};
1402: PetscFunctionBeginUser;
1403: PetscCall(PetscOptionsHasName(NULL, NULL, opts, &has));
1404: if (!has) PetscFunctionReturn(PETSC_SUCCESS);
1405: PetscCall(VecGetDM(u, &dm));
1406: PetscCall(DMGetDimension(dm, &dim));
1407: PetscCall(DMGetField(dm, C_FIELD_ID, NULL, (PetscObject *)&feC));
1408: PetscCall(DMConvert(dm, DMPLEX, &plex));
1409: PetscCall(DMPlexIsSimplex(plex, &simplex));
1410: PetscCall(DMDestroy(&plex));
1411: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, simplex, "flux_", -1, &feFluxC));
1412: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, "normc_", -1, &feNormC));
1413: PetscCall(PetscFECopyQuadrature(feC, feFluxC));
1414: PetscCall(PetscFECopyQuadrature(feC, feNormC));
1415: PetscCall(DMClone(dm, &dmFlux));
1416: PetscCall(DMSetNumFields(dmFlux, 1));
1417: PetscCall(DMSetField(dmFlux, 0, NULL, (PetscObject)feNormC));
1418: /* paraview segfaults! */
1419: //PetscCall(DMSetField(dmFlux, 1, NULL, (PetscObject)feFluxC));
1420: PetscCall(DMCreateDS(dmFlux));
1421: PetscCall(PetscFEDestroy(&feFluxC));
1422: PetscCall(PetscFEDestroy(&feNormC));
1424: PetscCall(DMGetGlobalVector(dmFlux, &fluxVec));
1425: PetscCall(DMProjectField(dmFlux, 0.0, u, funcs, INSERT_VALUES, fluxVec));
1426: PetscCall(VecViewFromOptions(fluxVec, NULL, opts));
1427: PetscCall(DMRestoreGlobalVector(dmFlux, &fluxVec));
1428: PetscCall(DMDestroy(&dmFlux));
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: static PetscErrorCode Run(MPI_Comm comm, AppCtx *ctx)
1433: {
1434: DM dm;
1435: TS ts;
1436: Vec u;
1437: AdaptCtx *actx;
1438: PetscBool flg;
1440: PetscFunctionBeginUser;
1441: if (ctx->test_restart) {
1442: PetscViewer subviewer;
1443: PetscMPIInt rank;
1445: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1446: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1447: if (ctx->load) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d loading from %s\n", rank, ctx->load_filename));
1448: if (ctx->save) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d saving to %s\n", rank, ctx->save_filename));
1449: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1450: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1451: } else {
1452: PetscCall(PetscPrintf(comm, "----------------------------\n"));
1453: PetscCall(PetscPrintf(comm, "Simulation parameters:\n"));
1454: PetscCall(PetscPrintf(comm, " r : %g\n", (double)ctx->r));
1455: PetscCall(PetscPrintf(comm, " eps : %g\n", (double)ctx->eps));
1456: PetscCall(PetscPrintf(comm, " alpha: %g\n", (double)ctx->alpha));
1457: PetscCall(PetscPrintf(comm, " gamma: %g\n", (double)ctx->gamma));
1458: PetscCall(PetscPrintf(comm, " D : %g\n", (double)ctx->D));
1459: PetscCall(PetscPrintf(comm, " c : %g\n", (double)ctx->c));
1460: if (ctx->load) PetscCall(PetscPrintf(comm, " load : %s\n", ctx->load_filename));
1461: else PetscCall(PetscPrintf(comm, " IC : %" PetscInt_FMT "\n", ctx->ic_num));
1462: PetscCall(PetscPrintf(comm, " S : %" PetscInt_FMT "\n", ctx->source_num));
1463: PetscCall(PetscPrintf(comm, " x0 : (%g,%g)\n", (double)ctx->x0[0], (double)ctx->x0[1]));
1464: PetscCall(PetscPrintf(comm, "----------------------------\n"));
1465: }
1467: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SetupStage));
1468: PetscCall(CreateMesh(comm, &dm, ctx));
1469: PetscCall(SetupDiscretization(dm, ctx));
1471: PetscCall(TSCreate(comm, &ts));
1472: PetscCall(TSSetApplicationContext(ts, ctx));
1474: PetscCall(TSSetDM(ts, dm));
1475: if (ctx->test_restart) {
1476: PetscViewer subviewer;
1477: PetscMPIInt rank;
1478: PetscInt level;
1480: PetscCall(DMGetRefineLevel(dm, &level));
1481: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1482: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1483: PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d DM refinement level %" PetscInt_FMT "\n", rank, level));
1484: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1485: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1486: }
1488: if (ctx->test_restart) PetscCall(TSSetMaxSteps(ts, 1));
1489: PetscCall(TSSetMaxTime(ts, 10.0));
1490: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1491: if (!ctx->test_restart) PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
1492: PetscCall(TSMonitorSet(ts, MonitorSave, ctx, NULL));
1493: PetscCall(PetscNew(&actx));
1494: if (ctx->amr) PetscCall(TSSetResize(ts, PETSC_TRUE, ResizeSetUp, ResizeTransfer, actx));
1495: PetscCall(TSSetPostStage(ts, PostStage));
1496: PetscCall(TSSetMaxSNESFailures(ts, -1));
1497: PetscCall(TSSetFromOptions(ts));
1499: PetscCall(DMCreateGlobalVector(dm, &u));
1500: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
1501: PetscCall(DMHasNamedGlobalVector(dm, "solution_", &flg));
1502: if (flg) {
1503: Vec ru;
1505: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &ru));
1506: PetscCall(VecCopy(ru, u));
1507: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &ru));
1508: }
1509: PetscCall(SetInitialConditionsAndTolerances(ts, 1, &u, PETSC_FALSE));
1510: PetscCall(TSSetSolution(ts, u));
1511: PetscCall(VecDestroy(&u));
1512: PetscCall(DMDestroy(&dm));
1513: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
1515: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SolveStage));
1516: PetscCall(TSSolve(ts, NULL));
1517: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
1519: PetscCall(TSGetSolution(ts, &u));
1520: PetscCall(VecViewFromOptions(u, NULL, "-final_view"));
1521: PetscCall(VecViewFlux(u, "-final_flux_view"));
1523: PetscCall(TSDestroy(&ts));
1524: PetscCall(VecTaggerDestroy(&actx->refineTag));
1525: PetscCall(PetscFree(actx));
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: int main(int argc, char **argv)
1530: {
1531: AppCtx ctx;
1533: PetscFunctionBeginUser;
1534: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1535: PetscCall(ProcessOptions(&ctx));
1536: PetscCall(PetscLogStageRegister("Setup", &SetupStage));
1537: PetscCall(PetscLogStageRegister("Solve", &SolveStage));
1538: if (ctx.test_restart) { /* Test sequences of save and loads */
1539: PetscMPIInt rank;
1541: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1543: /* test saving */
1544: ctx.load = PETSC_FALSE;
1545: ctx.save = PETSC_TRUE;
1546: /* sequential save */
1547: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential save\n"));
1548: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_seq_%d.h5", rank));
1549: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1550: /* parallel save */
1551: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel save\n"));
1552: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_par.h5"));
1553: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1555: /* test loading */
1556: ctx.load = PETSC_TRUE;
1557: ctx.save = PETSC_FALSE;
1558: /* sequential and parallel runs from sequential save */
1559: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_seq_0.h5"));
1560: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from sequential save\n"));
1561: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1562: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from sequential save\n"));
1563: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1564: /* sequential and parallel runs from parallel save */
1565: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_par.h5"));
1566: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from parallel save\n"));
1567: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1568: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from parallel save\n"));
1569: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1570: } else { /* Run the simulation */
1571: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1572: }
1573: PetscCall(PetscFinalize());
1574: return 0;
1575: }
1577: /*TEST
1579: testset:
1580: args: -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 1 -initial_snes_test_jacobian -snes_test_jacobian -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none
1582: test:
1583: suffix: 0
1584: nsize: {{1 2}}
1585: args: -dm_refine 1 -lump {{0 1}}
1587: test:
1588: suffix: 0_dirk
1589: nsize: {{1 2}}
1590: args: -dm_refine 1 -ts_type dirk
1592: test:
1593: suffix: 0_dirk_mg
1594: nsize: {{1 2}}
1595: args: -dm_refine_hierarchy 1 -ts_type dirk -pc_type mg -mg_levels_pc_type bjacobi -mg_levels_sub_pc_factor_levels 2 -mg_levels_sub_pc_factor_mat_ordering_type rcm -mg_levels_sub_pc_factor_reuse_ordering -mg_coarse_pc_type svd -lump {{0 1}}
1597: test:
1598: suffix: 0_dirk_fieldsplit
1599: nsize: {{1 2}}
1600: args: -dm_refine 1 -ts_type dirk -pc_type fieldsplit -pc_fieldsplit_type multiplicative -lump {{0 1}}
1602: test:
1603: requires: p4est
1604: suffix: 0_p4est
1605: nsize: {{1 2}}
1606: args: -dm_refine 1 -dm_plex_convert_type p4est -lump 0
1608: test:
1609: suffix: 0_periodic
1610: nsize: {{1 2}}
1611: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -lump {{0 1}}
1613: test:
1614: requires: p4est
1615: suffix: 0_p4est_periodic
1616: nsize: {{1 2}}
1617: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -dm_plex_convert_type p4est -lump 0
1619: test:
1620: requires: p4est
1621: suffix: 0_p4est_mg
1622: nsize: {{1 2}}
1623: args: -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_plex_convert_type p4est -pc_type mg -mg_coarse_pc_type svd -mg_levels_pc_type svd -lump 0
1625: testset:
1626: requires: hdf5
1627: args: -test_restart -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type mg -mg_levels_pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -petscpartitioner_type simple -test_restart
1629: test:
1630: requires: !single
1631: suffix: restart
1632: nsize: {{1 2}separate output}
1633: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 0
1635: test:
1636: requires: triangle
1637: suffix: restart_simplex
1638: nsize: {{1 2}separate output}
1639: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 1
1641: test:
1642: requires: !single
1643: suffix: restart_refonly
1644: nsize: {{1 2}separate output}
1645: args: -dm_refine 1 -dm_plex_simplex 0
1647: test:
1648: requires: triangle
1649: suffix: restart_simplex_refonly
1650: nsize: {{1 2}separate output}
1651: args: -dm_refine 1 -dm_plex_simplex 1
1653: TEST*/