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: NUM_CONSTANTS
50: } ConstantIdx;
52: PetscLogStage SetupStage, SolveStage;
54: #define NORM2C(c00, c01, c11) PetscSqr(c00) + 2 * PetscSqr(c01) + PetscSqr(c11)
56: /* residual for C when tested against basis functions */
57: 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[])
58: {
59: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
60: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
61: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
62: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
63: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
64: const PetscScalar crossp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[1] * gradp[1]};
65: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
66: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
67: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
68: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
69: const PetscScalar nexp = (gamma - 2.0) / 2.0;
70: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
72: 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];
73: }
75: /* Jacobian for C against C basis functions */
76: 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[])
77: {
78: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
79: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
80: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
81: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
82: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
83: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
84: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
85: const PetscScalar nexp = (gamma - 2.0) / 2.0;
86: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
87: const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0);
88: const PetscScalar dC[] = {2 * C00, 4 * C01, 2 * C11};
90: for (PetscInt k = 0; k < 3; k++) {
91: for (PetscInt j = 0; j < 3; j++) J[k * 3 + j] = alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k];
92: J[k * 3 + k] += alpha * fnorm + u_tShift;
93: }
94: }
96: /* Jacobian for C against C basis functions and gradients of P basis functions */
97: 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[])
98: {
99: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
100: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
102: J[0] = -c2 * 2 * gradp[0];
103: J[1] = 0.0;
104: J[2] = -c2 * gradp[1];
105: J[3] = -c2 * gradp[0];
106: J[4] = 0.0;
107: J[5] = -c2 * 2 * gradp[1];
108: }
110: /* residual for C when tested against gradients of basis functions */
111: 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[])
112: {
113: const PetscReal D = PetscRealPart(constants[D_ID]);
114: for (PetscInt k = 0; k < 3; k++)
115: for (PetscInt d = 0; d < 2; d++) f1[k * 2 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 2 + d];
116: }
118: /* Jacobian for C against gradients of C basis functions */
119: 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[])
120: {
121: const PetscReal D = PetscRealPart(constants[D_ID]);
122: for (PetscInt k = 0; k < 3; k++)
123: for (PetscInt d = 0; d < 2; d++) J[k * (3 + 1) * 2 * 2 + d * 2 + d] = PetscSqr(D);
124: }
126: /* residual for P when tested against basis functions.
127: The source term always comes from the auxiliary vec because it needs to have zero mean */
128: 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[])
129: {
130: PetscScalar S = a[aOff[P_FIELD_ID]];
132: f0[0] = -S;
133: }
135: /* residual for P when tested against gradients of basis functions */
136: 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[])
137: {
138: const PetscReal r = PetscRealPart(constants[R_ID]);
139: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
140: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
141: const PetscScalar C10 = C01;
142: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
143: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
145: f1[0] = (C00 + r) * gradp[0] + C01 * gradp[1];
146: f1[1] = C10 * gradp[0] + (C11 + r) * gradp[1];
147: }
149: /* Same as above for the P-only subproblem for initial conditions: the conductivity values come from the auxiliary vec */
150: 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[])
151: {
152: const PetscReal r = PetscRealPart(constants[R_ID]);
153: const PetscScalar C00 = a[aOff[C_FIELD_ID]];
154: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
155: const PetscScalar C10 = C01;
156: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2];
157: const PetscScalar gradp[] = {u_x[uOff_x[0]], u_x[uOff_x[0] + 1]};
159: f1[0] = (C00 + r) * gradp[0] + C01 * gradp[1];
160: f1[1] = C10 * gradp[0] + (C11 + r) * gradp[1];
161: }
163: /* Jacobian for P against gradients of P basis functions */
164: 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[])
165: {
166: const PetscReal r = PetscRealPart(constants[R_ID]);
167: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
168: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
169: const PetscScalar C10 = C01;
170: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
172: J[0] = C00 + r;
173: J[1] = C01;
174: J[2] = C10;
175: J[3] = C11 + r;
176: }
178: /* Same as above for the P-only subproblem for initial conditions */
179: 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[])
180: {
181: const PetscReal r = PetscRealPart(constants[R_ID]);
182: const PetscScalar C00 = a[aOff[C_FIELD_ID]];
183: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
184: const PetscScalar C10 = C01;
185: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2];
187: J[0] = C00 + r;
188: J[1] = C01;
189: J[2] = C10;
190: J[3] = C11 + r;
191: }
193: /* Jacobian for P against gradients of P basis functions and C basis functions */
194: 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[])
195: {
196: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
198: J[0] = gradp[0];
199: J[1] = 0;
200: J[2] = gradp[1];
201: J[3] = gradp[0];
202: J[4] = 0;
203: J[5] = gradp[1];
204: }
206: /* the source term S(x) = exp(-500*||x - x0||^2) */
207: static PetscErrorCode source_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
208: {
209: PetscReal *x0 = (PetscReal *)ctx;
210: PetscReal n = 0;
212: for (PetscInt d = 0; d < dim; ++d) n += (x[d] - x0[d]) * (x[d] - x0[d]);
213: u[0] = PetscExpReal(-500 * n);
214: return PETSC_SUCCESS;
215: }
217: static PetscErrorCode source_1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
218: {
219: PetscScalar ut[1];
220: const PetscReal x0[] = {0.25, 0.25};
221: const PetscReal x1[] = {0.75, 0.75};
223: PetscCall(source_0(dim, time, x, Nf, ut, (void *)x0));
224: PetscCall(source_0(dim, time, x, Nf, u, (void *)x1));
225: u[0] += ut[0];
226: return PETSC_SUCCESS;
227: }
229: /* functionals to be integrated: average -> \int_\Omega u dx */
230: 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[])
231: {
232: obj[0] = u[uOff[P_FIELD_ID]];
233: }
235: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
236: static inline void QuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal x[2])
237: {
238: PetscReal delta = PetscMax(b * b - 4 * a * c, 0); /* eigenvalues symmetric matrix */
239: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(delta));
241: x[0] = temp / a;
242: x[1] = c / temp;
243: }
245: /* functionals to be integrated: energy -> D^2/2 * ||\nabla C||^2 + c^2\nabla p * (r + C) * \nabla p + \alpha/ \gamma * ||C||^\gamma */
246: 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[])
247: {
248: const PetscReal D = PetscRealPart(constants[D_ID]);
249: const PetscReal c2 = PetscRealPart(constants[C2_ID]);
250: const PetscReal r = PetscRealPart(constants[R_ID]);
251: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
252: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
253: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
254: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
255: const PetscScalar C10 = C01;
256: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
257: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
258: const PetscScalar gradC00[] = {u_x[uOff_x[C_FIELD_ID] + 0], u_x[uOff_x[C_FIELD_ID] + 1]};
259: const PetscScalar gradC01[] = {u_x[uOff_x[C_FIELD_ID] + 2], u_x[uOff_x[C_FIELD_ID] + 3]};
260: const PetscScalar gradC11[] = {u_x[uOff_x[C_FIELD_ID] + 4], u_x[uOff_x[C_FIELD_ID] + 5]};
261: const PetscScalar normC = NORM2C(C00, C01, C11);
262: const PetscScalar normgradC = NORM2C(gradC00[0], gradC01[0], gradC11[0]) + NORM2C(gradC00[1], gradC01[1], gradC11[1]);
263: const PetscScalar nexp = gamma / 2.0;
265: const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC;
266: const PetscScalar t1 = c2 * (gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1]));
267: const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp);
269: obj[0] = t0 + t1 + t2;
270: }
272: /* 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 */
273: 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[])
274: {
275: const PetscReal r = PetscRealPart(constants[R_ID]);
276: const PetscReal C00 = PetscRealPart(u[uOff[C_FIELD_ID]] + r);
277: const PetscReal C01 = PetscRealPart(u[uOff[C_FIELD_ID] + 1]);
278: const PetscReal C11 = PetscRealPart(u[uOff[C_FIELD_ID] + 2] + r);
280: PetscReal eigs[2];
281: QuadraticRoots(1, -(C00 + C11), C00 * C11 - PetscSqr(C01), eigs);
282: if (eigs[0] < 0 || eigs[1] < 0) obj[0] = -PetscMin(eigs[0], eigs[1]);
283: else obj[0] = 0.0;
284: }
286: /* initial conditions for C: eq. 16 */
287: static PetscErrorCode initial_conditions_C_0(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
288: {
289: u[0] = 1;
290: u[1] = 0;
291: u[2] = 1;
292: return PETSC_SUCCESS;
293: }
295: /* initial conditions for C: eq. 17 */
296: static PetscErrorCode initial_conditions_C_1(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
297: {
298: const PetscReal x = xx[0];
299: const PetscReal y = xx[1];
301: u[0] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
302: u[1] = 0;
303: u[2] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
304: return PETSC_SUCCESS;
305: }
307: /* initial conditions for C: eq. 18 */
308: static PetscErrorCode initial_conditions_C_2(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
309: {
310: u[0] = 0;
311: u[1] = 0;
312: u[2] = 0;
313: return PETSC_SUCCESS;
314: }
316: /* functionals to be sampled: C * \grad p */
317: 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[])
318: {
319: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
320: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
321: const PetscScalar C10 = C01;
322: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
323: const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
325: f[0] = C00 * gradp[0] + C01 * gradp[1];
326: f[1] = C10 * gradp[0] + C11 * gradp[1];
327: }
329: /* functionals to be sampled: ||C|| */
330: 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[])
331: {
332: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
333: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
334: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
336: f[0] = PetscSqrtScalar(NORM2C(C00, C01, C11));
337: }
339: /* functionals to be sampled: zero */
340: 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[])
341: {
342: f[0] = 0.0;
343: }
345: /* functions to be sampled: zero function */
346: static PetscErrorCode zerof(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
347: {
348: for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
349: return PETSC_SUCCESS;
350: }
352: /* functions to be sampled: constant function */
353: static PetscErrorCode constantf(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
354: {
355: PetscInt d;
356: for (d = 0; d < Nc; ++d) u[d] = 1.0;
357: return PETSC_SUCCESS;
358: }
360: /* application context: customizable parameters */
361: typedef struct {
362: PetscReal r;
363: PetscReal eps;
364: PetscReal alpha;
365: PetscReal gamma;
366: PetscReal D;
367: PetscReal c;
368: PetscInt ic_num;
369: PetscInt source_num;
370: PetscReal x0[2];
371: PetscBool amr;
372: PetscBool load;
373: char load_filename[PETSC_MAX_PATH_LEN];
374: PetscBool save;
375: char save_filename[PETSC_MAX_PATH_LEN];
376: PetscInt save_every;
377: PetscBool test_restart;
378: PetscBool ellipticity;
379: } AppCtx;
381: /* process command line options */
382: static PetscErrorCode ProcessOptions(AppCtx *options)
383: {
384: PetscInt dim = PETSC_STATIC_ARRAY_LENGTH(options->x0);
386: PetscFunctionBeginUser;
387: options->r = 1.e-1;
388: options->eps = 1.e-3;
389: options->alpha = 0.75;
390: options->gamma = 0.75;
391: options->c = 5;
392: options->D = 1.e-2;
393: options->ic_num = 0;
394: options->source_num = 0;
395: options->x0[0] = 0.25;
396: options->x0[1] = 0.25;
397: options->amr = PETSC_FALSE;
398: options->load = PETSC_FALSE;
399: options->save = PETSC_FALSE;
400: options->save_every = -1;
401: options->test_restart = PETSC_FALSE;
402: options->ellipticity = PETSC_FALSE;
404: PetscOptionsBegin(PETSC_COMM_WORLD, "", __FILE__, "DMPLEX");
405: PetscCall(PetscOptionsReal("-alpha", "alpha", __FILE__, options->alpha, &options->alpha, NULL));
406: PetscCall(PetscOptionsReal("-gamma", "gamma", __FILE__, options->gamma, &options->gamma, NULL));
407: PetscCall(PetscOptionsReal("-c", "c", __FILE__, options->c, &options->c, NULL));
408: PetscCall(PetscOptionsReal("-d", "D", __FILE__, options->D, &options->D, NULL));
409: PetscCall(PetscOptionsReal("-eps", "eps", __FILE__, options->eps, &options->eps, NULL));
410: PetscCall(PetscOptionsReal("-r", "r", __FILE__, options->r, &options->r, NULL));
411: PetscCall(PetscOptionsRealArray("-x0", "x0", __FILE__, options->x0, &dim, NULL));
412: PetscCall(PetscOptionsInt("-ic_num", "ic_num", __FILE__, options->ic_num, &options->ic_num, NULL));
413: PetscCall(PetscOptionsInt("-source_num", "source_num", __FILE__, options->source_num, &options->source_num, NULL));
414: PetscCall(PetscOptionsBool("-amr", "use adaptive mesh refinement", __FILE__, options->amr, &options->amr, NULL));
415: PetscCall(PetscOptionsBool("-test_restart", "test restarting files", __FILE__, options->test_restart, &options->test_restart, NULL));
416: if (!options->test_restart) {
417: PetscCall(PetscOptionsString("-load", "filename with data to be loaded for restarting", __FILE__, options->load_filename, options->load_filename, PETSC_MAX_PATH_LEN, &options->load));
418: PetscCall(PetscOptionsString("-save", "filename with data to be saved for restarting", __FILE__, options->save_filename, options->save_filename, PETSC_MAX_PATH_LEN, &options->save));
419: if (options->save) PetscCall(PetscOptionsInt("-save_every", "save every n timestep (-1 saves only the last)", __FILE__, options->save_every, &options->save_every, NULL));
420: }
421: PetscCall(PetscOptionsBool("-monitor_ellipticity", "Dump locations of ellipticity violation", __FILE__, options->ellipticity, &options->ellipticity, NULL));
422: PetscOptionsEnd();
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode SaveToFile(DM dm, Vec u, const char *filename)
427: {
428: #if defined(PETSC_HAVE_HDF5)
429: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
430: PetscViewer viewer;
431: DM cdm = dm;
432: PetscInt numlevels = 0;
434: PetscFunctionBeginUser;
435: while (cdm) {
436: numlevels++;
437: PetscCall(DMGetCoarseDM(cdm, &cdm));
438: }
439: /* Cannot be set programmatically */
440: PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_view_hdf5_storage_version 3.0.0"));
441: PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), filename, FILE_MODE_WRITE, &viewer));
442: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numlevels", PETSC_INT, &numlevels));
443: PetscCall(PetscViewerPushFormat(viewer, format));
444: for (PetscInt level = numlevels - 1; level >= 0; level--) {
445: PetscInt cc, rr;
446: PetscBool isRegular, isUniform;
447: const char *dmname;
448: char groupname[PETSC_MAX_PATH_LEN];
450: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
451: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
452: PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
453: PetscCall(DMGetCoarsenLevel(dm, &cc));
454: PetscCall(DMGetRefineLevel(dm, &rr));
455: PetscCall(DMPlexGetRegularRefinement(dm, &isRegular));
456: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
457: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "meshname", PETSC_STRING, dmname));
458: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refinelevel", PETSC_INT, &rr));
459: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, &cc));
460: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refRegular", PETSC_BOOL, &isRegular));
461: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refUniform", PETSC_BOOL, &isUniform));
462: PetscCall(DMPlexTopologyView(dm, viewer));
463: PetscCall(DMPlexLabelsView(dm, viewer));
464: PetscCall(DMPlexCoordinatesView(dm, viewer));
465: PetscCall(DMPlexSectionView(dm, viewer, NULL));
466: if (level == numlevels - 1) {
467: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
468: PetscCall(DMPlexGlobalVectorView(dm, viewer, NULL, u));
469: }
470: if (level) {
471: PetscInt cStart, cEnd, ccStart, ccEnd, cpStart;
472: DMPolytopeType ct;
473: DMPlexTransform tr;
474: DM sdm;
475: PetscScalar *array;
476: PetscSection section;
477: Vec map;
478: IS gis;
479: const PetscInt *gidx;
481: PetscCall(DMGetCoarseDM(dm, &cdm));
482: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
483: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
484: PetscCall(PetscSectionSetChart(section, cStart, cEnd));
485: for (PetscInt c = cStart; c < cEnd; c++) PetscCall(PetscSectionSetDof(section, c, 1));
486: PetscCall(PetscSectionSetUp(section));
488: PetscCall(DMClone(dm, &sdm));
489: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
490: PetscCall(PetscObjectSetName((PetscObject)section, "pdm_section"));
491: PetscCall(DMSetLocalSection(sdm, section));
492: PetscCall(PetscSectionDestroy(§ion));
494: PetscCall(DMGetLocalVector(sdm, &map));
495: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
496: PetscCall(VecGetArray(map, &array));
497: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
498: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
499: PetscCall(DMPlexTransformSetDM(tr, cdm));
500: PetscCall(DMPlexTransformSetFromOptions(tr));
501: PetscCall(DMPlexTransformSetUp(tr));
502: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
503: PetscCall(DMPlexGetChart(cdm, &cpStart, NULL));
504: PetscCall(DMPlexCreatePointNumbering(cdm, &gis));
505: PetscCall(ISGetIndices(gis, &gidx));
506: for (PetscInt c = ccStart; c < ccEnd; c++) {
507: PetscInt *rsize, *rcone, *rornt, Nt;
508: DMPolytopeType *rct;
509: PetscInt gnum = gidx[c - cpStart] >= 0 ? gidx[c - cpStart] : -(gidx[c - cpStart] + 1);
511: PetscCall(DMPlexGetCellType(cdm, c, &ct));
512: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nt, &rct, &rsize, &rcone, &rornt));
513: for (PetscInt r = 0; r < rsize[Nt - 1]; ++r) {
514: PetscInt pNew;
516: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[Nt - 1], c, r, &pNew));
517: array[pNew - cStart] = gnum;
518: }
519: }
520: PetscCall(ISRestoreIndices(gis, &gidx));
521: PetscCall(ISDestroy(&gis));
522: PetscCall(VecRestoreArray(map, &array));
523: PetscCall(DMPlexTransformDestroy(&tr));
524: PetscCall(DMPlexSectionView(dm, viewer, sdm));
525: PetscCall(DMPlexLocalVectorView(dm, viewer, sdm, map));
526: PetscCall(DMRestoreLocalVector(sdm, &map));
527: PetscCall(DMDestroy(&sdm));
528: }
529: PetscCall(PetscViewerHDF5PopGroup(viewer));
530: PetscCall(DMGetCoarseDM(dm, &dm));
531: }
532: PetscCall(PetscViewerPopFormat(viewer));
533: PetscCall(PetscViewerDestroy(&viewer));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: #else
536: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
537: #endif
538: }
540: static PetscErrorCode LoadFromFile(MPI_Comm comm, const char *filename, DM *odm)
541: {
542: #if defined(PETSC_HAVE_HDF5)
543: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
544: PetscViewer viewer;
545: DM dm, cdm = NULL;
546: PetscSF sfXC = NULL;
547: PetscInt numlevels = -1;
549: PetscFunctionBeginUser;
550: PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer));
551: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numlevels", PETSC_INT, NULL, &numlevels));
552: PetscCall(PetscViewerPushFormat(viewer, format));
553: for (PetscInt level = 0; level < numlevels; level++) {
554: char groupname[PETSC_MAX_PATH_LEN], *dmname;
555: PetscSF sfXB, sfBC, sfG;
556: PetscPartitioner part;
557: PetscInt rr, cc;
558: PetscBool isRegular, isUniform;
560: PetscCall(DMCreate(comm, &dm));
561: PetscCall(DMSetType(dm, DMPLEX));
562: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
563: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
564: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "meshname", PETSC_STRING, NULL, &dmname));
565: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refinelevel", PETSC_INT, NULL, &rr));
566: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, NULL, &cc));
567: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refRegular", PETSC_BOOL, NULL, &isRegular));
568: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refUniform", PETSC_BOOL, NULL, &isUniform));
569: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
570: PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
571: PetscCall(DMPlexLabelsLoad(dm, viewer, sfXB));
572: PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXB));
573: PetscCall(DMPlexGetPartitioner(dm, &part));
574: if (!level) { /* partition the coarse level only */
575: PetscCall(PetscPartitionerSetFromOptions(part));
576: } else { /* propagate partitioning information from coarser to finer level */
577: DM sdm;
578: Vec map;
579: PetscSF sf;
580: PetscLayout clayout;
581: PetscScalar *array;
582: PetscInt *cranks_leaf, *cranks_root, *npoints, *points, *ranks, *starts, *gidxs;
583: PetscInt nparts, cStart, cEnd, nr, ccStart, ccEnd, cpStart, cpEnd;
584: PetscMPIInt size, rank;
586: PetscCall(DMClone(dm, &sdm));
587: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
588: PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXB, NULL, &sf));
589: PetscCall(DMGetLocalVector(sdm, &map));
590: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
591: PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, sf, map));
593: PetscCallMPI(MPI_Comm_size(comm, &size));
594: PetscCallMPI(MPI_Comm_rank(comm, &rank));
595: nparts = size;
596: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
597: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
598: PetscCall(DMPlexGetChart(cdm, &cpStart, &cpEnd));
599: PetscCall(PetscCalloc1(nparts, &npoints));
600: PetscCall(PetscMalloc4(cEnd - cStart, &points, cEnd - cStart, &ranks, nparts + 1, &starts, cEnd - cStart, &gidxs));
601: PetscCall(PetscSFGetGraph(sfXC, &nr, NULL, NULL, NULL));
602: PetscCall(PetscMalloc2(cpEnd - cpStart, &cranks_leaf, nr, &cranks_root));
603: for (PetscInt c = 0; c < cpEnd - cpStart; c++) cranks_leaf[c] = rank;
604: PetscCall(PetscSFReduceBegin(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
605: PetscCall(PetscSFReduceEnd(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
607: PetscCall(VecGetArray(map, &array));
608: for (PetscInt c = 0; c < cEnd - cStart; c++) gidxs[c] = (PetscInt)PetscRealPart(array[c]);
609: PetscCall(VecRestoreArray(map, &array));
611: PetscCall(PetscLayoutCreate(comm, &clayout));
612: PetscCall(PetscLayoutSetLocalSize(clayout, nr));
613: PetscCall(PetscSFSetGraphLayout(sf, clayout, cEnd - cStart, NULL, PETSC_OWN_POINTER, gidxs));
614: PetscCall(PetscLayoutDestroy(&clayout));
616: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
617: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
618: PetscCall(PetscSFDestroy(&sf));
619: PetscCall(PetscFree2(cranks_leaf, cranks_root));
620: for (PetscInt c = 0; c < cEnd - cStart; c++) npoints[ranks[c]]++;
622: starts[0] = 0;
623: for (PetscInt c = 0; c < nparts; c++) starts[c + 1] = starts[c] + npoints[c];
624: for (PetscInt c = 0; c < cEnd - cStart; c++) points[starts[ranks[c]]++] = c;
625: PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
626: PetscCall(PetscPartitionerShellSetPartition(part, nparts, npoints, points));
627: PetscCall(PetscFree(npoints));
628: PetscCall(PetscFree4(points, ranks, starts, gidxs));
629: PetscCall(DMRestoreLocalVector(sdm, &map));
630: PetscCall(DMDestroy(&sdm));
631: }
632: PetscCall(PetscSFDestroy(&sfXC));
633: PetscCall(DMPlexDistribute(dm, 0, &sfBC, odm));
634: if (*odm) {
635: PetscCall(DMDestroy(&dm));
636: dm = *odm;
637: *odm = NULL;
638: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
639: }
640: if (sfBC) PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
641: else {
642: PetscCall(PetscObjectReference((PetscObject)sfXB));
643: sfXC = sfXB;
644: }
645: PetscCall(PetscSFDestroy(&sfXB));
646: PetscCall(PetscSFDestroy(&sfBC));
647: PetscCall(DMSetCoarsenLevel(dm, cc));
648: PetscCall(DMSetRefineLevel(dm, rr));
649: PetscCall(DMPlexSetRegularRefinement(dm, isRegular));
650: PetscCall(DMPlexSetRefinementUniform(dm, isUniform));
651: PetscCall(DMPlexSectionLoad(dm, viewer, NULL, sfXC, &sfG, NULL));
652: if (level == numlevels - 1) {
653: Vec u;
655: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &u));
656: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
657: PetscCall(DMPlexGlobalVectorLoad(dm, viewer, NULL, sfG, u));
658: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &u));
659: }
660: PetscCall(PetscFree(dmname));
661: PetscCall(PetscSFDestroy(&sfG));
662: PetscCall(DMSetCoarseDM(dm, cdm));
663: PetscCall(DMDestroy(&cdm));
664: PetscCall(PetscViewerHDF5PopGroup(viewer));
665: cdm = dm;
666: }
667: *odm = dm;
668: PetscCall(PetscViewerPopFormat(viewer));
669: PetscCall(PetscViewerDestroy(&viewer));
670: PetscCall(PetscSFDestroy(&sfXC));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: #else
673: SETERRQ(comm, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
674: #endif
675: }
677: /* Project source function and make it zero-mean */
678: static PetscErrorCode ProjectSource(DM dm, PetscReal time, AppCtx *ctx)
679: {
680: PetscInt id = C_FIELD_ID;
681: DM dmAux;
682: Vec u, lu;
683: IS is;
684: void *ctxs[NUM_FIELDS];
685: PetscScalar vals[NUM_FIELDS];
686: PetscDS ds;
687: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
689: PetscFunctionBeginUser;
690: switch (ctx->source_num) {
691: case 0:
692: funcs[P_FIELD_ID] = source_0;
693: ctxs[P_FIELD_ID] = ctx->x0;
694: break;
695: case 1:
696: funcs[P_FIELD_ID] = source_1;
697: ctxs[P_FIELD_ID] = NULL;
698: break;
699: default:
700: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown source");
701: }
702: funcs[C_FIELD_ID] = zerof;
703: ctxs[C_FIELD_ID] = NULL;
704: PetscCall(DMGetDS(dm, &ds));
705: PetscCall(DMGetGlobalVector(dm, &u));
706: PetscCall(DMProjectFunction(dm, time, funcs, ctxs, INSERT_ALL_VALUES, u));
707: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
708: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
709: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
710: PetscCall(VecShift(u, -vals[P_FIELD_ID]));
711: PetscCall(DMCreateSubDM(dm, 1, &id, &is, NULL));
712: PetscCall(VecISSet(u, is, 0));
713: PetscCall(ISDestroy(&is));
715: /* Attach source vector as auxiliary vector:
716: Use a different DM to break ref cycles */
717: PetscCall(DMClone(dm, &dmAux));
718: PetscCall(DMCopyDisc(dm, dmAux));
719: PetscCall(DMCreateLocalVector(dmAux, &lu));
720: PetscCall(DMDestroy(&dmAux));
721: PetscCall(DMGlobalToLocal(dm, u, INSERT_VALUES, lu));
722: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, lu));
723: PetscCall(VecViewFromOptions(lu, NULL, "-aux_view"));
724: PetscCall(VecDestroy(&lu));
725: PetscCall(DMRestoreGlobalVector(dm, &u));
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /* callback for the creation of the potential null space */
730: static PetscErrorCode CreatePotentialNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
731: {
732: Vec vec;
733: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zerof};
735: PetscFunctionBeginUser;
736: funcs[nfield] = constantf;
737: PetscCall(DMCreateGlobalVector(dm, &vec));
738: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
739: PetscCall(VecNormalize(vec, NULL));
740: PetscCall(PetscObjectSetName((PetscObject)vec, "Potential Null Space"));
741: PetscCall(VecViewFromOptions(vec, NULL, "-potential_nullspace_view"));
742: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
743: /* break ref cycles */
744: PetscCall(VecSetDM(vec, NULL));
745: PetscCall(VecDestroy(&vec));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /* customize residuals and Jacobians */
750: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
751: {
752: PetscDS ds;
753: PetscInt cdim, dim;
754: PetscScalar constants[NUM_CONSTANTS];
756: PetscFunctionBeginUser;
757: constants[R_ID] = ctx->r;
758: constants[EPS_ID] = ctx->eps;
759: constants[ALPHA_ID] = ctx->alpha;
760: constants[GAMMA_ID] = ctx->gamma;
761: constants[D_ID] = ctx->D;
762: constants[C2_ID] = ctx->c * ctx->c;
764: PetscCall(DMGetDimension(dm, &dim));
765: PetscCall(DMGetCoordinateDim(dm, &cdim));
766: PetscCheck(dim == 2 && cdim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2D meshes");
767: PetscCall(DMGetDS(dm, &ds));
768: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
769: PetscCall(PetscDSSetImplicit(ds, C_FIELD_ID, PETSC_TRUE));
770: PetscCall(PetscDSSetImplicit(ds, P_FIELD_ID, PETSC_TRUE));
771: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
772: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
773: PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0, C_1));
774: PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1));
775: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0, NULL, NULL, JC_1_c1c1));
776: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1, NULL, NULL));
777: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0, NULL));
778: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1));
780: /* Attach potential nullspace */
781: PetscCall(DMSetNullSpaceConstructor(dm, P_FIELD_ID, CreatePotentialNullSpace));
783: /* Attach source function as auxiliary vector */
784: PetscCall(ProjectSource(dm, 0, ctx));
786: /* Add callbacks */
787: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL));
788: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, NULL));
789: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, NULL));
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /* setup discrete spaces and residuals */
794: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
795: {
796: DM plex, cdm = dm;
797: PetscFE feC, feP;
798: PetscBool simplex;
799: PetscInt dim;
800: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
801: MatNullSpace nsp;
803: PetscFunctionBeginUser;
804: PetscCall(DMGetDimension(dm, &dim));
806: PetscCall(DMConvert(dm, DMPLEX, &plex));
807: PetscCall(DMPlexIsSimplex(plex, &simplex));
808: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &simplex, 1, MPIU_BOOL, MPI_LOR, comm));
809: PetscCall(DMDestroy(&plex));
811: /* We model Cij in H^1 with Cij = Cji -> dim*(dim+1)/2 components */
812: PetscCall(PetscFECreateDefault(comm, dim, (dim * (dim + 1)) / 2, simplex, "c_", -1, &feC));
813: PetscCall(PetscObjectSetName((PetscObject)feC, "conductivity"));
814: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "p_", -1, &feP));
815: PetscCall(PetscObjectSetName((PetscObject)feP, "potential"));
816: PetscCall(MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, &nsp));
817: PetscCall(PetscObjectCompose((PetscObject)feP, "nullspace", (PetscObject)nsp));
818: PetscCall(MatNullSpaceDestroy(&nsp));
819: PetscCall(PetscFECopyQuadrature(feC, feP));
821: PetscCall(DMSetNumFields(dm, 2));
822: PetscCall(DMSetField(dm, C_FIELD_ID, NULL, (PetscObject)feC));
823: PetscCall(DMSetField(dm, P_FIELD_ID, NULL, (PetscObject)feP));
824: PetscCall(PetscFEDestroy(&feC));
825: PetscCall(PetscFEDestroy(&feP));
826: PetscCall(DMCreateDS(dm));
827: while (cdm) {
828: PetscCall(DMCopyDisc(dm, cdm));
829: PetscCall(SetupProblem(cdm, ctx));
830: PetscCall(DMGetCoarseDM(cdm, &cdm));
831: }
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /* Create mesh by command line options */
836: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
837: {
838: PetscFunctionBeginUser;
839: if (ctx->load) {
840: PetscInt refine = 0;
841: PetscBool isHierarchy;
842: DM *dms;
843: char typeName[256];
844: PetscBool flg;
846: PetscCall(LoadFromFile(comm, ctx->load_filename, dm));
847: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
848: PetscCall(PetscOptionsFList("-dm_mat_type", "Matrix type used for created matrices", "DMSetMatType", MatList, MATAIJ, typeName, sizeof(typeName), &flg));
849: if (flg) PetscCall(DMSetMatType(*dm, typeName));
850: PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0));
851: PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0));
852: PetscOptionsEnd();
853: if (refine) {
854: PetscCall(SetupDiscretization(*dm, ctx));
855: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
856: }
857: PetscCall(PetscCalloc1(refine, &dms));
858: if (isHierarchy) PetscCall(DMRefineHierarchy(*dm, refine, dms));
859: for (PetscInt r = 0; r < refine; r++) {
860: Mat M;
861: DM dmr = dms[r];
862: Vec u, ur;
864: if (!isHierarchy) {
865: PetscCall(DMRefine(*dm, PetscObjectComm((PetscObject)dm), &dmr));
866: PetscCall(DMSetCoarseDM(dmr, *dm));
867: }
868: PetscCall(DMCreateInterpolation(*dm, dmr, &M, NULL));
869: PetscCall(DMGetNamedGlobalVector(*dm, "solution_", &u));
870: PetscCall(DMGetNamedGlobalVector(dmr, "solution_", &ur));
871: PetscCall(MatInterpolate(M, u, ur));
872: PetscCall(DMRestoreNamedGlobalVector(*dm, "solution_", &u));
873: PetscCall(DMRestoreNamedGlobalVector(dmr, "solution_", &ur));
874: PetscCall(MatDestroy(&M));
875: if (!isHierarchy) PetscCall(DMSetCoarseDM(dmr, NULL));
876: PetscCall(DMDestroy(dm));
877: *dm = dmr;
878: }
879: if (refine && !isHierarchy) PetscCall(DMSetRefineLevel(*dm, 0));
880: PetscCall(PetscFree(dms));
881: } else {
882: PetscCall(DMCreate(comm, dm));
883: PetscCall(DMSetType(*dm, DMPLEX));
884: PetscCall(DMSetFromOptions(*dm));
885: PetscCall(DMLocalizeCoordinates(*dm));
886: {
887: char convType[256];
888: PetscBool flg;
889: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
890: PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg));
891: PetscOptionsEnd();
892: if (flg) {
893: DM dmConv;
894: PetscCall(DMConvert(*dm, convType, &dmConv));
895: if (dmConv) {
896: PetscCall(DMDestroy(dm));
897: *dm = dmConv;
898: PetscCall(DMSetFromOptions(*dm));
899: PetscCall(DMSetUp(*dm));
900: }
901: }
902: }
903: }
904: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: /* Make potential field zero mean */
909: static PetscErrorCode ZeroMeanPotential(DM dm, Vec u)
910: {
911: PetscScalar vals[NUM_FIELDS];
912: PetscDS ds;
913: IS is;
915: PetscFunctionBeginUser;
916: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
917: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
918: PetscCall(DMGetDS(dm, &ds));
919: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
920: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
921: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
922: PetscCall(VecISShift(u, is, -vals[P_FIELD_ID]));
923: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: /* Compute initial conditions and exclude potential from local truncation error
928: Since we are solving a DAE, once the initial conditions for the differential
929: variables are set, we need to compute the corresponding value for the
930: algebraic variables. We do so by creating a subDM for the potential only
931: and solve a static problem with SNES */
932: static PetscErrorCode SetInitialConditionsAndTolerances(TS ts, PetscInt nv, Vec vecs[], PetscBool valid)
933: {
934: DM dm;
935: Vec tu, u, p, lsource, subaux, vatol, vrtol;
936: PetscReal t, atol, rtol;
937: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
938: IS isp;
939: DM dmp;
940: VecScatter sctp = NULL;
941: PetscDS ds;
942: SNES snes;
943: KSP ksp;
944: PC pc;
945: AppCtx *ctx;
947: PetscFunctionBeginUser;
948: PetscCall(TSGetDM(ts, &dm));
949: PetscCall(TSGetApplicationContext(ts, &ctx));
950: if (nv > 1) {
951: PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, NULL));
952: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
953: PetscCall(DMCreateGlobalVector(dm, &vatol));
954: PetscCall(DMCreateGlobalVector(dm, &vrtol));
955: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
956: PetscCall(VecSet(vatol, atol));
957: PetscCall(VecISSet(vatol, isp, -1));
958: PetscCall(VecSet(vrtol, rtol));
959: PetscCall(VecISSet(vrtol, isp, -1));
960: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
961: PetscCall(VecDestroy(&vatol));
962: PetscCall(VecDestroy(&vrtol));
963: PetscCall(ISDestroy(&isp));
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
966: PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, &dmp));
967: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
968: PetscCall(DMSetMatType(dmp, MATAIJ));
969: PetscCall(DMGetDS(dmp, &ds));
970: //PetscCall(PetscDSSetResidual(ds, 0, P_0, P_1_aux));
971: PetscCall(PetscDSSetResidual(ds, 0, 0, P_1_aux));
972: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux));
973: PetscCall(DMPlexSetSNESLocalFEM(dmp, PETSC_FALSE, NULL));
975: PetscCall(DMCreateGlobalVector(dmp, &p));
977: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dmp), &snes));
978: PetscCall(SNESSetDM(snes, dmp));
979: PetscCall(SNESSetOptionsPrefix(snes, "initial_"));
980: PetscCall(SNESSetErrorIfNotConverged(snes, PETSC_TRUE));
981: PetscCall(SNESGetKSP(snes, &ksp));
982: PetscCall(KSPGetPC(ksp, &pc));
983: PetscCall(PCSetType(pc, PCGAMG));
984: PetscCall(SNESSetFromOptions(snes));
985: PetscCall(SNESSetUp(snes));
987: /* Loop over input vectors and compute corresponding potential */
988: for (PetscInt i = 0; i < nv; i++) {
989: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
991: u = vecs[i];
992: if (!valid) { /* Assumes entries in u are not valid */
993: PetscCall(TSGetTime(ts, &t));
994: switch (ctx->ic_num) {
995: case 0:
996: funcs[C_FIELD_ID] = initial_conditions_C_0;
997: break;
998: case 1:
999: funcs[C_FIELD_ID] = initial_conditions_C_1;
1000: break;
1001: case 2:
1002: funcs[C_FIELD_ID] = initial_conditions_C_2;
1003: break;
1004: default:
1005: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unknown IC");
1006: }
1007: funcs[P_FIELD_ID] = zerof;
1008: PetscCall(DMProjectFunction(dm, t, funcs, NULL, INSERT_ALL_VALUES, u));
1009: }
1011: /* pass conductivity and source information via auxiliary data */
1012: PetscCall(DMGetGlobalVector(dm, &tu));
1013: PetscCall(VecCopy(u, tu));
1014: PetscCall(VecISSet(tu, isp, 0.0));
1015: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &lsource));
1016: PetscCall(DMCreateLocalVector(dm, &subaux));
1017: PetscCall(DMGlobalToLocal(dm, tu, INSERT_VALUES, subaux));
1018: PetscCall(DMRestoreGlobalVector(dm, &tu));
1019: PetscCall(VecAXPY(subaux, 1.0, lsource));
1020: PetscCall(VecViewFromOptions(subaux, NULL, "-initial_aux_view"));
1021: PetscCall(DMSetAuxiliaryVec(dmp, NULL, 0, 0, subaux));
1022: PetscCall(VecDestroy(&subaux));
1024: /* solve the subproblem */
1025: if (!sctp) PetscCall(VecScatterCreate(u, isp, p, NULL, &sctp));
1026: PetscCall(VecScatterBegin(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1027: PetscCall(VecScatterEnd(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1028: PetscCall(SNESSolve(snes, NULL, p));
1030: /* scatter from potential only to full space */
1031: PetscCall(VecScatterBegin(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1032: PetscCall(VecScatterEnd(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1033: PetscCall(ZeroMeanPotential(dm, u));
1034: }
1035: PetscCall(VecDestroy(&p));
1036: PetscCall(DMDestroy(&dmp));
1037: PetscCall(SNESDestroy(&snes));
1038: PetscCall(VecScatterDestroy(&sctp));
1040: /* exclude potential from computation of the LTE */
1041: PetscCall(DMCreateGlobalVector(dm, &vatol));
1042: PetscCall(DMCreateGlobalVector(dm, &vrtol));
1043: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
1044: PetscCall(VecSet(vatol, atol));
1045: PetscCall(VecISSet(vatol, isp, -1));
1046: PetscCall(VecSet(vrtol, rtol));
1047: PetscCall(VecISSet(vrtol, isp, -1));
1048: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
1049: PetscCall(VecDestroy(&vatol));
1050: PetscCall(VecDestroy(&vrtol));
1051: PetscCall(ISDestroy(&isp));
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: /* mesh adaption context */
1056: typedef struct {
1057: VecTagger refineTag;
1058: DMLabel adaptLabel;
1059: PetscInt cnt;
1060: } AdaptCtx;
1062: static PetscErrorCode ResizeSetUp(TS ts, PetscInt nstep, PetscReal time, Vec u, PetscBool *resize, void *vctx)
1063: {
1064: AdaptCtx *ctx = (AdaptCtx *)vctx;
1065: Vec ellVecCells, ellVecCellsF;
1066: DM dm, plex;
1067: PetscDS ds;
1068: PetscReal norm;
1069: PetscInt cStart, cEnd;
1071: PetscFunctionBeginUser;
1072: PetscCall(TSGetDM(ts, &dm));
1073: PetscCall(DMConvert(dm, DMPLEX, &plex));
1074: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
1075: PetscCall(DMDestroy(&plex));
1076: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), NUM_FIELDS * (cEnd - cStart), PETSC_DECIDE, &ellVecCellsF));
1077: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), cEnd - cStart, PETSC_DECIDE, &ellVecCells));
1078: PetscCall(VecSetBlockSize(ellVecCellsF, NUM_FIELDS));
1079: PetscCall(DMGetDS(dm, &ds));
1080: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
1081: PetscCall(DMPlexComputeCellwiseIntegralFEM(dm, u, ellVecCellsF, NULL));
1082: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1083: PetscCall(VecStrideGather(ellVecCellsF, C_FIELD_ID, ellVecCells, INSERT_VALUES));
1084: PetscCall(VecDestroy(&ellVecCellsF));
1085: PetscCall(VecNorm(ellVecCells, NORM_1, &norm));
1086: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "STEP %d norm %g\n", (int)nstep, (double)norm));
1087: if (norm && !ctx->cnt) {
1088: IS refineIS;
1090: *resize = PETSC_TRUE;
1091: if (!ctx->refineTag) {
1092: VecTaggerBox refineBox;
1093: refineBox.min = PETSC_MACHINE_EPSILON;
1094: refineBox.max = PETSC_MAX_REAL;
1096: PetscCall(VecTaggerCreate(PETSC_COMM_SELF, &ctx->refineTag));
1097: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->refineTag, "refine_"));
1098: PetscCall(VecTaggerSetType(ctx->refineTag, VECTAGGERABSOLUTE));
1099: PetscCall(VecTaggerAbsoluteSetBox(ctx->refineTag, &refineBox));
1100: PetscCall(VecTaggerSetFromOptions(ctx->refineTag));
1101: PetscCall(VecTaggerSetUp(ctx->refineTag));
1102: PetscCall(PetscObjectViewFromOptions((PetscObject)ctx->refineTag, NULL, "-tag_view"));
1103: }
1104: PetscCall(DMLabelDestroy(&ctx->adaptLabel));
1105: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)ts), "adapt", &ctx->adaptLabel));
1106: PetscCall(VecTaggerComputeIS(ctx->refineTag, ellVecCells, &refineIS, NULL));
1107: PetscCall(DMLabelSetStratumIS(ctx->adaptLabel, DM_ADAPT_REFINE, refineIS));
1108: PetscCall(ISDestroy(&refineIS));
1109: #if 0
1110: 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[]);
1111: Vec ellVec;
1113: funcs[P_FIELD_ID] = ellipticity_fail;
1114: funcs[C_FIELD_ID] = NULL;
1116: PetscCall(DMGetGlobalVector(dm, &ellVec));
1117: PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
1118: PetscCall(VecViewFromOptions(ellVec,NULL,"-view_amr_ell"));
1119: PetscCall(DMRestoreGlobalVector(dm, &ellVec));
1120: #endif
1121: ctx->cnt++;
1122: } else {
1123: ctx->cnt = 0;
1124: }
1125: PetscCall(VecDestroy(&ellVecCells));
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: static PetscErrorCode ResizeTransfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *vctx)
1130: {
1131: AdaptCtx *actx = (AdaptCtx *)vctx;
1132: AppCtx *ctx;
1133: DM dm, adm;
1134: PetscReal time;
1136: PetscFunctionBeginUser;
1137: PetscCall(TSGetDM(ts, &dm));
1138: PetscCheck(actx->adaptLabel, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptLabel");
1139: PetscCall(DMAdaptLabel(dm, actx->adaptLabel, &adm));
1140: PetscCheck(adm, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adapted DM");
1141: PetscCall(TSGetTime(ts, &time));
1142: PetscCall(DMLabelDestroy(&actx->adaptLabel));
1143: for (PetscInt i = 0; i < nv; i++) {
1144: PetscCall(DMCreateGlobalVector(adm, &vecsout[i]));
1145: PetscCall(DMForestTransferVec(dm, vecsin[i], adm, vecsout[i], PETSC_TRUE, time));
1146: }
1147: PetscCall(DMForestSetAdaptivityForest(adm, NULL));
1148: PetscCall(DMSetCoarseDM(adm, NULL));
1149: PetscCall(DMSetLocalSection(adm, NULL));
1150: PetscCall(TSSetDM(ts, adm));
1151: PetscCall(TSGetTime(ts, &time));
1152: PetscCall(TSGetApplicationContext(ts, &ctx));
1153: PetscCall(ProjectSource(adm, time, ctx));
1154: PetscCall(SetInitialConditionsAndTolerances(ts, nv, vecsout, PETSC_TRUE));
1155: PetscCall(DMDestroy(&adm));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1160: {
1161: PetscFunctionBeginUser;
1162: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1163: PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
1164: PetscCall(PetscViewerFileSetName(*viewer, filename));
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: /* Monitor relevant functionals */
1169: static PetscErrorCode Monitor(TS ts, PetscInt stepnum, PetscReal time, Vec u, void *vctx)
1170: {
1171: PetscScalar vals[2 * NUM_FIELDS];
1172: DM dm;
1173: PetscDS ds;
1174: SNES snes;
1175: PetscInt nits, lits;
1176: AppCtx *ctx;
1178: PetscFunctionBeginUser;
1179: PetscCall(TSGetDM(ts, &dm));
1180: PetscCall(TSGetApplicationContext(ts, &ctx));
1181: PetscCall(DMGetDS(dm, &ds));
1183: /* monitor energy and potential average */
1184: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
1185: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1186: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1188: /* monitor ellipticity_fail */
1189: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
1190: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals + NUM_FIELDS, NULL));
1191: if (ctx->ellipticity) {
1192: 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[]);
1193: Vec ellVec;
1194: PetscViewer viewer;
1195: char filename[PETSC_MAX_PATH_LEN];
1197: funcs[P_FIELD_ID] = ellipticity_fail;
1198: funcs[C_FIELD_ID] = NULL;
1200: PetscCall(DMGetGlobalVector(dm, &ellVec));
1201: PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
1202: PetscCall(PetscSNPrintf(filename, sizeof filename, "ellipticity_fail-%03" PetscInt_FMT ".vtu", stepnum));
1203: PetscCall(OutputVTK(dm, filename, &viewer));
1204: PetscCall(VecView(ellVec, viewer));
1205: PetscCall(PetscViewerDestroy(&viewer));
1206: PetscCall(DMRestoreGlobalVector(dm, &ellVec));
1207: }
1208: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1210: /* monitor linear and nonlinear iterations */
1211: PetscCall(TSGetSNES(ts, &snes));
1212: PetscCall(SNESGetIterationNumber(snes, &nits));
1213: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1215: 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])));
1216: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), " nonlinear its %" PetscInt_FMT ", linear its %" PetscInt_FMT "\n", nits, lits));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: /* Save restart information */
1221: static PetscErrorCode MonitorSave(TS ts, PetscInt steps, PetscReal time, Vec u, void *vctx)
1222: {
1223: DM dm;
1224: AppCtx *ctx = (AppCtx *)vctx;
1225: PetscInt save_every = ctx->save_every;
1226: TSConvergedReason reason;
1228: PetscFunctionBeginUser;
1229: if (!ctx->save) PetscFunctionReturn(PETSC_SUCCESS);
1230: PetscCall(TSGetDM(ts, &dm));
1231: PetscCall(TSGetConvergedReason(ts, &reason));
1232: if ((save_every > 0 && steps % save_every == 0) || (save_every == -1 && reason) || save_every < -1) PetscCall(SaveToFile(dm, u, ctx->save_filename));
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /* Make potential zero mean after SNES solve */
1237: static PetscErrorCode PostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
1238: {
1239: DM dm;
1240: Vec u = Y[stageindex];
1242: PetscFunctionBeginUser;
1243: PetscCall(TSGetDM(ts, &dm));
1244: PetscCall(ZeroMeanPotential(dm, u));
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: static PetscErrorCode VecViewFlux(Vec u, const char *opts)
1249: {
1250: Vec fluxVec;
1251: DM dmFlux, dm, plex;
1252: PetscInt dim;
1253: PetscFE feC, feFluxC, feNormC;
1254: PetscBool simplex, has;
1256: 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};
1258: PetscFunctionBeginUser;
1259: PetscCall(PetscOptionsHasName(NULL, NULL, opts, &has));
1260: if (!has) PetscFunctionReturn(PETSC_SUCCESS);
1261: PetscCall(VecGetDM(u, &dm));
1262: PetscCall(DMGetDimension(dm, &dim));
1263: PetscCall(DMGetField(dm, C_FIELD_ID, NULL, (PetscObject *)&feC));
1264: PetscCall(DMConvert(dm, DMPLEX, &plex));
1265: PetscCall(DMPlexIsSimplex(plex, &simplex));
1266: PetscCall(DMDestroy(&plex));
1267: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, simplex, "flux_", -1, &feFluxC));
1268: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, "normc_", -1, &feNormC));
1269: PetscCall(PetscFECopyQuadrature(feC, feFluxC));
1270: PetscCall(PetscFECopyQuadrature(feC, feNormC));
1271: PetscCall(DMClone(dm, &dmFlux));
1272: PetscCall(DMSetNumFields(dmFlux, 1));
1273: PetscCall(DMSetField(dmFlux, 0, NULL, (PetscObject)feNormC));
1274: /* paraview segfaults! */
1275: //PetscCall(DMSetField(dmFlux, 1, NULL, (PetscObject)feFluxC));
1276: PetscCall(DMCreateDS(dmFlux));
1277: PetscCall(PetscFEDestroy(&feFluxC));
1278: PetscCall(PetscFEDestroy(&feNormC));
1280: PetscCall(DMGetGlobalVector(dmFlux, &fluxVec));
1281: PetscCall(DMProjectField(dmFlux, 0.0, u, funcs, INSERT_VALUES, fluxVec));
1282: PetscCall(VecViewFromOptions(fluxVec, NULL, opts));
1283: PetscCall(DMRestoreGlobalVector(dmFlux, &fluxVec));
1284: PetscCall(DMDestroy(&dmFlux));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: static PetscErrorCode Run(MPI_Comm comm, AppCtx *ctx)
1289: {
1290: DM dm;
1291: TS ts;
1292: Vec u;
1293: AdaptCtx *actx;
1294: PetscBool flg;
1296: PetscFunctionBeginUser;
1297: if (ctx->test_restart) {
1298: PetscViewer subviewer;
1299: PetscMPIInt rank;
1301: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1302: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1303: if (ctx->load) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d loading from %s\n", rank, ctx->load_filename));
1304: if (ctx->save) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d saving to %s\n", rank, ctx->save_filename));
1305: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1306: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1307: } else {
1308: PetscCall(PetscPrintf(comm, "----------------------------\n"));
1309: PetscCall(PetscPrintf(comm, "Simulation parameters:\n"));
1310: PetscCall(PetscPrintf(comm, " r : %g\n", (double)ctx->r));
1311: PetscCall(PetscPrintf(comm, " eps : %g\n", (double)ctx->eps));
1312: PetscCall(PetscPrintf(comm, " alpha: %g\n", (double)ctx->alpha));
1313: PetscCall(PetscPrintf(comm, " gamma: %g\n", (double)ctx->gamma));
1314: PetscCall(PetscPrintf(comm, " D : %g\n", (double)ctx->D));
1315: PetscCall(PetscPrintf(comm, " c : %g\n", (double)ctx->c));
1316: if (ctx->load) PetscCall(PetscPrintf(comm, " load : %s\n", ctx->load_filename));
1317: else PetscCall(PetscPrintf(comm, " IC : %" PetscInt_FMT "\n", ctx->ic_num));
1318: PetscCall(PetscPrintf(comm, " S : %" PetscInt_FMT "\n", ctx->source_num));
1319: PetscCall(PetscPrintf(comm, " x0 : (%g,%g)\n", (double)ctx->x0[0], (double)ctx->x0[1]));
1320: PetscCall(PetscPrintf(comm, "----------------------------\n"));
1321: }
1323: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SetupStage));
1324: PetscCall(CreateMesh(comm, &dm, ctx));
1325: PetscCall(SetupDiscretization(dm, ctx));
1327: PetscCall(TSCreate(comm, &ts));
1328: PetscCall(TSSetApplicationContext(ts, ctx));
1330: PetscCall(TSSetDM(ts, dm));
1331: if (ctx->test_restart) {
1332: PetscViewer subviewer;
1333: PetscMPIInt rank;
1334: PetscInt level;
1336: PetscCall(DMGetRefineLevel(dm, &level));
1337: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1338: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1339: PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d DM refinement level %" PetscInt_FMT "\n", rank, level));
1340: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
1341: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1342: }
1344: if (ctx->test_restart) PetscCall(TSSetMaxSteps(ts, 1));
1345: PetscCall(TSSetMaxTime(ts, 10.0));
1346: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1347: if (!ctx->test_restart) PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
1348: PetscCall(TSMonitorSet(ts, MonitorSave, ctx, NULL));
1349: PetscCall(PetscNew(&actx));
1350: if (ctx->amr) PetscCall(TSSetResize(ts, PETSC_TRUE, ResizeSetUp, ResizeTransfer, actx));
1351: PetscCall(TSSetPostStage(ts, PostStage));
1352: PetscCall(TSSetMaxSNESFailures(ts, -1));
1353: PetscCall(TSSetFromOptions(ts));
1355: PetscCall(DMCreateGlobalVector(dm, &u));
1356: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
1357: PetscCall(DMHasNamedGlobalVector(dm, "solution_", &flg));
1358: if (flg) {
1359: Vec ru;
1361: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &ru));
1362: PetscCall(VecCopy(ru, u));
1363: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &ru));
1364: }
1365: PetscCall(SetInitialConditionsAndTolerances(ts, 1, &u, PETSC_FALSE));
1366: PetscCall(TSSetSolution(ts, u));
1367: PetscCall(VecDestroy(&u));
1368: PetscCall(DMDestroy(&dm));
1369: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
1371: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SolveStage));
1372: PetscCall(TSSolve(ts, NULL));
1373: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
1375: PetscCall(TSGetSolution(ts, &u));
1376: PetscCall(VecViewFromOptions(u, NULL, "-final_view"));
1377: PetscCall(VecViewFlux(u, "-final_flux_view"));
1379: PetscCall(TSDestroy(&ts));
1380: PetscCall(VecTaggerDestroy(&actx->refineTag));
1381: PetscCall(PetscFree(actx));
1382: PetscFunctionReturn(PETSC_SUCCESS);
1383: }
1385: int main(int argc, char **argv)
1386: {
1387: AppCtx ctx;
1389: PetscFunctionBeginUser;
1390: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1391: PetscCall(ProcessOptions(&ctx));
1392: PetscCall(PetscLogStageRegister("Setup", &SetupStage));
1393: PetscCall(PetscLogStageRegister("Solve", &SolveStage));
1394: if (ctx.test_restart) { /* Test sequences of save and loads */
1395: PetscMPIInt rank;
1397: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1399: /* test saving */
1400: ctx.load = PETSC_FALSE;
1401: ctx.save = PETSC_TRUE;
1402: /* sequential save */
1403: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential save\n"));
1404: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_seq_%d.h5", rank));
1405: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1406: /* parallel save */
1407: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel save\n"));
1408: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_par.h5"));
1409: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1411: /* test loading */
1412: ctx.load = PETSC_TRUE;
1413: ctx.save = PETSC_FALSE;
1414: /* sequential and parallel runs from sequential save */
1415: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_seq_0.h5"));
1416: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from sequential save\n"));
1417: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1418: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from sequential save\n"));
1419: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1420: /* sequential and parallel runs from parallel save */
1421: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_par.h5"));
1422: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from parallel save\n"));
1423: PetscCall(Run(PETSC_COMM_SELF, &ctx));
1424: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from parallel save\n"));
1425: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1426: } else { /* Run the simulation */
1427: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
1428: }
1429: PetscCall(PetscFinalize());
1430: return 0;
1431: }
1433: /*TEST
1435: testset:
1436: 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
1438: test:
1439: suffix: 0
1440: nsize: {{1 2}}
1441: args: -dm_refine 1
1443: test:
1444: suffix: 0_dirk
1445: nsize: {{1 2}}
1446: args: -dm_refine 1 -ts_type dirk
1448: test:
1449: suffix: 0_dirk_mg
1450: nsize: {{1 2}}
1451: 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
1453: test:
1454: requires: p4est
1455: suffix: 0_p4est
1456: nsize: {{1 2}}
1457: args: -dm_refine 1 -dm_plex_convert_type p4est
1459: test:
1460: suffix: 0_periodic
1461: nsize: {{1 2}}
1462: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1
1464: test:
1465: requires: p4est
1466: suffix: 0_p4est_periodic
1467: nsize: {{1 2}}
1468: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -dm_plex_convert_type p4est
1470: test:
1471: requires: p4est
1472: suffix: 0_p4est_mg
1473: nsize: {{1 2}}
1474: 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
1476: testset:
1477: requires: hdf5
1478: 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
1480: test:
1481: requires: !single
1482: suffix: restart
1483: nsize: {{1 2}separate output}
1484: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 0
1486: test:
1487: requires: triangle
1488: suffix: restart_simplex
1489: nsize: {{1 2}separate output}
1490: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 1
1492: test:
1493: requires: !single
1494: suffix: restart_refonly
1495: nsize: {{1 2}separate output}
1496: args: -dm_refine 1 -dm_plex_simplex 0
1498: test:
1499: requires: triangle
1500: suffix: restart_simplex_refonly
1501: nsize: {{1 2}separate output}
1502: args: -dm_refine 1 -dm_plex_simplex 1
1504: TEST*/