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^d:
14: * dC/dt - D^2 \Delta C - \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 dxd 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: \alpha = metabolic coefficient
29: \gamma = metabolic exponent
30: r, eps are regularization parameters
32: We use Lagrange elements for C_ij and P.
33: Equations are rescaled to obtain a symmetric Jacobian.
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: FIXC_ID,
49: SPLIT_ID,
50: NUM_CONSTANTS
51: } ConstantIdx;
53: PetscLogStage SetupStage, SolveStage;
55: #define NORM2C(c00, c01, c11) (PetscSqr(c00) + 2 * PetscSqr(c01) + PetscSqr(c11))
56: #define NORM2C_3d(c00, c01, c02, c11, c12, c22) (PetscSqr(c00) + 2 * PetscSqr(c01) + 2 * PetscSqr(c02) + PetscSqr(c11) + 2 * PetscSqr(c12) + PetscSqr(c22))
58: /* Eigenvalues real 3x3 symmetric matrix https://onlinelibrary.wiley.com/doi/full/10.1002/nme.7153 */
59: #if !PetscDefined(USE_COMPLEX)
60: static inline void Eigenvalues_Sym3x3(PetscReal a11, PetscReal a12, PetscReal a13, PetscReal a22, PetscReal a23, PetscReal a33, PetscReal x[3])
61: {
62: const PetscBool td = (PetscBool)(a13 == 0 && a23 == 0);
63: if (td) { /* two-dimensional case */
64: PetscReal a12_2 = PetscSqr(a12);
65: PetscReal delta = PetscSqr(a11 - a22) + 4 * a12_2;
66: PetscReal b = -(a11 + a22);
67: PetscReal c = a11 * a22 - a12_2;
68: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(delta));
70: x[0] = temp;
71: x[1] = c / temp;
72: x[2] = a33;
73: } else {
74: const PetscReal I1 = a11 + a22 + a33;
75: const PetscReal J2 = (PetscSqr(a11 - a22) + PetscSqr(a22 - a33) + PetscSqr(a33 - a11)) / 6 + PetscSqr(a12) + PetscSqr(a23) + PetscSqr(a13);
76: const PetscReal s = PetscSqrtReal(J2 / 3);
77: const PetscReal tol = PETSC_MACHINE_EPSILON;
79: if (s < tol) {
80: x[0] = x[1] = x[2] = 0.0;
81: } else {
82: const PetscReal S[] = {a11 - I1 / 3, a12, a13, a22 - I1 / 3, a23, a33 - I1 / 3};
84: /* T = S^2 */
85: PetscReal T[6];
86: T[0] = S[0] * S[0] + S[1] * S[1] + S[2] * S[2];
87: T[1] = S[0] * S[1] + S[1] * S[3] + S[2] * S[4];
88: T[2] = S[0] * S[2] + S[1] * S[4] + S[2] * S[5];
89: T[3] = S[1] * S[1] + S[3] * S[3] + S[4] * S[4];
90: T[4] = S[1] * S[2] + S[3] * S[4] + S[4] * S[5];
91: T[5] = S[2] * S[2] + S[4] * S[4] + S[5] * S[5];
93: T[0] = T[0] - 2.0 * J2 / 3.0;
94: T[3] = T[3] - 2.0 * J2 / 3.0;
95: T[5] = T[5] - 2.0 * J2 / 3.0;
97: const PetscReal aa = NORM2C_3d(T[0] - s * S[0], T[1] - s * S[1], T[2] - s * S[2], T[3] - s * S[3], T[4] - s * S[4], T[5] - s * S[5]);
98: const PetscReal bb = NORM2C_3d(T[0] + s * S[0], T[1] + s * S[1], T[2] + s * S[2], T[3] + s * S[3], T[4] + s * S[4], T[5] + s * S[5]);
99: const PetscReal d = PetscSqrtReal(aa / bb);
100: const PetscBool sj = (PetscBool)(1.0 - d > 0.0);
102: if (PetscAbsReal(1 - d) < tol) {
103: x[0] = -PetscSqrtReal(J2);
104: x[1] = 0.0;
105: x[2] = PetscSqrtReal(J2);
106: } else {
107: const PetscReal sjn = sj ? 1.0 : -1.0;
108: //const PetscReal atanarg = sj ? d : 1.0 / d;
109: //const PetscReal alpha = 2.0 * PetscAtanReal(atanarg) / 3.0;
110: const PetscReal atanval = d > tol ? 2.0 * PetscAtanReal(sj ? d : 1.0 / d) : (sj ? 0.0 : PETSC_PI);
111: const PetscReal alpha = atanval / 3.0;
112: const PetscReal cd = s * PetscCosReal(alpha) * sjn;
113: const PetscReal sd = PetscSqrtReal(J2) * PetscSinReal(alpha);
115: x[0] = 2 * cd;
116: x[1] = -cd + sd;
117: x[2] = -cd - sd;
118: }
119: }
120: x[0] += I1 / 3.0;
121: x[1] += I1 / 3.0;
122: x[2] += I1 / 3.0;
123: }
124: }
125: #endif
127: /* compute shift to make C positive definite */
128: static inline PetscReal FIX_C_3d(PetscScalar C00, PetscScalar C01, PetscScalar C02, PetscScalar C11, PetscScalar C12, PetscScalar C22)
129: {
130: #if !PetscDefined(USE_COMPLEX)
131: PetscReal eigs[3], s = 0.0;
132: PetscBool twod = (PetscBool)(C02 == 0 && C12 == 0 && C22 == 0);
133: Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs);
134: if (twod) eigs[2] = 1.0;
135: if (eigs[0] <= 0 || eigs[1] <= 0 || eigs[2] <= 0) s = -PetscMin(eigs[0], PetscMin(eigs[1], eigs[2])) + PETSC_SMALL;
136: return s;
137: #else
138: return 0.0;
139: #endif
140: }
142: static inline PetscReal FIX_C(PetscScalar C00, PetscScalar C01, PetscScalar C11)
143: {
144: return FIX_C_3d(C00, C01, 0, C11, 0, 0);
145: }
147: /* residual for C when tested against basis functions */
148: 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[])
149: {
150: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
151: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
152: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
153: const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0);
154: const PetscScalar *gradp = split ? a_x + aOff_x[P_FIELD_ID] : u_x + uOff_x[P_FIELD_ID];
155: const PetscScalar outerp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[1] * gradp[1]};
156: const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]];
157: const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1];
158: const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2];
159: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
160: const PetscScalar nexp = (gamma - 2.0) / 2.0;
161: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
162: const PetscScalar eqss[] = {0.5, 1., 0.5}; /* equations rescaling for a symmetric Jacobian */
164: for (PetscInt k = 0; k < 3; k++) f0[k] = eqss[k] * (u_t[uOff[C_FIELD_ID] + k] - outerp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k]);
165: }
167: /* 3D version */
168: static void C_0_3d(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[])
169: {
170: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
171: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
172: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
173: const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0);
174: const PetscScalar *gradp = split ? a_x + aOff_x[P_FIELD_ID] : u_x + uOff_x[P_FIELD_ID];
175: const PetscScalar outerp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[0] * gradp[2], gradp[1] * gradp[1], gradp[1] * gradp[2], gradp[2] * gradp[2]};
176: const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]];
177: const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1];
178: const PetscScalar C02 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2];
179: const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 3] : u[uOff[C_FIELD_ID] + 3];
180: const PetscScalar C12 = split ? a[aOff[C_FIELD_ID] + 4] : u[uOff[C_FIELD_ID] + 4];
181: const PetscScalar C22 = split ? a[aOff[C_FIELD_ID] + 5] : u[uOff[C_FIELD_ID] + 5];
182: const PetscScalar norm = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps;
183: const PetscScalar nexp = (gamma - 2.0) / 2.0;
184: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
185: const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5};
187: for (PetscInt k = 0; k < 6; k++) f0[k] = eqss[k] * (u_t[uOff[C_FIELD_ID] + k] - outerp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k]);
188: }
190: /* Jacobian for C against C basis functions */
191: 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[])
192: {
193: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
194: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
195: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
196: const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0);
197: const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]];
198: const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1];
199: const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2];
200: const PetscScalar norm = NORM2C(C00, C01, C11) + eps;
201: const PetscScalar nexp = (gamma - 2.0) / 2.0;
202: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
203: const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0);
204: const PetscScalar dC[] = {2 * C00, 4 * C01, 2 * C11};
205: const PetscScalar eqss[] = {0.5, 1., 0.5};
207: for (PetscInt k = 0; k < 3; k++) {
208: if (!split) {
209: for (PetscInt j = 0; j < 3; j++) J[k * 3 + j] = eqss[k] * (alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k]);
210: }
211: J[k * 3 + k] += eqss[k] * (alpha * fnorm + u_tShift);
212: }
213: }
215: static void JC_0_c0c0_3d(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[])
216: {
217: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
218: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
219: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
220: const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0);
221: const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]];
222: const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1];
223: const PetscScalar C02 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2];
224: const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 3] : u[uOff[C_FIELD_ID] + 3];
225: const PetscScalar C12 = split ? a[aOff[C_FIELD_ID] + 4] : u[uOff[C_FIELD_ID] + 4];
226: const PetscScalar C22 = split ? a[aOff[C_FIELD_ID] + 5] : u[uOff[C_FIELD_ID] + 5];
227: const PetscScalar norm = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps;
228: const PetscScalar nexp = (gamma - 2.0) / 2.0;
229: const PetscScalar fnorm = PetscPowScalar(norm, nexp);
230: const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0);
231: const PetscScalar dC[] = {2 * C00, 4 * C01, 4 * C02, 2 * C11, 4 * C12, 2 * C22};
232: const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5};
234: for (PetscInt k = 0; k < 6; k++) {
235: if (!split) {
236: for (PetscInt j = 0; j < 6; j++) J[k * 6 + j] = eqss[k] * (alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k]);
237: }
238: J[k * 6 + k] += eqss[k] * (alpha * fnorm + u_tShift);
239: }
240: }
242: /* Jacobian for C against C basis functions and gradients of P basis functions */
243: 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[])
244: {
245: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
246: const PetscScalar eqss[] = {0.5, 1., 0.5};
248: J[0] = -2 * gradp[0] * eqss[0];
249: J[1] = 0.0;
251: J[2] = -gradp[1] * eqss[1];
252: J[3] = -gradp[0] * eqss[1];
254: J[4] = 0.0;
255: J[5] = -2 * gradp[1] * eqss[2];
256: }
258: static void JC_0_c0p1_3d(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[])
259: {
260: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
261: const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5};
263: J[0] = -2 * gradp[0] * eqss[0];
264: J[1] = 0.0;
265: J[2] = 0.0;
267: J[3] = -gradp[1] * eqss[1];
268: J[4] = -gradp[0] * eqss[1];
269: J[5] = 0.0;
271: J[6] = -gradp[2] * eqss[2];
272: J[7] = 0.0;
273: J[8] = -gradp[0] * eqss[2];
275: J[9] = 0.0;
276: J[10] = -2 * gradp[1] * eqss[3];
277: J[11] = 0.0;
279: J[12] = 0.0;
280: J[13] = -gradp[2] * eqss[4];
281: J[14] = -gradp[1] * eqss[4];
283: J[15] = 0.0;
284: J[16] = 0.0;
285: J[17] = -2 * gradp[2] * eqss[5];
286: }
288: /* residual for C when tested against gradients of basis functions */
289: 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[])
290: {
291: const PetscReal D = PetscRealPart(constants[D_ID]);
292: for (PetscInt k = 0; k < 3; k++)
293: for (PetscInt d = 0; d < 2; d++) f1[k * 2 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 2 + d];
294: }
296: static void C_1_3d(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[])
297: {
298: const PetscReal D = PetscRealPart(constants[D_ID]);
299: for (PetscInt k = 0; k < 6; k++)
300: for (PetscInt d = 0; d < 3; d++) f1[k * 3 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 3 + d];
301: }
303: /* Jacobian for C against gradients of C basis functions */
304: 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[])
305: {
306: const PetscReal D = PetscRealPart(constants[D_ID]);
307: for (PetscInt k = 0; k < 3; k++)
308: for (PetscInt d = 0; d < 2; d++) J[k * (3 + 1) * 2 * 2 + d * 2 + d] = PetscSqr(D);
309: }
311: static void JC_1_c1c1_3d(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[])
312: {
313: const PetscReal D = PetscRealPart(constants[D_ID]);
314: for (PetscInt k = 0; k < 6; k++)
315: for (PetscInt d = 0; d < 3; d++) J[k * (6 + 1) * 3 * 3 + d * 3 + d] = PetscSqr(D);
316: }
318: /* residual for P when tested against basis functions.
319: The source term always comes from the auxiliary data because it must be zero mean (algebraically) */
320: 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[])
321: {
322: PetscScalar S = a[aOff[NUM_FIELDS]];
324: f0[0] = S;
325: }
327: /* residual for P when tested against basis functions for the initial condition problem
328: here we don't impose symmetry, and we thus flip the sign of the source function */
329: static void P_0_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 f0[])
330: {
331: PetscScalar S = a[aOff[NUM_FIELDS]];
333: f0[0] = -S;
334: }
336: /* residual for P when tested against gradients of basis functions */
337: 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[])
338: {
339: const PetscReal r = PetscRealPart(constants[R_ID]);
340: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
341: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
342: const PetscScalar C10 = C01;
343: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r;
344: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
345: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
346: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
348: f1[0] = -((C00 + s) * gradp[0] + C01 * gradp[1]);
349: f1[1] = -(C10 * gradp[0] + (C11 + s) * gradp[1]);
350: }
352: static void P_1_3d(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[])
353: {
354: const PetscReal r = PetscRealPart(constants[R_ID]);
355: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
356: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
357: const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2];
358: const PetscScalar C10 = C01;
359: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r;
360: const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4];
361: const PetscScalar C20 = C02;
362: const PetscScalar C21 = C12;
363: const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r;
364: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
365: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
366: const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0;
368: f1[0] = -((C00 + s) * gradp[0] + C01 * gradp[1] + C02 * gradp[2]);
369: f1[1] = -(C10 * gradp[0] + (C11 + s) * gradp[1] + C12 * gradp[2]);
370: f1[2] = -(C20 * gradp[0] + C21 * gradp[1] + (C22 + s) * gradp[2]);
371: }
373: /* Same as above except that the conductivity values come from the auxiliary vec */
374: 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[])
375: {
376: const PetscReal r = PetscRealPart(constants[R_ID]);
377: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
378: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
379: const PetscScalar C10 = C01;
380: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r;
381: const PetscScalar *gradp = u_x + uOff_x[Nf > 1 ? P_FIELD_ID : 0];
382: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
383: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
385: f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1];
386: f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1];
387: }
389: static void P_1_aux_3d(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[])
390: {
391: const PetscReal r = PetscRealPart(constants[R_ID]);
392: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
393: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
394: const PetscScalar C02 = a[aOff[C_FIELD_ID] + 2];
395: const PetscScalar C10 = C01;
396: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 3] + r;
397: const PetscScalar C12 = a[aOff[C_FIELD_ID] + 4];
398: const PetscScalar C20 = C02;
399: const PetscScalar C21 = C12;
400: const PetscScalar C22 = a[aOff[C_FIELD_ID] + 5] + r;
401: const PetscScalar *gradp = u_x + uOff_x[Nf > 1 ? P_FIELD_ID : 0];
402: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
403: const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0;
405: f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1] + C02 * gradp[2];
406: f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1] + C12 * gradp[2];
407: f1[2] = C20 * gradp[0] + C21 * gradp[1] + (C22 + s) * gradp[2];
408: }
410: /* Jacobian for P against gradients of P basis functions */
411: 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[])
412: {
413: const PetscReal r = PetscRealPart(constants[R_ID]);
414: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
415: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
416: const PetscScalar C10 = C01;
417: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r;
418: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
419: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
421: J[0] = -(C00 + s);
422: J[1] = -C01;
423: J[2] = -C10;
424: J[3] = -(C11 + s);
425: }
427: static void JP_1_p1p1_3d(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[])
428: {
429: const PetscReal r = PetscRealPart(constants[R_ID]);
430: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
431: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
432: const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2];
433: const PetscScalar C10 = C01;
434: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r;
435: const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4];
436: const PetscScalar C20 = C02;
437: const PetscScalar C21 = C12;
438: const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r;
439: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
440: const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0;
442: J[0] = -(C00 + s);
443: J[1] = -C01;
444: J[2] = -C02;
445: J[3] = -C10;
446: J[4] = -(C11 + s);
447: J[5] = -C12;
448: J[6] = -C20;
449: J[7] = -C21;
450: J[8] = -(C22 + s);
451: }
453: 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[])
454: {
455: const PetscReal r = PetscRealPart(constants[R_ID]);
456: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
457: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
458: const PetscScalar C10 = C01;
459: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r;
460: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
461: const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0;
463: J[0] = C00 + s;
464: J[1] = C01;
465: J[2] = C10;
466: J[3] = C11 + s;
467: }
469: static void JP_1_p1p1_aux_3d(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[])
470: {
471: const PetscReal r = PetscRealPart(constants[R_ID]);
472: const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r;
473: const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1];
474: const PetscScalar C02 = a[aOff[C_FIELD_ID] + 2];
475: const PetscScalar C10 = C01;
476: const PetscScalar C11 = a[aOff[C_FIELD_ID] + 3] + r;
477: const PetscScalar C12 = a[aOff[C_FIELD_ID] + 4];
478: const PetscScalar C20 = C02;
479: const PetscScalar C21 = C12;
480: const PetscScalar C22 = a[aOff[C_FIELD_ID] + 5] + r;
481: const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
482: const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0;
484: J[0] = C00 + s;
485: J[1] = C01;
486: J[2] = C02;
487: J[3] = C10;
488: J[4] = C11 + s;
489: J[5] = C12;
490: J[6] = C20;
491: J[7] = C21;
492: J[8] = C22 + s;
493: }
495: /* Jacobian for P against gradients of P basis functions and C basis functions */
496: 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[])
497: {
498: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
500: J[0] = -gradp[0];
501: J[1] = 0;
503: J[2] = -gradp[1];
504: J[3] = -gradp[0];
506: J[4] = 0;
507: J[5] = -gradp[1];
508: }
510: static void JP_1_p1c0_3d(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[])
511: {
512: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
514: J[0] = -gradp[0];
515: J[1] = 0;
516: J[2] = 0;
518: J[3] = -gradp[1];
519: J[4] = -gradp[0];
520: J[5] = 0;
522: J[6] = -gradp[2];
523: J[7] = 0;
524: J[8] = -gradp[0];
526: J[9] = 0;
527: J[10] = -gradp[1];
528: J[11] = 0;
530: J[12] = 0;
531: J[13] = -gradp[2];
532: J[14] = -gradp[1];
534: J[15] = 0;
535: J[16] = 0;
536: J[17] = -gradp[2];
537: }
539: /* a collection of gaussian, Dirac-like, source term S(x) = scale * cos(2*pi*(frequency*t + phase)) * exp(-w*||x - x0||^2) */
540: typedef struct {
541: PetscInt n;
542: PetscReal *x0;
543: PetscReal *w;
544: PetscReal *k;
545: PetscReal *p;
546: PetscReal *r;
547: } MultiSourceCtx;
549: typedef struct {
550: PetscReal x0[3];
551: PetscReal w;
552: PetscReal k;
553: PetscReal p;
554: PetscReal r;
555: } SingleSourceCtx;
557: static PetscErrorCode gaussian(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
558: {
559: SingleSourceCtx *s = (SingleSourceCtx *)ctx;
560: const PetscReal *x0 = s->x0;
561: const PetscReal w = s->w;
562: const PetscReal k = s->k; /* frequency */
563: const PetscReal p = s->p; /* phase */
564: const PetscReal r = s->r; /* scale */
565: PetscReal n = 0;
567: for (PetscInt d = 0; d < dim; ++d) n += (x[d] - x0[d]) * (x[d] - x0[d]);
568: u[0] = r * PetscCosReal(2 * PETSC_PI * (k * time + p)) * PetscExpReal(-w * n);
569: return PETSC_SUCCESS;
570: }
572: static PetscErrorCode source(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
573: {
574: MultiSourceCtx *s = (MultiSourceCtx *)ctx;
576: u[0] = 0.0;
577: for (PetscInt i = 0; i < s->n; i++) {
578: SingleSourceCtx sctx;
579: PetscScalar ut[1];
581: sctx.x0[0] = s->x0[dim * i];
582: sctx.x0[1] = s->x0[dim * i + 1];
583: sctx.x0[2] = dim > 2 ? s->x0[dim * i + 2] : 0.0;
584: sctx.w = s->w[i];
585: sctx.k = s->k[i];
586: sctx.p = s->p[i];
587: sctx.r = s->r[i];
589: PetscCall(gaussian(dim, time, x, Nf, ut, &sctx));
591: u[0] += ut[0];
592: }
593: return PETSC_SUCCESS;
594: }
596: /* functionals to be integrated: average -> \int_\Omega u dx */
597: 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[])
598: {
599: const PetscInt fid = (PetscInt)PetscRealPart(constants[numConstants]);
600: obj[0] = u[uOff[fid]];
601: }
603: /* functionals to be integrated: volume -> \int_\Omega dx */
604: static void volume(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[])
605: {
606: obj[0] = 1;
607: }
609: /* functionals to be integrated: energy -> D^2/2 * ||\nabla C||^2 + c^2\nabla p * (r + C) * \nabla p + \alpha/ \gamma * ||C||^\gamma */
610: 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[])
611: {
612: const PetscReal D = PetscRealPart(constants[D_ID]);
613: const PetscReal r = PetscRealPart(constants[R_ID]);
614: const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]);
615: const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]);
616: const PetscReal eps = PetscRealPart(constants[EPS_ID]);
618: if (dim == 2) {
619: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
620: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
621: const PetscScalar C10 = C01;
622: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
623: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
624: const PetscScalar *gradC00 = u_x + uOff_x[C_FIELD_ID];
625: const PetscScalar *gradC01 = u_x + uOff_x[C_FIELD_ID] + 2;
626: const PetscScalar *gradC11 = u_x + uOff_x[C_FIELD_ID] + 4;
627: const PetscScalar normC = NORM2C(C00, C01, C11) + eps;
628: const PetscScalar normgradC = NORM2C(gradC00[0], gradC01[0], gradC11[0]) + NORM2C(gradC00[1], gradC01[1], gradC11[1]);
629: const PetscScalar nexp = gamma / 2.0;
631: const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC;
632: const PetscScalar t1 = gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1]);
633: const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp);
635: obj[0] = t0 + t1 + t2;
636: } else {
637: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
638: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
639: const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2];
640: const PetscScalar C10 = C01;
641: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3];
642: const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4];
643: const PetscScalar C20 = C02;
644: const PetscScalar C21 = C12;
645: const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5];
646: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
647: const PetscScalar *gradC00 = u_x + uOff_x[C_FIELD_ID];
648: const PetscScalar *gradC01 = u_x + uOff_x[C_FIELD_ID] + 3;
649: const PetscScalar *gradC02 = u_x + uOff_x[C_FIELD_ID] + 6;
650: const PetscScalar *gradC11 = u_x + uOff_x[C_FIELD_ID] + 9;
651: const PetscScalar *gradC12 = u_x + uOff_x[C_FIELD_ID] + 12;
652: const PetscScalar *gradC22 = u_x + uOff_x[C_FIELD_ID] + 15;
653: const PetscScalar normC = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps;
654: const PetscScalar normgradC = NORM2C_3d(gradC00[0], gradC01[0], gradC02[0], gradC11[0], gradC12[0], gradC22[0]) + NORM2C_3d(gradC00[1], gradC01[1], gradC02[1], gradC11[1], gradC12[1], gradC22[1]) + NORM2C_3d(gradC00[2], gradC01[2], gradC02[2], gradC11[2], gradC12[2], gradC22[2]);
655: const PetscScalar nexp = gamma / 2.0;
657: const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC;
658: const PetscScalar t1 = gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1] + C02 * gradp[2]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1] + C12 * gradp[2]) + gradp[2] * (C20 * gradp[0] + C21 * gradp[1] + (C22 + r) * gradp[2]);
659: const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp);
661: obj[0] = t0 + t1 + t2;
662: }
663: }
665: /* functionals to be integrated: ellipticity_fail_private -> see below */
666: static void ellipticity_fail_private(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[], PetscBool add_reg)
667: {
668: #if !PetscDefined(USE_COMPLEX)
669: PetscReal eigs[3];
670: PetscScalar C00, C01, C02 = 0.0, C11, C12 = 0.0, C22 = 0.0;
671: const PetscReal r = add_reg ? PetscRealPart(constants[R_ID]) : 0.0;
673: if (dim == 2) {
674: C00 = u[uOff[C_FIELD_ID]] + r;
675: C01 = u[uOff[C_FIELD_ID] + 1];
676: C11 = u[uOff[C_FIELD_ID] + 2] + r;
677: } else {
678: C00 = u[uOff[C_FIELD_ID]] + r;
679: C01 = u[uOff[C_FIELD_ID] + 1];
680: C02 = u[uOff[C_FIELD_ID] + 2];
681: C11 = u[uOff[C_FIELD_ID] + 3] + r;
682: C12 = u[uOff[C_FIELD_ID] + 4];
683: C22 = u[uOff[C_FIELD_ID] + 5] + r;
684: }
685: Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs);
686: if (eigs[0] < 0 || eigs[1] < 0 || eigs[2] < 0) obj[0] = -PetscMin(eigs[0], PetscMin(eigs[1], eigs[2]));
687: else obj[0] = 0.0;
688: #else
689: obj[0] = 0.0;
690: #endif
691: }
693: /* functionals to be integrated: ellipticity_fail -> 0 means C is elliptic at quadrature point, otherwise it returns 1 */
694: 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[])
695: {
696: ellipticity_fail_private(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, obj, PETSC_FALSE);
697: if (PetscAbsScalar(obj[0]) > 0.0) obj[0] = 1.0;
698: }
700: /* functionals to be integrated: jacobian_fail -> 0 means C + r is elliptic at quadrature point, otherwise it returns the absolute value of the most negative eigenvalue */
701: static void jacobian_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[])
702: {
703: ellipticity_fail_private(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, obj, PETSC_TRUE);
704: }
706: /* initial conditions for C: eq. 16 */
707: static PetscErrorCode initial_conditions_C_0(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
708: {
709: if (dim == 2) {
710: u[0] = 1;
711: u[1] = 0;
712: u[2] = 1;
713: } else {
714: u[0] = 1;
715: u[1] = 0;
716: u[2] = 0;
717: u[3] = 1;
718: u[4] = 0;
719: u[5] = 1;
720: }
721: return PETSC_SUCCESS;
722: }
724: /* initial conditions for C: eq. 17 */
725: static PetscErrorCode initial_conditions_C_1(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
726: {
727: const PetscReal x = xx[0];
728: const PetscReal y = xx[1];
730: if (dim == 3) return PETSC_ERR_SUP;
731: u[0] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
732: u[1] = 0;
733: u[2] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
734: return PETSC_SUCCESS;
735: }
737: /* initial conditions for C: eq. 18 */
738: static PetscErrorCode initial_conditions_C_2(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
739: {
740: u[0] = 0;
741: u[1] = 0;
742: u[2] = 0;
743: if (dim == 3) {
744: u[3] = 0;
745: u[4] = 0;
746: u[5] = 0;
747: }
748: return PETSC_SUCCESS;
749: }
751: /* random initial conditions for the diagonal of C */
752: static PetscErrorCode initial_conditions_C_random(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
753: {
754: PetscScalar vals[3];
755: PetscRandom r = (PetscRandom)ctx;
757: PetscCall(PetscRandomGetValues(r, dim, vals));
758: if (dim == 2) {
759: u[0] = vals[0];
760: u[1] = 0;
761: u[2] = vals[1];
762: } else {
763: u[0] = vals[0];
764: u[1] = 0;
765: u[2] = 0;
766: u[3] = vals[1];
767: u[4] = 0;
768: u[5] = vals[2];
769: }
770: return PETSC_SUCCESS;
771: }
773: /* functionals to be sampled: zero */
774: 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[])
775: {
776: f[0] = 0.0;
777: }
779: /* functionals to be sampled: - (C + r) * \grad p */
780: 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[])
781: {
782: const PetscReal r = PetscRealPart(constants[R_ID]);
783: const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID];
785: if (dim == 2) {
786: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
787: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
788: const PetscScalar C10 = C01;
789: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r;
791: f[0] = -C00 * gradp[0] - C01 * gradp[1];
792: f[1] = -C10 * gradp[0] - C11 * gradp[1];
793: } else {
794: const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r;
795: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
796: const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2];
797: const PetscScalar C10 = C01;
798: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r;
799: const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4];
800: const PetscScalar C20 = C02;
801: const PetscScalar C21 = C12;
802: const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r;
804: f[0] = -C00 * gradp[0] - C01 * gradp[1] - C02 * gradp[2];
805: f[1] = -C10 * gradp[0] - C11 * gradp[1] - C12 * gradp[2];
806: f[2] = -C20 * gradp[0] - C21 * gradp[1] - C22 * gradp[2];
807: }
808: }
810: /* functionals to be sampled: ||C|| */
811: 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[])
812: {
813: if (dim == 2) {
814: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
815: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
816: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];
818: f[0] = PetscSqrtReal(PetscRealPart(NORM2C(C00, C01, C11)));
819: } else {
820: const PetscScalar C00 = u[uOff[C_FIELD_ID]];
821: const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
822: const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2];
823: const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3];
824: const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4];
825: const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5];
827: f[0] = PetscSqrtReal(PetscRealPart(NORM2C_3d(C00, C01, C02, C11, C12, C22)));
828: }
829: }
831: /* functionals to be sampled: eigenvalues of C */
832: static void eigsc(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[])
833: {
834: #if !PetscDefined(USE_COMPLEX)
835: PetscReal eigs[3];
836: PetscScalar C00, C01, C02 = 0.0, C11, C12 = 0.0, C22 = 0.0;
837: if (dim == 2) {
838: C00 = u[uOff[C_FIELD_ID]];
839: C01 = u[uOff[C_FIELD_ID] + 1];
840: C11 = u[uOff[C_FIELD_ID] + 2];
841: } else {
842: C00 = u[uOff[C_FIELD_ID]];
843: C01 = u[uOff[C_FIELD_ID] + 1];
844: C02 = u[uOff[C_FIELD_ID] + 2];
845: C11 = u[uOff[C_FIELD_ID] + 3];
846: C12 = u[uOff[C_FIELD_ID] + 4];
847: C22 = u[uOff[C_FIELD_ID] + 5];
848: }
849: Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs);
850: PetscCallVoid(PetscSortReal(dim, eigs));
851: for (PetscInt d = 0; d < dim; d++) f[d] = eigs[d];
852: #else
853: for (PetscInt d = 0; d < dim; d++) f[d] = 0;
854: #endif
855: }
857: /* functions to be sampled: zero function */
858: static PetscErrorCode zerof(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
859: {
860: for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
861: return PETSC_SUCCESS;
862: }
864: /* functions to be sampled: constant function */
865: static PetscErrorCode constantf(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
866: {
867: for (PetscInt d = 0; d < Nc; ++d) u[d] = 1.0;
868: return PETSC_SUCCESS;
869: }
871: /* application context: customizable parameters */
872: typedef struct {
873: PetscInt dim;
874: PetscBool simplex;
875: PetscReal r;
876: PetscReal eps;
877: PetscReal alpha;
878: PetscReal gamma;
879: PetscReal D;
880: PetscReal domain_volume;
881: PetscInt ic_num;
882: PetscBool split;
883: PetscBool lump;
884: PetscBool amr;
885: PetscBool load;
886: char load_filename[PETSC_MAX_PATH_LEN];
887: PetscBool save;
888: char save_filename[PETSC_MAX_PATH_LEN];
889: PetscInt save_every;
890: PetscBool test_restart;
891: PetscInt fix_c;
892: MultiSourceCtx *source_ctx;
893: DM view_dm;
894: TSMonitorVTKCtx view_vtk_ctx;
895: PetscViewerAndFormat *view_hdf5_ctx;
896: PetscInt diagnostic_num;
897: PetscReal view_times[64];
898: PetscInt view_times_n, view_times_k;
899: PetscReal function_domain_error_tol;
900: VecScatter subsct[NUM_FIELDS];
901: Vec subvec[NUM_FIELDS];
902: PetscBool monitor_norms;
903: PetscBool exclude_potential_lte;
905: /* hack: need some more plumbing in the library */
906: SNES snes;
907: } AppCtx;
909: /* process command line options */
910: #include <petsc/private/tsimpl.h>
911: static PetscErrorCode ProcessOptions(AppCtx *options)
912: {
913: char vtkmonfilename[PETSC_MAX_PATH_LEN];
914: char hdf5monfilename[PETSC_MAX_PATH_LEN];
915: PetscInt tmp;
916: PetscBool flg;
918: PetscFunctionBeginUser;
919: options->dim = 2;
920: options->r = 1.e-1;
921: options->eps = 1.e-3;
922: options->alpha = 0.75;
923: options->gamma = 0.75;
924: options->D = 1.e-2;
925: options->ic_num = 0;
926: options->split = PETSC_FALSE;
927: options->lump = PETSC_FALSE;
928: options->amr = PETSC_FALSE;
929: options->load = PETSC_FALSE;
930: options->save = PETSC_FALSE;
931: options->save_every = -1;
932: options->test_restart = PETSC_FALSE;
933: options->fix_c = 1; /* 1 means only Jac, 2 means function and Jac, < 0 means raise SNESFunctionDomainError when C+r is not posdef */
934: options->view_vtk_ctx = NULL;
935: options->view_hdf5_ctx = NULL;
936: options->view_dm = NULL;
937: options->diagnostic_num = 1;
938: options->function_domain_error_tol = -1;
939: options->monitor_norms = PETSC_FALSE;
940: options->exclude_potential_lte = PETSC_FALSE;
941: for (PetscInt i = 0; i < NUM_FIELDS; i++) {
942: options->subsct[i] = NULL;
943: options->subvec[i] = NULL;
944: }
945: for (PetscInt i = 0; i < 64; i++) options->view_times[i] = PETSC_MAX_REAL;
947: PetscOptionsBegin(PETSC_COMM_WORLD, "", __FILE__, "DMPLEX");
948: PetscCall(PetscOptionsInt("-dim", "space dimension", __FILE__, options->dim, &options->dim, NULL));
949: PetscCall(PetscOptionsReal("-alpha", "alpha", __FILE__, options->alpha, &options->alpha, NULL));
950: PetscCall(PetscOptionsReal("-gamma", "gamma", __FILE__, options->gamma, &options->gamma, NULL));
951: PetscCall(PetscOptionsReal("-d", "D", __FILE__, options->D, &options->D, NULL));
952: PetscCall(PetscOptionsReal("-eps", "eps", __FILE__, options->eps, &options->eps, NULL));
953: PetscCall(PetscOptionsReal("-r", "r", __FILE__, options->r, &options->r, NULL));
954: PetscCall(PetscOptionsInt("-ic_num", "ic_num", __FILE__, options->ic_num, &options->ic_num, NULL));
955: PetscCall(PetscOptionsBool("-split", "Operator splitting", __FILE__, options->split, &options->split, NULL));
956: PetscCall(PetscOptionsBool("-lump", "use mass lumping", __FILE__, options->lump, &options->lump, NULL));
957: PetscCall(PetscOptionsInt("-fix_c", "Fix conductivity: shift to always be positive semi-definite or raise domain error", __FILE__, options->fix_c, &options->fix_c, NULL));
958: PetscCall(PetscOptionsBool("-amr", "use adaptive mesh refinement", __FILE__, options->amr, &options->amr, NULL));
959: PetscCall(PetscOptionsReal("-domain_error_tol", "domain error tolerance", __FILE__, options->function_domain_error_tol, &options->function_domain_error_tol, NULL));
960: PetscCall(PetscOptionsBool("-test_restart", "test restarting files", __FILE__, options->test_restart, &options->test_restart, NULL));
961: if (!options->test_restart) {
962: PetscCall(PetscOptionsString("-load", "filename with data to be loaded for restarting", __FILE__, options->load_filename, options->load_filename, PETSC_MAX_PATH_LEN, &options->load));
963: PetscCall(PetscOptionsString("-save", "filename with data to be saved for restarting", __FILE__, options->save_filename, options->save_filename, PETSC_MAX_PATH_LEN, &options->save));
964: if (options->save) PetscCall(PetscOptionsInt("-save_every", "save every n timestep (-1 saves only the last)", __FILE__, options->save_every, &options->save_every, NULL));
965: }
966: PetscCall(PetscOptionsBool("-exclude_potential_lte", "exclude potential from LTE", __FILE__, options->exclude_potential_lte, &options->exclude_potential_lte, NULL));
967: options->view_times_k = 0;
968: options->view_times_n = 0;
969: PetscCall(PetscOptionsRealArray("-monitor_times", "Save at specific times", NULL, options->view_times, (tmp = 64, &tmp), &flg));
970: if (flg) options->view_times_n = tmp;
972: PetscCall(PetscOptionsString("-monitor_vtk", "Dump VTK file for diagnostic", "TSMonitorSolutionVTK", NULL, vtkmonfilename, sizeof(vtkmonfilename), &flg));
973: if (flg) {
974: PetscCall(TSMonitorSolutionVTKCtxCreate(vtkmonfilename, &options->view_vtk_ctx));
975: PetscCall(PetscOptionsInt("-monitor_vtk_interval", "Save every interval time steps", NULL, options->view_vtk_ctx->interval, &options->view_vtk_ctx->interval, NULL));
976: }
977: PetscCall(PetscOptionsString("-monitor_hdf5", "Dump HDF5 file for diagnostic", "TSMonitorSolution", NULL, hdf5monfilename, sizeof(hdf5monfilename), &flg));
978: PetscCall(PetscOptionsInt("-monitor_diagnostic_num", "Number of diagnostics to be computed", __FILE__, options->diagnostic_num, &options->diagnostic_num, NULL));
980: if (flg) {
981: #if defined(PETSC_HAVE_HDF5)
982: PetscViewer viewer;
984: PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, hdf5monfilename, FILE_MODE_WRITE, &viewer));
985: PetscCall(PetscViewerAndFormatCreate(viewer, PETSC_VIEWER_HDF5_VIZ, &options->view_hdf5_ctx));
986: options->view_hdf5_ctx->view_interval = 1;
987: PetscCall(PetscOptionsInt("-monitor_hdf5_interval", "Save every interval time steps", NULL, options->view_hdf5_ctx->view_interval, &options->view_hdf5_ctx->view_interval, NULL));
988: PetscCall(PetscViewerDestroy(&viewer));
989: #else
990: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
991: #endif
992: }
993: PetscCall(PetscOptionsBool("-monitor_norms", "Monitor separate SNES norms", __FILE__, options->monitor_norms, &options->monitor_norms, NULL));
995: /* source options */
996: PetscCall(PetscNew(&options->source_ctx));
997: options->source_ctx->n = 1;
999: PetscCall(PetscOptionsInt("-source_num", "number of sources", __FILE__, options->source_ctx->n, &options->source_ctx->n, NULL));
1000: tmp = options->source_ctx->n;
1001: PetscCall(PetscMalloc5(options->dim * tmp, &options->source_ctx->x0, tmp, &options->source_ctx->w, tmp, &options->source_ctx->k, tmp, &options->source_ctx->p, tmp, &options->source_ctx->r));
1002: for (PetscInt i = 0; i < options->source_ctx->n; i++) {
1003: for (PetscInt d = 0; d < options->dim; d++) options->source_ctx->x0[options->dim * i + d] = 0.25;
1004: options->source_ctx->w[i] = 500;
1005: options->source_ctx->k[i] = 0;
1006: options->source_ctx->p[i] = 0;
1007: options->source_ctx->r[i] = 1;
1008: }
1009: tmp = options->dim * options->source_ctx->n;
1010: PetscCall(PetscOptionsRealArray("-source_x0", "source location", __FILE__, options->source_ctx->x0, &tmp, NULL));
1011: tmp = options->source_ctx->n;
1012: PetscCall(PetscOptionsRealArray("-source_w", "source factor", __FILE__, options->source_ctx->w, &tmp, NULL));
1013: tmp = options->source_ctx->n;
1014: PetscCall(PetscOptionsRealArray("-source_k", "source frequency", __FILE__, options->source_ctx->k, &tmp, NULL));
1015: tmp = options->source_ctx->n;
1016: PetscCall(PetscOptionsRealArray("-source_p", "source phase", __FILE__, options->source_ctx->p, &tmp, NULL));
1017: tmp = options->source_ctx->n;
1018: PetscCall(PetscOptionsRealArray("-source_r", "source scaling", __FILE__, options->source_ctx->r, &tmp, NULL));
1019: PetscOptionsEnd();
1020: PetscFunctionReturn(PETSC_SUCCESS);
1021: }
1023: static PetscErrorCode SaveToFile(DM dm, Vec u, const char *filename)
1024: {
1025: #if defined(PETSC_HAVE_HDF5)
1026: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
1027: PetscViewer viewer;
1028: DM cdm = dm;
1029: PetscInt numlevels = 0;
1031: PetscFunctionBeginUser;
1032: while (cdm) {
1033: numlevels++;
1034: PetscCall(DMGetCoarseDM(cdm, &cdm));
1035: }
1036: /* Cannot be set programmatically */
1037: PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_view_hdf5_storage_version 3.0.0"));
1038: PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), filename, FILE_MODE_WRITE, &viewer));
1039: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numlevels", PETSC_INT, &numlevels));
1040: PetscCall(PetscViewerPushFormat(viewer, format));
1041: for (PetscInt level = numlevels - 1; level >= 0; level--) {
1042: PetscInt cc, rr;
1043: PetscBool isRegular, isUniform;
1044: const char *dmname;
1045: char groupname[PETSC_MAX_PATH_LEN];
1047: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
1048: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
1049: PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
1050: PetscCall(DMGetCoarsenLevel(dm, &cc));
1051: PetscCall(DMGetRefineLevel(dm, &rr));
1052: PetscCall(DMPlexGetRegularRefinement(dm, &isRegular));
1053: PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
1054: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "meshname", PETSC_STRING, dmname));
1055: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refinelevel", PETSC_INT, &rr));
1056: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, &cc));
1057: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refRegular", PETSC_BOOL, &isRegular));
1058: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refUniform", PETSC_BOOL, &isUniform));
1059: PetscCall(DMPlexTopologyView(dm, viewer));
1060: PetscCall(DMPlexLabelsView(dm, viewer));
1061: PetscCall(DMPlexCoordinatesView(dm, viewer));
1062: PetscCall(DMPlexSectionView(dm, viewer, NULL));
1063: if (level == numlevels - 1) {
1064: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
1065: PetscCall(DMPlexGlobalVectorView(dm, viewer, NULL, u));
1066: }
1067: if (level) {
1068: PetscInt cStart, cEnd, ccStart, ccEnd, cpStart;
1069: DMPolytopeType ct;
1070: DMPlexTransform tr;
1071: DM sdm;
1072: PetscScalar *array;
1073: PetscSection section;
1074: Vec map;
1075: IS gis;
1076: const PetscInt *gidx;
1078: PetscCall(DMGetCoarseDM(dm, &cdm));
1079: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1080: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
1081: PetscCall(PetscSectionSetChart(section, cStart, cEnd));
1082: for (PetscInt c = cStart; c < cEnd; c++) PetscCall(PetscSectionSetDof(section, c, 1));
1083: PetscCall(PetscSectionSetUp(section));
1085: PetscCall(DMClone(dm, &sdm));
1086: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
1087: PetscCall(PetscObjectSetName((PetscObject)section, "pdm_section"));
1088: PetscCall(DMSetLocalSection(sdm, section));
1089: PetscCall(PetscSectionDestroy(§ion));
1091: PetscCall(DMGetLocalVector(sdm, &map));
1092: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
1093: PetscCall(VecGetArray(map, &array));
1094: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1095: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1096: PetscCall(DMPlexTransformSetDM(tr, cdm));
1097: PetscCall(DMPlexTransformSetFromOptions(tr));
1098: PetscCall(DMPlexTransformSetUp(tr));
1099: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
1100: PetscCall(DMPlexGetChart(cdm, &cpStart, NULL));
1101: PetscCall(DMPlexCreatePointNumbering(cdm, &gis));
1102: PetscCall(ISGetIndices(gis, &gidx));
1103: for (PetscInt c = ccStart; c < ccEnd; c++) {
1104: PetscInt *rsize, *rcone, *rornt, Nt;
1105: DMPolytopeType *rct;
1106: PetscInt gnum = gidx[c - cpStart] >= 0 ? gidx[c - cpStart] : -(gidx[c - cpStart] + 1);
1108: PetscCall(DMPlexGetCellType(cdm, c, &ct));
1109: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nt, &rct, &rsize, &rcone, &rornt));
1110: for (PetscInt r = 0; r < rsize[Nt - 1]; ++r) {
1111: PetscInt pNew;
1113: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[Nt - 1], c, r, &pNew));
1114: array[pNew - cStart] = gnum;
1115: }
1116: }
1117: PetscCall(ISRestoreIndices(gis, &gidx));
1118: PetscCall(ISDestroy(&gis));
1119: PetscCall(VecRestoreArray(map, &array));
1120: PetscCall(DMPlexTransformDestroy(&tr));
1121: PetscCall(DMPlexSectionView(dm, viewer, sdm));
1122: PetscCall(DMPlexLocalVectorView(dm, viewer, sdm, map));
1123: PetscCall(DMRestoreLocalVector(sdm, &map));
1124: PetscCall(DMDestroy(&sdm));
1125: }
1126: PetscCall(PetscViewerHDF5PopGroup(viewer));
1127: PetscCall(DMGetCoarseDM(dm, &dm));
1128: }
1129: PetscCall(PetscViewerPopFormat(viewer));
1130: PetscCall(PetscViewerDestroy(&viewer));
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: #else
1133: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
1134: #endif
1135: }
1137: static PetscErrorCode LoadFromFile(MPI_Comm comm, const char *filename, DM *odm)
1138: {
1139: #if defined(PETSC_HAVE_HDF5)
1140: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
1141: PetscViewer viewer;
1142: DM dm, cdm = NULL;
1143: PetscSF sfXC = NULL;
1144: PetscInt numlevels = -1;
1146: PetscFunctionBeginUser;
1147: PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer));
1148: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numlevels", PETSC_INT, NULL, &numlevels));
1149: PetscCall(PetscViewerPushFormat(viewer, format));
1150: for (PetscInt level = 0; level < numlevels; level++) {
1151: char groupname[PETSC_MAX_PATH_LEN], *dmname;
1152: PetscSF sfXB, sfBC, sfG;
1153: PetscPartitioner part;
1154: PetscInt rr, cc;
1155: PetscBool isRegular, isUniform;
1157: PetscCall(DMCreate(comm, &dm));
1158: PetscCall(DMSetType(dm, DMPLEX));
1159: PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
1160: PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
1161: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "meshname", PETSC_STRING, NULL, &dmname));
1162: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refinelevel", PETSC_INT, NULL, &rr));
1163: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, NULL, &cc));
1164: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refRegular", PETSC_BOOL, NULL, &isRegular));
1165: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refUniform", PETSC_BOOL, NULL, &isUniform));
1166: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
1167: PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
1168: PetscCall(DMPlexLabelsLoad(dm, viewer, sfXB));
1169: PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXB));
1170: PetscCall(DMPlexGetPartitioner(dm, &part));
1171: if (!level) { /* partition the coarse level only */
1172: PetscCall(PetscPartitionerSetFromOptions(part));
1173: } else { /* propagate partitioning information from coarser to finer level */
1174: DM sdm;
1175: Vec map;
1176: PetscSF sf;
1177: PetscLayout clayout;
1178: PetscScalar *array;
1179: PetscInt *cranks_leaf, *cranks_root, *npoints, *points, *ranks, *starts, *gidxs;
1180: PetscInt nparts, cStart, cEnd, nr, ccStart, ccEnd, cpStart, cpEnd;
1181: PetscMPIInt size, rank;
1183: PetscCall(DMClone(dm, &sdm));
1184: PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
1185: PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXB, NULL, &sf));
1186: PetscCall(DMGetLocalVector(sdm, &map));
1187: PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
1188: PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, sf, map));
1190: PetscCallMPI(MPI_Comm_size(comm, &size));
1191: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1192: nparts = size;
1193: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1194: PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
1195: PetscCall(DMPlexGetChart(cdm, &cpStart, &cpEnd));
1196: PetscCall(PetscCalloc1(nparts, &npoints));
1197: PetscCall(PetscMalloc4(cEnd - cStart, &points, cEnd - cStart, &ranks, nparts + 1, &starts, cEnd - cStart, &gidxs));
1198: PetscCall(PetscSFGetGraph(sfXC, &nr, NULL, NULL, NULL));
1199: PetscCall(PetscMalloc2(cpEnd - cpStart, &cranks_leaf, nr, &cranks_root));
1200: for (PetscInt c = 0; c < cpEnd - cpStart; c++) cranks_leaf[c] = rank;
1201: PetscCall(PetscSFReduceBegin(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
1202: PetscCall(PetscSFReduceEnd(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
1204: PetscCall(VecGetArray(map, &array));
1205: for (PetscInt c = 0; c < cEnd - cStart; c++) gidxs[c] = (PetscInt)PetscRealPart(array[c]);
1206: PetscCall(VecRestoreArray(map, &array));
1208: PetscCall(PetscLayoutCreate(comm, &clayout));
1209: PetscCall(PetscLayoutSetLocalSize(clayout, nr));
1210: PetscCall(PetscSFSetGraphLayout(sf, clayout, cEnd - cStart, NULL, PETSC_OWN_POINTER, gidxs));
1211: PetscCall(PetscLayoutDestroy(&clayout));
1213: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
1214: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
1215: PetscCall(PetscSFDestroy(&sf));
1216: PetscCall(PetscFree2(cranks_leaf, cranks_root));
1217: for (PetscInt c = 0; c < cEnd - cStart; c++) npoints[ranks[c]]++;
1219: starts[0] = 0;
1220: for (PetscInt c = 0; c < nparts; c++) starts[c + 1] = starts[c] + npoints[c];
1221: for (PetscInt c = 0; c < cEnd - cStart; c++) points[starts[ranks[c]]++] = c;
1222: PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
1223: PetscCall(PetscPartitionerShellSetPartition(part, nparts, npoints, points));
1224: PetscCall(PetscFree(npoints));
1225: PetscCall(PetscFree4(points, ranks, starts, gidxs));
1226: PetscCall(DMRestoreLocalVector(sdm, &map));
1227: PetscCall(DMDestroy(&sdm));
1228: }
1229: PetscCall(PetscSFDestroy(&sfXC));
1230: PetscCall(DMPlexDistribute(dm, 0, &sfBC, odm));
1231: if (*odm) {
1232: PetscCall(DMDestroy(&dm));
1233: dm = *odm;
1234: *odm = NULL;
1235: PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
1236: }
1237: if (sfBC) PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
1238: else {
1239: PetscCall(PetscObjectReference((PetscObject)sfXB));
1240: sfXC = sfXB;
1241: }
1242: PetscCall(PetscSFDestroy(&sfXB));
1243: PetscCall(PetscSFDestroy(&sfBC));
1244: PetscCall(DMSetCoarsenLevel(dm, cc));
1245: PetscCall(DMSetRefineLevel(dm, rr));
1246: PetscCall(DMPlexSetRegularRefinement(dm, isRegular));
1247: PetscCall(DMPlexSetRefinementUniform(dm, isUniform));
1248: PetscCall(DMPlexSectionLoad(dm, viewer, NULL, sfXC, &sfG, NULL));
1249: if (level == numlevels - 1) {
1250: Vec u;
1252: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &u));
1253: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
1254: PetscCall(DMPlexGlobalVectorLoad(dm, viewer, NULL, sfG, u));
1255: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &u));
1256: }
1257: PetscCall(PetscFree(dmname));
1258: PetscCall(PetscSFDestroy(&sfG));
1259: PetscCall(DMSetCoarseDM(dm, cdm));
1260: PetscCall(DMDestroy(&cdm));
1261: PetscCall(PetscViewerHDF5PopGroup(viewer));
1262: cdm = dm;
1263: }
1264: *odm = dm;
1265: PetscCall(PetscViewerPopFormat(viewer));
1266: PetscCall(PetscViewerDestroy(&viewer));
1267: PetscCall(PetscSFDestroy(&sfXC));
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: #else
1270: SETERRQ(comm, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
1271: #endif
1272: }
1274: /*
1275: Setup AuxDM:
1276: - project source function and make it zero-mean
1277: - sample frozen fields for operator splitting
1278: */
1279: static PetscErrorCode ProjectAuxDM(DM dm, PetscReal time, Vec u, AppCtx *ctx)
1280: {
1281: DM dmAux;
1282: Vec la, a;
1283: void *ctxs[NUM_FIELDS + 1];
1284: PetscScalar vals[NUM_FIELDS + 1];
1285: VecScatter sctAux;
1286: PetscDS ds;
1287: PetscErrorCode (*funcs[NUM_FIELDS + 1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1289: PetscFunctionBeginUser;
1290: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &la));
1291: if (!la) {
1292: PetscFE field;
1293: PetscInt fields[NUM_FIELDS];
1294: IS is;
1295: Vec tu, ta;
1296: PetscInt dim;
1298: PetscCall(DMClone(dm, &dmAux));
1299: PetscCall(DMSetNumFields(dmAux, NUM_FIELDS + 1));
1300: for (PetscInt i = 0; i < NUM_FIELDS; i++) {
1301: PetscCall(DMGetField(dm, i, NULL, (PetscObject *)&field));
1302: PetscCall(DMSetField(dmAux, i, NULL, (PetscObject)field));
1303: fields[i] = i;
1304: }
1305: /* PetscFEDuplicate? */
1306: PetscCall(DMGetDimension(dm, &dim));
1307: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, ctx->simplex, "p_", -1, &field));
1308: PetscCall(PetscObjectSetName((PetscObject)field, "source"));
1309: PetscCall(DMSetField(dmAux, NUM_FIELDS, NULL, (PetscObject)field));
1310: PetscCall(PetscFEDestroy(&field));
1311: PetscCall(DMCreateDS(dmAux));
1312: PetscCall(DMCreateSubDM(dmAux, NUM_FIELDS, fields, &is, NULL));
1313: PetscCall(DMGetGlobalVector(dm, &tu));
1314: PetscCall(DMGetGlobalVector(dmAux, &ta));
1315: PetscCall(VecScatterCreate(tu, NULL, ta, is, &sctAux));
1316: PetscCall(DMRestoreGlobalVector(dm, &tu));
1317: PetscCall(DMRestoreGlobalVector(dmAux, &ta));
1318: PetscCall(PetscObjectCompose((PetscObject)dmAux, "scatterAux", (PetscObject)sctAux));
1319: PetscCall(VecScatterDestroy(&sctAux));
1320: PetscCall(ISDestroy(&is));
1321: PetscCall(DMCreateLocalVector(dmAux, &la));
1322: PetscCall(PetscObjectSetName((PetscObject)la, "auxiliary_"));
1323: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, la));
1324: PetscCall(DMDestroy(&dmAux));
1325: PetscCall(VecDestroy(&la));
1326: }
1327: if (time == PETSC_MIN_REAL) PetscFunctionReturn(PETSC_SUCCESS);
1329: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &la));
1330: PetscCall(VecGetDM(la, &dmAux));
1331: PetscCall(DMGetDS(dmAux, &ds));
1332: PetscCall(DMGetGlobalVector(dmAux, &a));
1333: funcs[C_FIELD_ID] = zerof;
1334: ctxs[C_FIELD_ID] = NULL;
1335: funcs[P_FIELD_ID] = zerof;
1336: ctxs[P_FIELD_ID] = NULL;
1337: funcs[NUM_FIELDS] = source;
1338: ctxs[NUM_FIELDS] = ctx->source_ctx;
1339: PetscCall(DMProjectFunction(dmAux, time, funcs, ctxs, INSERT_ALL_VALUES, a));
1340: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1341: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, zero));
1342: PetscCall(PetscDSSetObjective(ds, NUM_FIELDS, average));
1343: PetscCall(DMPlexComputeIntegralFEM(dmAux, a, vals, NULL));
1344: PetscCall(VecShift(a, -vals[NUM_FIELDS] / ctx->domain_volume));
1345: if (u) {
1346: PetscCall(PetscObjectQuery((PetscObject)dmAux, "scatterAux", (PetscObject *)&sctAux));
1347: PetscCheck(sctAux, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing scatterAux");
1348: PetscCall(VecScatterBegin(sctAux, u, a, INSERT_VALUES, SCATTER_FORWARD));
1349: PetscCall(VecScatterEnd(sctAux, u, a, INSERT_VALUES, SCATTER_FORWARD));
1350: }
1351: PetscCall(DMGlobalToLocal(dmAux, a, INSERT_VALUES, la));
1352: PetscCall(VecViewFromOptions(la, NULL, "-aux_view"));
1353: PetscCall(DMRestoreGlobalVector(dmAux, &a));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /* callback for the creation of the potential null space */
1358: static PetscErrorCode CreatePotentialNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
1359: {
1360: Vec vec;
1361: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zerof};
1363: PetscFunctionBeginUser;
1364: funcs[nfield] = constantf;
1365: PetscCall(DMCreateGlobalVector(dm, &vec));
1366: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
1367: PetscCall(VecNormalize(vec, NULL));
1368: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
1369: /* break ref cycles */
1370: PetscCall(VecSetDM(vec, NULL));
1371: PetscCall(VecDestroy(&vec));
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }
1375: static PetscErrorCode DMGetLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
1376: {
1377: PetscBool has;
1379: PetscFunctionBeginUser;
1380: if (local) {
1381: PetscCall(DMHasNamedLocalVector(dm, "lumped_mass", &has));
1382: PetscCall(DMGetNamedLocalVector(dm, "lumped_mass", lumped_mass));
1383: } else {
1384: PetscCall(DMHasNamedGlobalVector(dm, "lumped_mass", &has));
1385: PetscCall(DMGetNamedGlobalVector(dm, "lumped_mass", lumped_mass));
1386: }
1387: if (!has) {
1388: Vec w;
1389: IS is;
1391: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
1392: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
1393: if (local) {
1394: Vec w2, wg;
1396: PetscCall(DMCreateMassMatrixLumped(dm, &w, NULL));
1397: PetscCall(DMGetGlobalVector(dm, &wg));
1398: PetscCall(DMGetLocalVector(dm, &w2));
1399: PetscCall(VecSet(w2, 0.0));
1400: PetscCall(VecSet(wg, 1.0));
1401: PetscCall(VecISSet(wg, is, 0.0));
1402: PetscCall(DMGlobalToLocal(dm, wg, INSERT_VALUES, w2));
1403: PetscCall(VecPointwiseMult(w, w, w2));
1404: PetscCall(DMRestoreGlobalVector(dm, &wg));
1405: PetscCall(DMRestoreLocalVector(dm, &w2));
1406: } else {
1407: PetscCall(DMCreateMassMatrixLumped(dm, NULL, &w));
1408: PetscCall(VecISSet(w, is, 0.0));
1409: }
1410: PetscCall(VecCopy(w, *lumped_mass));
1411: PetscCall(VecDestroy(&w));
1412: }
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: static PetscErrorCode DMRestoreLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
1417: {
1418: PetscFunctionBeginUser;
1419: if (local) PetscCall(DMRestoreNamedLocalVector(dm, "lumped_mass", lumped_mass));
1420: else PetscCall(DMRestoreNamedGlobalVector(dm, "lumped_mass", lumped_mass));
1421: PetscFunctionReturn(PETSC_SUCCESS);
1422: }
1424: /* callbacks for residual and Jacobian */
1425: static PetscErrorCode DMPlexTSComputeIFunctionFEM_Private(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
1426: {
1427: Vec work, local_lumped_mass;
1428: AppCtx *ctx = (AppCtx *)user;
1430: PetscFunctionBeginUser;
1431: if (ctx->fix_c < 0) {
1432: PetscReal vals[NUM_FIELDS];
1433: PetscDS ds;
1435: PetscCall(DMGetDS(dm, &ds));
1436: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, jacobian_fail));
1437: PetscCall(DMPlexSNESComputeObjectiveFEM(dm, locX, vals, NULL));
1438: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1439: if (vals[C_FIELD_ID]) PetscCall(SNESSetFunctionDomainError(ctx->snes));
1440: }
1441: if (ctx->lump) {
1442: PetscCall(DMGetLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
1443: PetscCall(DMGetLocalVector(dm, &work));
1444: PetscCall(VecSet(work, 0.0));
1445: PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, work, locF, user));
1446: PetscCall(VecPointwiseMult(work, locX_t, local_lumped_mass));
1447: PetscCall(VecAXPY(locF, 1.0, work));
1448: PetscCall(DMRestoreLocalVector(dm, &work));
1449: PetscCall(DMRestoreLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
1450: } else {
1451: PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, locX_t, locF, user));
1452: }
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: static PetscErrorCode DMPlexTSComputeIJacobianFEM_Private(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
1457: {
1458: Vec lumped_mass, work;
1459: AppCtx *ctx = (AppCtx *)user;
1461: PetscFunctionBeginUser;
1462: if (ctx->lump) {
1463: PetscCall(DMGetLumpedMass(dm, PETSC_FALSE, &lumped_mass));
1464: PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, 0.0, Jac, JacP, user));
1465: PetscCall(DMGetGlobalVector(dm, &work));
1466: PetscCall(VecAXPBY(work, X_tShift, 0.0, lumped_mass));
1467: PetscCall(MatDiagonalSet(JacP, work, ADD_VALUES));
1468: PetscCall(DMRestoreGlobalVector(dm, &work));
1469: PetscCall(DMRestoreLumpedMass(dm, PETSC_FALSE, &lumped_mass));
1470: } else {
1471: PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, X_tShift, Jac, JacP, user));
1472: }
1473: PetscFunctionReturn(PETSC_SUCCESS);
1474: }
1476: /* customize residuals and Jacobians */
1477: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
1478: {
1479: PetscDS ds;
1480: PetscInt cdim, dim;
1481: PetscScalar constants[NUM_CONSTANTS];
1482: PetscScalar vals[NUM_FIELDS];
1483: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
1484: Vec dummy;
1485: IS is;
1487: PetscFunctionBeginUser;
1488: constants[R_ID] = ctx->r;
1489: constants[EPS_ID] = ctx->eps;
1490: constants[ALPHA_ID] = ctx->alpha;
1491: constants[GAMMA_ID] = ctx->gamma;
1492: constants[D_ID] = ctx->D;
1493: constants[FIXC_ID] = ctx->fix_c;
1494: constants[SPLIT_ID] = ctx->split;
1496: PetscCall(DMGetDimension(dm, &dim));
1497: PetscCall(DMGetCoordinateDim(dm, &cdim));
1498: PetscCheck(dim == 2 || dim == 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Topological dimension must be 2 or 3");
1499: PetscCheck(dim == ctx->dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Topological dimension mismatch: expected %" PetscInt_FMT ", found %" PetscInt_FMT, dim, ctx->dim);
1500: PetscCheck(cdim == ctx->dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Geometrical dimension mismatch: expected %" PetscInt_FMT ", found %" PetscInt_FMT, cdim, ctx->dim);
1501: PetscCall(DMGetDS(dm, &ds));
1502: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
1503: PetscCall(PetscDSSetImplicit(ds, C_FIELD_ID, PETSC_TRUE));
1504: PetscCall(PetscDSSetImplicit(ds, P_FIELD_ID, PETSC_TRUE));
1505: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1506: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1507: if (ctx->dim == 2) {
1508: PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0, C_1));
1509: PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1));
1510: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0, NULL, NULL, ctx->D > 0 ? JC_1_c1c1 : NULL));
1511: if (!ctx->split) { /* we solve a block diagonal system to mimic operator splitting, jacobians will not be correct */
1512: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1, NULL, NULL));
1513: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0, NULL));
1514: }
1515: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1));
1516: } else {
1517: PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0_3d, C_1_3d));
1518: PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1_3d));
1519: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0_3d, NULL, NULL, ctx->D > 0 ? JC_1_c1c1_3d : NULL));
1520: if (!ctx->split) {
1521: PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1_3d, NULL, NULL));
1522: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0_3d, NULL));
1523: }
1524: PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1_3d));
1525: }
1526: /* Attach potential nullspace */
1527: PetscCall(DMSetNullSpaceConstructor(dm, P_FIELD_ID, CreatePotentialNullSpace));
1529: /* Compute domain volume */
1530: PetscCall(DMGetGlobalVector(dm, &dummy));
1531: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, volume));
1532: PetscCall(DMPlexComputeIntegralFEM(dm, dummy, vals, NULL));
1533: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1534: PetscCall(DMRestoreGlobalVector(dm, &dummy));
1535: ctx->domain_volume = PetscRealPart(vals[P_FIELD_ID]);
1537: /* Index sets for potential and conductivities */
1538: PetscCall(DMCreateSubDM(dm, 1, fields, &is, NULL));
1539: PetscCall(PetscObjectCompose((PetscObject)dm, "IS conductivity", (PetscObject)is));
1540: PetscCall(PetscObjectSetName((PetscObject)is, "C"));
1541: PetscCall(ISViewFromOptions(is, NULL, "-is_conductivity_view"));
1542: PetscCall(ISDestroy(&is));
1543: PetscCall(DMCreateSubDM(dm, 1, fields + 1, &is, NULL));
1544: PetscCall(PetscObjectSetName((PetscObject)is, "P"));
1545: PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)is));
1546: PetscCall(ISViewFromOptions(is, NULL, "-is_potential_view"));
1547: PetscCall(ISDestroy(&is));
1549: /* Create auxiliary DM */
1550: PetscCall(ProjectAuxDM(dm, PETSC_MIN_REAL, NULL, ctx));
1552: /* Add callbacks */
1553: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM_Private, ctx));
1554: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM_Private, ctx));
1555: /* DMPlexTSComputeBoundary is not needed because we use natural boundary conditions */
1556: PetscCall(DMTSSetBoundaryLocal(dm, NULL, NULL));
1558: /* handle nullspace in residual (move it to TSComputeIFunction_DMLocal?) */
1559: {
1560: MatNullSpace nullsp;
1562: PetscCall(CreatePotentialNullSpace(dm, P_FIELD_ID, P_FIELD_ID, &nullsp));
1563: PetscCall(PetscObjectCompose((PetscObject)dm, "__dmtsnullspace", (PetscObject)nullsp));
1564: PetscCall(MatNullSpaceDestroy(&nullsp));
1565: }
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: /* setup discrete spaces and residuals */
1570: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
1571: {
1572: DM cdm = dm;
1573: PetscFE feC, feP;
1574: PetscInt dim;
1575: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1577: PetscFunctionBeginUser;
1578: PetscCall(DMGetDimension(dm, &dim));
1580: /* We model Cij with Cij = Cji -> dim*(dim+1)/2 components */
1581: PetscCall(PetscFECreateDefault(comm, dim, (dim * (dim + 1)) / 2, ctx->simplex, "c_", -1, &feC));
1582: PetscCall(PetscObjectSetName((PetscObject)feC, "conductivity"));
1583: PetscCall(PetscFECreateDefault(comm, dim, 1, ctx->simplex, "p_", -1, &feP));
1584: PetscCall(PetscObjectSetName((PetscObject)feP, "potential"));
1585: PetscCall(PetscFECopyQuadrature(feP, feC));
1586: PetscCall(PetscFEViewFromOptions(feP, NULL, "-view_fe"));
1587: PetscCall(PetscFEViewFromOptions(feC, NULL, "-view_fe"));
1589: PetscCall(DMSetNumFields(dm, 2));
1590: PetscCall(DMSetField(dm, C_FIELD_ID, NULL, (PetscObject)feC));
1591: PetscCall(DMSetField(dm, P_FIELD_ID, NULL, (PetscObject)feP));
1592: PetscCall(PetscFEDestroy(&feC));
1593: PetscCall(PetscFEDestroy(&feP));
1594: PetscCall(DMCreateDS(dm));
1595: while (cdm) {
1596: PetscCall(DMCopyDisc(dm, cdm));
1597: PetscCall(SetupProblem(cdm, ctx));
1598: PetscCall(DMGetCoarseDM(cdm, &cdm));
1599: }
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: /* Create mesh by command line options */
1604: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
1605: {
1606: DM plex;
1608: PetscFunctionBeginUser;
1609: if (ctx->load) {
1610: PetscInt refine = 0;
1611: PetscBool isHierarchy;
1612: DM *dms;
1613: char typeName[256];
1614: PetscBool flg;
1616: PetscCall(LoadFromFile(comm, ctx->load_filename, dm));
1617: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
1618: PetscCall(PetscOptionsFList("-dm_mat_type", "Matrix type used for created matrices", "DMSetMatType", MatList, MATAIJ, typeName, sizeof(typeName), &flg));
1619: if (flg) PetscCall(DMSetMatType(*dm, typeName));
1620: PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0));
1621: PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0));
1622: PetscOptionsEnd();
1623: if (refine) {
1624: PetscCall(SetupDiscretization(*dm, ctx));
1625: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
1626: }
1627: PetscCall(PetscCalloc1(refine, &dms));
1628: if (isHierarchy) PetscCall(DMRefineHierarchy(*dm, refine, dms));
1629: for (PetscInt r = 0; r < refine; r++) {
1630: Mat M;
1631: DM dmr = dms[r];
1632: Vec u, ur;
1634: if (!isHierarchy) {
1635: PetscCall(DMRefine(*dm, PetscObjectComm((PetscObject)*dm), &dmr));
1636: PetscCall(DMSetCoarseDM(dmr, *dm));
1637: }
1638: PetscCall(DMCreateInterpolation(*dm, dmr, &M, NULL));
1639: PetscCall(DMGetNamedGlobalVector(*dm, "solution_", &u));
1640: PetscCall(DMGetNamedGlobalVector(dmr, "solution_", &ur));
1641: PetscCall(MatInterpolate(M, u, ur));
1642: PetscCall(DMRestoreNamedGlobalVector(*dm, "solution_", &u));
1643: PetscCall(DMRestoreNamedGlobalVector(dmr, "solution_", &ur));
1644: PetscCall(MatDestroy(&M));
1645: if (!isHierarchy) PetscCall(DMSetCoarseDM(dmr, NULL));
1646: PetscCall(DMDestroy(dm));
1647: *dm = dmr;
1648: }
1649: if (refine && !isHierarchy) PetscCall(DMSetRefineLevel(*dm, 0));
1650: PetscCall(PetscFree(dms));
1651: } else {
1652: PetscCall(DMCreate(comm, dm));
1653: PetscCall(DMSetType(*dm, DMPLEX));
1654: PetscCall(DMSetFromOptions(*dm));
1655: PetscCall(DMLocalizeCoordinates(*dm));
1656: {
1657: char convType[256];
1658: PetscBool flg;
1659: PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
1660: PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg));
1661: PetscOptionsEnd();
1662: if (flg) {
1663: DM dmConv;
1664: PetscCall(DMConvert(*dm, convType, &dmConv));
1665: if (dmConv) {
1666: PetscCall(DMDestroy(dm));
1667: *dm = dmConv;
1668: PetscCall(DMSetFromOptions(*dm));
1669: PetscCall(DMSetUp(*dm));
1670: }
1671: }
1672: }
1673: }
1674: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "ref_"));
1675: PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
1676: PetscCall(DMSetFromOptions(*dm));
1677: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
1679: PetscCall(DMConvert(*dm, DMPLEX, &plex));
1680: PetscCall(DMPlexIsSimplex(plex, &ctx->simplex));
1681: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->simplex, 1, MPIU_BOOL, MPI_LOR, comm));
1682: PetscCall(DMDestroy(&plex));
1684: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /* Make potential field zero mean */
1689: static PetscErrorCode ZeroMeanPotential(DM dm, Vec u, PetscScalar domain_volume)
1690: {
1691: PetscScalar vals[NUM_FIELDS];
1692: PetscDS ds;
1693: IS is;
1695: PetscFunctionBeginUser;
1696: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
1697: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
1698: PetscCall(DMGetDS(dm, &ds));
1699: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
1700: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1701: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
1702: PetscCall(VecISShift(u, is, -vals[P_FIELD_ID] / domain_volume));
1703: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /* Compute initial conditions and exclude potential from local truncation error
1708: Since we are solving a DAE, once the initial conditions for the differential
1709: variables are set, we need to compute the corresponding value for the
1710: algebraic variables. We do so by creating a subDM for the potential only
1711: and solve a static problem with SNES */
1712: static PetscErrorCode SetInitialConditionsAndTolerances(TS ts, PetscInt nv, Vec vecs[], PetscBool valid)
1713: {
1714: DM dm;
1715: Vec u, p, subaux, vatol, vrtol;
1716: PetscReal t, atol, rtol;
1717: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
1718: IS isp;
1719: DM dmp;
1720: VecScatter sctp = NULL;
1721: PetscDS ds;
1722: SNES snes;
1723: KSP ksp;
1724: PC pc;
1725: AppCtx *ctx;
1727: PetscFunctionBeginUser;
1728: PetscCall(TSGetDM(ts, &dm));
1729: PetscCall(TSGetApplicationContext(ts, &ctx));
1730: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&isp));
1731: PetscCheck(isp, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
1732: if (valid) {
1733: if (ctx->exclude_potential_lte) {
1734: PetscCall(DMCreateGlobalVector(dm, &vatol));
1735: PetscCall(DMCreateGlobalVector(dm, &vrtol));
1736: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
1737: PetscCall(VecSet(vatol, atol));
1738: PetscCall(VecISSet(vatol, isp, -1));
1739: PetscCall(VecSet(vrtol, rtol));
1740: PetscCall(VecISSet(vrtol, isp, -1));
1741: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
1742: PetscCall(VecDestroy(&vatol));
1743: PetscCall(VecDestroy(&vrtol));
1744: }
1745: for (PetscInt i = 0; i < nv; i++) PetscCall(ZeroMeanPotential(dm, vecs[i], ctx->domain_volume));
1746: PetscFunctionReturn(PETSC_SUCCESS);
1747: }
1748: PetscCall(DMCreateSubDM(dm, 1, fields + 1, NULL, &dmp));
1749: PetscCall(DMSetMatType(dmp, MATAIJ));
1750: PetscCall(DMGetDS(dmp, &ds));
1751: if (ctx->dim == 2) {
1752: PetscCall(PetscDSSetResidual(ds, 0, P_0_aux, P_1_aux));
1753: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux));
1754: } else {
1755: PetscCall(PetscDSSetResidual(ds, 0, P_0_aux, P_1_aux_3d));
1756: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux_3d));
1757: }
1758: PetscCall(DMPlexSetSNESLocalFEM(dmp, PETSC_FALSE, NULL));
1760: PetscCall(DMCreateGlobalVector(dmp, &p));
1762: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dmp), &snes));
1763: PetscCall(SNESSetDM(snes, dmp));
1764: PetscCall(SNESSetOptionsPrefix(snes, "initial_"));
1765: PetscCall(SNESSetErrorIfNotConverged(snes, PETSC_TRUE));
1766: PetscCall(SNESGetKSP(snes, &ksp));
1767: PetscCall(KSPSetType(ksp, KSPFGMRES));
1768: PetscCall(KSPGetPC(ksp, &pc));
1769: PetscCall(PCSetType(pc, PCGAMG));
1770: PetscCall(SNESSetFromOptions(snes));
1771: PetscCall(SNESSetUp(snes));
1773: /* Loop over input vectors and compute corresponding potential */
1774: for (PetscInt i = 0; i < nv; i++) {
1775: u = vecs[i];
1776: if (!valid) { /* Assumes entries in u are not valid */
1777: PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1778: void *ctxs[NUM_FIELDS] = {NULL};
1779: PetscRandom r = NULL;
1781: PetscCall(TSGetTime(ts, &t));
1782: switch (ctx->ic_num) {
1783: case 0:
1784: funcs[C_FIELD_ID] = initial_conditions_C_0;
1785: break;
1786: case 1:
1787: funcs[C_FIELD_ID] = initial_conditions_C_1;
1788: break;
1789: case 2:
1790: funcs[C_FIELD_ID] = initial_conditions_C_2;
1791: break;
1792: case 3:
1793: funcs[C_FIELD_ID] = initial_conditions_C_random;
1794: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ts), &r));
1795: PetscCall(PetscRandomSetOptionsPrefix(r, "ic_"));
1796: PetscCall(PetscRandomSetFromOptions(r));
1797: ctxs[C_FIELD_ID] = r;
1798: break;
1799: default:
1800: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unknown IC");
1801: }
1802: funcs[P_FIELD_ID] = zerof;
1803: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_ALL_VALUES, u));
1804: PetscCall(ProjectAuxDM(dm, t, u, ctx));
1805: PetscCall(PetscRandomDestroy(&r));
1806: }
1808: /* pass conductivity information via auxiliary data */
1809: PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &subaux));
1810: PetscCall(DMSetAuxiliaryVec(dmp, NULL, 0, 0, subaux));
1812: /* solve the subproblem */
1813: if (!sctp) PetscCall(VecScatterCreate(u, isp, p, NULL, &sctp));
1814: PetscCall(VecScatterBegin(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1815: PetscCall(VecScatterEnd(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
1816: PetscCall(SNESSolve(snes, NULL, p));
1818: /* scatter from potential only to full space */
1819: PetscCall(VecScatterBegin(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1820: PetscCall(VecScatterEnd(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
1821: PetscCall(ZeroMeanPotential(dm, u, ctx->domain_volume));
1822: }
1823: PetscCall(VecDestroy(&p));
1824: PetscCall(DMDestroy(&dmp));
1825: PetscCall(SNESDestroy(&snes));
1826: PetscCall(VecScatterDestroy(&sctp));
1828: /* exclude potential from computation of the LTE */
1829: if (ctx->exclude_potential_lte) {
1830: PetscCall(DMCreateGlobalVector(dm, &vatol));
1831: PetscCall(DMCreateGlobalVector(dm, &vrtol));
1832: PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
1833: PetscCall(VecSet(vatol, atol));
1834: PetscCall(VecISSet(vatol, isp, -1));
1835: PetscCall(VecSet(vrtol, rtol));
1836: PetscCall(VecISSet(vrtol, isp, -1));
1837: PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
1838: PetscCall(VecDestroy(&vatol));
1839: PetscCall(VecDestroy(&vrtol));
1840: }
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: /* mesh adaption context */
1845: typedef struct {
1846: VecTagger refineTag;
1847: DMLabel adaptLabel;
1848: PetscInt cnt;
1849: } AdaptCtx;
1851: static PetscErrorCode ResizeSetUp(TS ts, PetscInt nstep, PetscReal time, Vec u, PetscBool *resize, void *vctx)
1852: {
1853: AdaptCtx *ctx = (AdaptCtx *)vctx;
1854: Vec ellVecCells, ellVecCellsF;
1855: DM dm, plex;
1856: PetscDS ds;
1857: PetscReal norm;
1858: PetscInt cStart, cEnd;
1860: PetscFunctionBeginUser;
1861: PetscCall(TSGetDM(ts, &dm));
1862: PetscCall(DMConvert(dm, DMPLEX, &plex));
1863: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
1864: PetscCall(DMDestroy(&plex));
1865: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), NUM_FIELDS * (cEnd - cStart), PETSC_DECIDE, &ellVecCellsF));
1866: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), cEnd - cStart, PETSC_DECIDE, &ellVecCells));
1867: PetscCall(VecSetBlockSize(ellVecCellsF, NUM_FIELDS));
1868: PetscCall(DMGetDS(dm, &ds));
1869: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
1870: PetscCall(DMPlexComputeCellwiseIntegralFEM(dm, u, ellVecCellsF, NULL));
1871: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
1872: PetscCall(VecStrideGather(ellVecCellsF, C_FIELD_ID, ellVecCells, INSERT_VALUES));
1873: PetscCall(VecDestroy(&ellVecCellsF));
1874: PetscCall(VecNorm(ellVecCells, NORM_1, &norm));
1875: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "STEP %d norm %g\n", (int)nstep, (double)norm));
1876: if (norm && !ctx->cnt) {
1877: IS refineIS;
1879: *resize = PETSC_TRUE;
1880: if (!ctx->refineTag) {
1881: VecTaggerBox refineBox;
1882: refineBox.min = PETSC_MACHINE_EPSILON;
1883: refineBox.max = PETSC_MAX_REAL;
1885: PetscCall(VecTaggerCreate(PETSC_COMM_SELF, &ctx->refineTag));
1886: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->refineTag, "refine_"));
1887: PetscCall(VecTaggerSetType(ctx->refineTag, VECTAGGERABSOLUTE));
1888: PetscCall(VecTaggerAbsoluteSetBox(ctx->refineTag, &refineBox));
1889: PetscCall(VecTaggerSetFromOptions(ctx->refineTag));
1890: PetscCall(VecTaggerSetUp(ctx->refineTag));
1891: PetscCall(PetscObjectViewFromOptions((PetscObject)ctx->refineTag, NULL, "-tag_view"));
1892: }
1893: PetscCall(DMLabelDestroy(&ctx->adaptLabel));
1894: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)ts), "adapt", &ctx->adaptLabel));
1895: PetscCall(VecTaggerComputeIS(ctx->refineTag, ellVecCells, &refineIS, NULL));
1896: PetscCall(DMLabelSetStratumIS(ctx->adaptLabel, DM_ADAPT_REFINE, refineIS));
1897: PetscCall(ISDestroy(&refineIS));
1898: #if 0
1899: 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[]);
1900: Vec ellVec;
1902: funcs[P_FIELD_ID] = ellipticity_fail;
1903: funcs[C_FIELD_ID] = NULL;
1905: PetscCall(DMGetGlobalVector(dm, &ellVec));
1906: PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
1907: PetscCall(VecViewFromOptions(ellVec,NULL,"-view_amr_ell"));
1908: PetscCall(DMRestoreGlobalVector(dm, &ellVec));
1909: #endif
1910: ctx->cnt++;
1911: } else {
1912: ctx->cnt = 0;
1913: }
1914: PetscCall(VecDestroy(&ellVecCells));
1915: PetscFunctionReturn(PETSC_SUCCESS);
1916: }
1918: static PetscErrorCode ResizeTransfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *vctx)
1919: {
1920: AdaptCtx *actx = (AdaptCtx *)vctx;
1921: AppCtx *ctx;
1922: DM dm, adm;
1923: PetscReal time;
1924: PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
1925: IS is;
1927: PetscFunctionBeginUser;
1928: PetscCall(TSGetDM(ts, &dm));
1929: PetscCheck(actx->adaptLabel, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptLabel");
1930: PetscCall(DMAdaptLabel(dm, actx->adaptLabel, &adm));
1931: PetscCheck(adm, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adapted DM");
1932: PetscCall(TSGetTime(ts, &time));
1933: PetscCall(DMLabelDestroy(&actx->adaptLabel));
1934: for (PetscInt i = 0; i < nv; i++) {
1935: PetscCall(DMCreateGlobalVector(adm, &vecsout[i]));
1936: PetscCall(DMForestTransferVec(dm, vecsin[i], adm, vecsout[i], PETSC_TRUE, time));
1937: }
1938: PetscCall(DMForestSetAdaptivityForest(adm, NULL));
1939: PetscCall(DMSetCoarseDM(adm, NULL));
1940: PetscCall(DMSetLocalSection(adm, NULL));
1941: PetscCall(TSSetDM(ts, adm));
1942: PetscCall(TSGetTime(ts, &time));
1943: PetscCall(TSGetApplicationContext(ts, &ctx));
1944: PetscCall(DMSetNullSpaceConstructor(adm, P_FIELD_ID, CreatePotentialNullSpace));
1945: PetscCall(DMCreateSubDM(adm, 1, fields, &is, NULL));
1946: PetscCall(PetscObjectCompose((PetscObject)adm, "IS conductivity", (PetscObject)is));
1947: PetscCall(ISDestroy(&is));
1948: PetscCall(DMCreateSubDM(adm, 1, fields + 1, &is, NULL));
1949: PetscCall(PetscObjectCompose((PetscObject)adm, "IS potential", (PetscObject)is));
1950: PetscCall(ISDestroy(&is));
1951: PetscCall(ProjectAuxDM(adm, time, NULL, ctx));
1952: {
1953: MatNullSpace nullsp;
1955: PetscCall(CreatePotentialNullSpace(adm, P_FIELD_ID, P_FIELD_ID, &nullsp));
1956: PetscCall(PetscObjectCompose((PetscObject)adm, "__dmtsnullspace", (PetscObject)nullsp));
1957: PetscCall(MatNullSpaceDestroy(&nullsp));
1958: }
1959: PetscCall(SetInitialConditionsAndTolerances(ts, nv, vecsout, PETSC_TRUE));
1960: PetscCall(DMDestroy(&ctx->view_dm));
1961: for (PetscInt i = 0; i < NUM_FIELDS; i++) {
1962: PetscCall(VecScatterDestroy(&ctx->subsct[i]));
1963: PetscCall(VecDestroy(&ctx->subvec[i]));
1964: }
1965: PetscFunctionReturn(PETSC_SUCCESS);
1966: }
1968: static PetscErrorCode ComputeDiagnostic(Vec u, AppCtx *ctx, Vec *out)
1969: {
1970: DM dm;
1971: PetscInt dim;
1972: PetscFE feFluxC, feNormC, feEigsC;
1974: 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, eigsc, flux};
1976: PetscFunctionBeginUser;
1977: if (!ctx->view_dm) {
1978: PetscFE feP;
1979: PetscInt nf = PetscMax(PetscMin(ctx->diagnostic_num, 3), 1);
1981: PetscCall(VecGetDM(u, &dm));
1982: PetscCall(DMGetDimension(dm, &dim));
1983: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, ctx->simplex, "normc_", -1, &feNormC));
1984: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, ctx->simplex, "eigsc_", -1, &feEigsC));
1985: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, ctx->simplex, "flux_", -1, &feFluxC));
1986: PetscCall(DMGetField(dm, P_FIELD_ID, NULL, (PetscObject *)&feP));
1987: PetscCall(PetscFECopyQuadrature(feP, feNormC));
1988: PetscCall(PetscFECopyQuadrature(feP, feEigsC));
1989: PetscCall(PetscFECopyQuadrature(feP, feFluxC));
1990: PetscCall(PetscObjectSetName((PetscObject)feNormC, "normC"));
1991: PetscCall(PetscObjectSetName((PetscObject)feEigsC, "eigsC"));
1992: PetscCall(PetscObjectSetName((PetscObject)feFluxC, "flux"));
1994: PetscCall(DMClone(dm, &ctx->view_dm));
1995: PetscCall(DMSetNumFields(ctx->view_dm, nf));
1996: PetscCall(DMSetField(ctx->view_dm, 0, NULL, (PetscObject)feNormC));
1997: if (nf > 1) PetscCall(DMSetField(ctx->view_dm, 1, NULL, (PetscObject)feEigsC));
1998: if (nf > 2) PetscCall(DMSetField(ctx->view_dm, 2, NULL, (PetscObject)feFluxC));
1999: PetscCall(DMCreateDS(ctx->view_dm));
2000: PetscCall(PetscFEDestroy(&feFluxC));
2001: PetscCall(PetscFEDestroy(&feNormC));
2002: PetscCall(PetscFEDestroy(&feEigsC));
2003: }
2004: PetscCall(DMCreateGlobalVector(ctx->view_dm, out));
2005: PetscCall(DMProjectField(ctx->view_dm, 0.0, u, funcs, INSERT_VALUES, *out));
2006: PetscFunctionReturn(PETSC_SUCCESS);
2007: }
2009: static PetscErrorCode MakeScatterAndVec(Vec X, IS is, Vec *Y, VecScatter *sct)
2010: {
2011: PetscInt n;
2013: PetscFunctionBeginUser;
2014: PetscCall(ISGetLocalSize(is, &n));
2015: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)X), n, PETSC_DECIDE, Y));
2016: PetscCall(VecScatterCreate(X, is, *Y, NULL, sct));
2017: PetscFunctionReturn(PETSC_SUCCESS);
2018: }
2020: static PetscErrorCode FunctionDomainError(TS ts, PetscReal time, Vec X, PetscBool *accept)
2021: {
2022: AppCtx *ctx;
2023: PetscScalar vals[NUM_FIELDS];
2024: DM dm;
2025: PetscDS ds;
2027: PetscFunctionBeginUser;
2028: *accept = PETSC_TRUE;
2029: PetscCall(TSGetApplicationContext(ts, &ctx));
2030: if (ctx->function_domain_error_tol < 0) PetscFunctionReturn(PETSC_SUCCESS);
2031: PetscCall(TSGetDM(ts, &dm));
2032: PetscCall(DMGetDS(dm, &ds));
2033: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
2034: PetscCall(DMPlexComputeIntegralFEM(dm, X, vals, NULL));
2035: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
2036: if (PetscAbsScalar(vals[C_FIELD_ID]) > ctx->function_domain_error_tol) *accept = PETSC_FALSE;
2037: if (!*accept) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Domain error value %g > %g\n", (double)PetscAbsScalar(vals[C_FIELD_ID]), (double)ctx->function_domain_error_tol));
2038: PetscFunctionReturn(PETSC_SUCCESS);
2039: }
2041: /* Monitor relevant functionals */
2042: static PetscErrorCode Monitor(TS ts, PetscInt stepnum, PetscReal time, Vec u, void *vctx)
2043: {
2044: PetscScalar vals[2 * NUM_FIELDS];
2045: DM dm;
2046: PetscDS ds;
2047: AppCtx *ctx;
2048: PetscBool need_hdf5, need_vtk;
2050: PetscFunctionBeginUser;
2051: PetscCall(TSGetDM(ts, &dm));
2052: PetscCall(TSGetApplicationContext(ts, &ctx));
2053: PetscCall(DMGetDS(dm, &ds));
2055: /* monitor energy and potential average */
2056: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
2057: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
2058: PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
2060: /* monitor ellipticity_fail */
2061: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
2062: PetscCall(DMPlexComputeIntegralFEM(dm, u, vals + NUM_FIELDS, NULL));
2063: PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
2064: vals[C_FIELD_ID] /= ctx->domain_volume;
2066: 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])));
2068: /* monitor diagnostic */
2069: need_hdf5 = (PetscBool)(ctx->view_hdf5_ctx && ((ctx->view_hdf5_ctx->view_interval > 0 && !(stepnum % ctx->view_hdf5_ctx->view_interval)) || (ctx->view_hdf5_ctx->view_interval && ts->reason)));
2070: need_vtk = (PetscBool)(ctx->view_vtk_ctx && ((ctx->view_vtk_ctx->interval > 0 && !(stepnum % ctx->view_vtk_ctx->interval)) || (ctx->view_vtk_ctx->interval && ts->reason)));
2071: if (ctx->view_times_k < ctx->view_times_n && time >= ctx->view_times[ctx->view_times_k] && time < ctx->view_times[ctx->view_times_k + 1]) {
2072: if (ctx->view_hdf5_ctx) need_hdf5 = PETSC_TRUE;
2073: if (ctx->view_vtk_ctx) need_vtk = PETSC_TRUE;
2074: ctx->view_times_k++;
2075: }
2076: if (need_vtk || need_hdf5) {
2077: Vec diagnostic;
2078: PetscBool dump_dm = (PetscBool)(!!ctx->view_dm);
2080: PetscCall(ComputeDiagnostic(u, ctx, &diagnostic));
2081: if (need_vtk) {
2082: PetscCall(PetscObjectSetName((PetscObject)diagnostic, ""));
2083: PetscCall(TSMonitorSolutionVTK(ts, stepnum, time, diagnostic, ctx->view_vtk_ctx));
2084: }
2085: if (need_hdf5) {
2086: DM vdm;
2087: PetscInt seqnum;
2089: PetscCall(VecGetDM(diagnostic, &vdm));
2090: if (!dump_dm) {
2091: PetscCall(PetscViewerPushFormat(ctx->view_hdf5_ctx->viewer, ctx->view_hdf5_ctx->format));
2092: PetscCall(DMView(vdm, ctx->view_hdf5_ctx->viewer));
2093: PetscCall(PetscViewerPopFormat(ctx->view_hdf5_ctx->viewer));
2094: }
2095: PetscCall(DMGetOutputSequenceNumber(vdm, &seqnum, NULL));
2096: PetscCall(DMSetOutputSequenceNumber(vdm, seqnum + 1, time));
2097: PetscCall(PetscObjectSetName((PetscObject)diagnostic, "diagnostic"));
2098: PetscCall(PetscViewerPushFormat(ctx->view_hdf5_ctx->viewer, ctx->view_hdf5_ctx->format));
2099: PetscCall(VecView(diagnostic, ctx->view_hdf5_ctx->viewer));
2100: if (ctx->diagnostic_num > 3 || ctx->diagnostic_num < 0) {
2101: PetscCall(DMSetOutputSequenceNumber(dm, seqnum + 1, time));
2102: PetscCall(VecView(u, ctx->view_hdf5_ctx->viewer));
2103: }
2104: PetscCall(PetscViewerPopFormat(ctx->view_hdf5_ctx->viewer));
2105: }
2106: PetscCall(VecDestroy(&diagnostic));
2107: }
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: /* Save restart information */
2112: static PetscErrorCode MonitorSave(TS ts, PetscInt steps, PetscReal time, Vec u, void *vctx)
2113: {
2114: DM dm;
2115: AppCtx *ctx = (AppCtx *)vctx;
2116: PetscInt save_every = ctx->save_every;
2117: TSConvergedReason reason;
2119: PetscFunctionBeginUser;
2120: if (!ctx->save) PetscFunctionReturn(PETSC_SUCCESS);
2121: PetscCall(TSGetDM(ts, &dm));
2122: PetscCall(TSGetConvergedReason(ts, &reason));
2123: if ((save_every > 0 && steps % save_every == 0) || (save_every == -1 && reason) || save_every < -1) PetscCall(SaveToFile(dm, u, ctx->save_filename));
2124: PetscFunctionReturn(PETSC_SUCCESS);
2125: }
2127: /* Resample source if time dependent */
2128: static PetscErrorCode PreStage(TS ts, PetscReal stagetime)
2129: {
2130: AppCtx *ctx;
2131: PetscBool resample, ismatis;
2132: Mat A, P;
2134: PetscFunctionBeginUser;
2135: PetscCall(TSGetApplicationContext(ts, &ctx));
2136: /* in case we need to call SNESSetFunctionDomainError */
2137: PetscCall(TSGetSNES(ts, &ctx->snes));
2139: resample = ctx->split ? PETSC_TRUE : PETSC_FALSE;
2140: for (PetscInt i = 0; i < ctx->source_ctx->n; i++) {
2141: if (ctx->source_ctx->k[i] != 0.0) {
2142: resample = PETSC_TRUE;
2143: break;
2144: }
2145: }
2146: if (resample) {
2147: DM dm;
2148: Vec u = NULL;
2150: PetscCall(TSGetDM(ts, &dm));
2151: /* In case of a multistage method, we always use the frozen values at the previous time-step */
2152: if (ctx->split) PetscCall(TSGetSolution(ts, &u));
2153: PetscCall(ProjectAuxDM(dm, stagetime, u, ctx));
2154: }
2156: /* element matrices are sparse, ignore zero entries */
2157: PetscCall(TSGetIJacobian(ts, &A, &P, NULL, NULL));
2158: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
2159: if (!ismatis) PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2160: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis));
2161: if (!ismatis) PetscCall(MatSetOption(P, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2163: /* Set symmetric flag */
2164: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
2165: PetscCall(MatSetOption(P, MAT_SYMMETRIC, PETSC_TRUE));
2166: PetscCall(MatSetOption(A, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
2167: PetscCall(MatSetOption(P, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
2168: PetscFunctionReturn(PETSC_SUCCESS);
2169: }
2171: /* Make potential zero mean after SNES solve */
2172: static PetscErrorCode PostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
2173: {
2174: DM dm;
2175: Vec u = Y[stageindex];
2176: SNES snes;
2177: PetscInt nits, lits, stepnum;
2178: AppCtx *ctx;
2180: PetscFunctionBeginUser;
2181: PetscCall(TSGetDM(ts, &dm));
2182: PetscCall(TSGetApplicationContext(ts, &ctx));
2184: PetscCall(ZeroMeanPotential(dm, u, ctx->domain_volume));
2186: if (ctx->test_restart) PetscFunctionReturn(PETSC_SUCCESS);
2188: /* monitor linear and nonlinear iterations */
2189: PetscCall(TSGetStepNumber(ts, &stepnum));
2190: PetscCall(TSGetSNES(ts, &snes));
2191: PetscCall(SNESGetIterationNumber(snes, &nits));
2192: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
2194: /* if function evals in TSDIRK are zero in the first stage, it is FSAL */
2195: if (stageindex == 0) {
2196: PetscBool dirk;
2197: PetscInt nf;
2199: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2200: PetscCall(SNESGetNumberFunctionEvals(snes, &nf));
2201: if (dirk && nf == 0) nits = 0;
2202: }
2203: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), " step %" PetscInt_FMT " stage %" PetscInt_FMT " nonlinear its %" PetscInt_FMT ", linear its %" PetscInt_FMT "\n", stepnum, stageindex, nits, lits));
2204: PetscFunctionReturn(PETSC_SUCCESS);
2205: }
2207: PetscErrorCode MonitorNorms(SNES snes, PetscInt its, PetscReal f, void *vctx)
2208: {
2209: AppCtx *ctx = (AppCtx *)vctx;
2210: Vec F;
2211: DM dm;
2212: PetscReal subnorm[NUM_FIELDS];
2214: PetscFunctionBeginUser;
2215: PetscCall(SNESGetDM(snes, &dm));
2216: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
2217: if (!ctx->subsct[C_FIELD_ID]) {
2218: IS is;
2220: PetscCall(PetscObjectQuery((PetscObject)dm, "IS conductivity", (PetscObject *)&is));
2221: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing conductivity IS");
2222: PetscCall(MakeScatterAndVec(F, is, &ctx->subvec[C_FIELD_ID], &ctx->subsct[C_FIELD_ID]));
2223: }
2224: if (!ctx->subsct[P_FIELD_ID]) {
2225: IS is;
2227: PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
2228: PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
2229: PetscCall(MakeScatterAndVec(F, is, &ctx->subvec[P_FIELD_ID], &ctx->subsct[P_FIELD_ID]));
2230: }
2231: PetscCall(VecScatterBegin(ctx->subsct[C_FIELD_ID], F, ctx->subvec[C_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD));
2232: PetscCall(VecScatterEnd(ctx->subsct[C_FIELD_ID], F, ctx->subvec[C_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD));
2233: PetscCall(VecScatterBegin(ctx->subsct[P_FIELD_ID], F, ctx->subvec[P_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD));
2234: PetscCall(VecScatterEnd(ctx->subsct[P_FIELD_ID], F, ctx->subvec[P_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD));
2235: PetscCall(VecNorm(ctx->subvec[C_FIELD_ID], NORM_2, &subnorm[C_FIELD_ID]));
2236: PetscCall(VecNorm(ctx->subvec[P_FIELD_ID], NORM_2, &subnorm[P_FIELD_ID]));
2237: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " %3" PetscInt_FMT " SNES Function norms %14.12e, %14.12e\n", its, (double)subnorm[C_FIELD_ID], (double)subnorm[P_FIELD_ID]));
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: }
2241: static PetscErrorCode Run(MPI_Comm comm, AppCtx *ctx)
2242: {
2243: DM dm;
2244: TS ts;
2245: Vec u;
2246: AdaptCtx *actx;
2247: PetscBool flg;
2249: PetscFunctionBeginUser;
2250: if (ctx->test_restart) {
2251: PetscViewer subviewer;
2252: PetscMPIInt rank;
2254: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2255: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
2256: if (ctx->load) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d loading from %s\n", rank, ctx->load_filename));
2257: if (ctx->save) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d saving to %s\n", rank, ctx->save_filename));
2258: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
2259: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
2260: } else {
2261: PetscCall(PetscPrintf(comm, "----------------------------\n"));
2262: PetscCall(PetscPrintf(comm, "Simulation parameters:\n"));
2263: PetscCall(PetscPrintf(comm, " dim : %" PetscInt_FMT "\n", ctx->dim));
2264: PetscCall(PetscPrintf(comm, " r : %g\n", (double)ctx->r));
2265: PetscCall(PetscPrintf(comm, " eps : %g\n", (double)ctx->eps));
2266: PetscCall(PetscPrintf(comm, " alpha: %g\n", (double)ctx->alpha));
2267: PetscCall(PetscPrintf(comm, " gamma: %g\n", (double)ctx->gamma));
2268: PetscCall(PetscPrintf(comm, " D : %g\n", (double)ctx->D));
2269: if (ctx->load) PetscCall(PetscPrintf(comm, " load : %s\n", ctx->load_filename));
2270: else PetscCall(PetscPrintf(comm, " IC : %" PetscInt_FMT "\n", ctx->ic_num));
2271: PetscCall(PetscPrintf(comm, " snum : %" PetscInt_FMT "\n", ctx->source_ctx->n));
2272: for (PetscInt i = 0; i < ctx->source_ctx->n; i++) {
2273: const PetscReal *x0 = ctx->source_ctx->x0 + ctx->dim * i;
2274: const PetscReal w = ctx->source_ctx->w[i];
2275: const PetscReal k = ctx->source_ctx->k[i];
2276: const PetscReal p = ctx->source_ctx->p[i];
2277: const PetscReal r = ctx->source_ctx->r[i];
2279: if (ctx->dim == 2) {
2280: PetscCall(PetscPrintf(comm, " x0 : (%g,%g)\n", (double)x0[0], (double)x0[1]));
2281: } else {
2282: PetscCall(PetscPrintf(comm, " x0 : (%g,%g,%g)\n", (double)x0[0], (double)x0[1], (double)x0[2]));
2283: }
2284: PetscCall(PetscPrintf(comm, " scals: (%g,%g,%g,%g)\n", (double)w, (double)k, (double)p, (double)r));
2285: }
2286: PetscCall(PetscPrintf(comm, "----------------------------\n"));
2287: }
2289: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SetupStage));
2290: PetscCall(CreateMesh(comm, &dm, ctx));
2291: PetscCall(SetupDiscretization(dm, ctx));
2293: PetscCall(TSCreate(comm, &ts));
2294: PetscCall(TSSetApplicationContext(ts, ctx));
2296: PetscCall(TSSetDM(ts, dm));
2297: if (ctx->test_restart) {
2298: PetscViewer subviewer;
2299: PetscMPIInt rank;
2300: PetscInt level;
2302: PetscCall(DMGetRefineLevel(dm, &level));
2303: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2304: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
2305: PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d DM refinement level %" PetscInt_FMT "\n", rank, level));
2306: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
2307: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
2308: }
2310: if (ctx->test_restart) PetscCall(TSSetMaxSteps(ts, 1));
2311: PetscCall(TSSetMaxTime(ts, 10.0));
2312: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
2313: if (!ctx->test_restart) PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
2314: PetscCall(TSMonitorSet(ts, MonitorSave, ctx, NULL));
2315: PetscCall(PetscNew(&actx));
2316: if (ctx->amr) PetscCall(TSSetResize(ts, PETSC_TRUE, ResizeSetUp, ResizeTransfer, actx));
2317: PetscCall(TSSetPreStage(ts, PreStage));
2318: PetscCall(TSSetPostStage(ts, PostStage));
2319: PetscCall(TSSetMaxSNESFailures(ts, -1));
2320: PetscCall(TSSetFunctionDomainError(ts, FunctionDomainError));
2321: PetscCall(TSSetFromOptions(ts));
2322: if (ctx->monitor_norms) {
2323: SNES snes;
2325: PetscCall(TSGetSNES(ts, &snes));
2326: PetscCall(SNESMonitorSet(snes, MonitorNorms, ctx, NULL));
2327: }
2329: PetscCall(DMCreateGlobalVector(dm, &u));
2330: PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
2331: PetscCall(DMHasNamedGlobalVector(dm, "solution_", &flg));
2332: if (flg) { /* load from restart file */
2333: Vec ru;
2335: PetscCall(DMGetNamedGlobalVector(dm, "solution_", &ru));
2336: PetscCall(VecCopy(ru, u));
2337: PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &ru));
2338: }
2339: PetscCall(SetInitialConditionsAndTolerances(ts, 1, &u, flg));
2340: PetscCall(TSSetSolution(ts, u));
2341: PetscCall(VecDestroy(&u));
2342: PetscCall(DMDestroy(&dm));
2343: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
2345: if (!ctx->test_restart) PetscCall(PetscLogStagePush(SolveStage));
2346: PetscCall(TSSolve(ts, NULL));
2347: if (!ctx->test_restart) PetscCall(PetscLogStagePop());
2348: if (ctx->view_vtk_ctx) PetscCall(TSMonitorSolutionVTKDestroy(&ctx->view_vtk_ctx));
2349: if (ctx->view_hdf5_ctx) PetscCall(PetscViewerAndFormatDestroy(&ctx->view_hdf5_ctx));
2350: PetscCall(DMDestroy(&ctx->view_dm));
2351: for (PetscInt i = 0; i < NUM_FIELDS; i++) {
2352: PetscCall(VecScatterDestroy(&ctx->subsct[i]));
2353: PetscCall(VecDestroy(&ctx->subvec[i]));
2354: }
2356: PetscCall(TSDestroy(&ts));
2357: PetscCall(VecTaggerDestroy(&actx->refineTag));
2358: PetscCall(PetscFree(actx));
2359: PetscFunctionReturn(PETSC_SUCCESS);
2360: }
2362: int main(int argc, char **argv)
2363: {
2364: AppCtx ctx;
2366: PetscFunctionBeginUser;
2367: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2368: PetscCall(ProcessOptions(&ctx));
2369: PetscCall(PetscLogStageRegister("Setup", &SetupStage));
2370: PetscCall(PetscLogStageRegister("Solve", &SolveStage));
2371: if (ctx.test_restart) { /* Test sequences of save and loads */
2372: PetscMPIInt rank;
2374: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2376: /* test saving */
2377: ctx.load = PETSC_FALSE;
2378: ctx.save = PETSC_TRUE;
2379: /* sequential save */
2380: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential save\n"));
2381: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_seq_%d.h5", rank));
2382: PetscCall(Run(PETSC_COMM_SELF, &ctx));
2383: /* parallel save */
2384: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel save\n"));
2385: PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_par.h5"));
2386: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
2388: /* test loading */
2389: ctx.load = PETSC_TRUE;
2390: ctx.save = PETSC_FALSE;
2391: /* sequential and parallel runs from sequential save */
2392: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_seq_0.h5"));
2393: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from sequential save\n"));
2394: PetscCall(Run(PETSC_COMM_SELF, &ctx));
2395: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from sequential save\n"));
2396: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
2397: /* sequential and parallel runs from parallel save */
2398: PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_par.h5"));
2399: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from parallel save\n"));
2400: PetscCall(Run(PETSC_COMM_SELF, &ctx));
2401: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from parallel save\n"));
2402: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
2403: } else { /* Run the simulation */
2404: PetscCall(Run(PETSC_COMM_WORLD, &ctx));
2405: }
2406: PetscCall(PetscFree5(ctx.source_ctx->x0, ctx.source_ctx->w, ctx.source_ctx->k, ctx.source_ctx->p, ctx.source_ctx->r));
2407: PetscCall(PetscFree(ctx.source_ctx));
2408: PetscCall(PetscFinalize());
2409: return 0;
2410: }
2412: /*TEST
2414: testset:
2415: 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 -ic_num 3
2417: test:
2418: suffix: 0
2419: nsize: {{1 2}}
2420: args: -dm_refine 1 -lump {{0 1}} -exclude_potential_lte
2422: test:
2423: suffix: 0_split
2424: nsize: {{1 2}}
2425: args: -dm_refine 1 -split
2427: test:
2428: suffix: 0_3d
2429: nsize: {{1 2}}
2430: args: -dm_plex_box_faces 3,3,3 -dim 3 -dm_plex_dim 3 -lump {{0 1}}
2432: test:
2433: suffix: 0_dirk
2434: nsize: {{1 2}}
2435: args: -dm_refine 1 -ts_type dirk
2437: test:
2438: suffix: 0_dirk_mg
2439: nsize: {{1 2}}
2440: 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}}
2442: test:
2443: suffix: 0_dirk_fieldsplit
2444: nsize: {{1 2}}
2445: args: -dm_refine 1 -ts_type dirk -pc_type fieldsplit -pc_fieldsplit_type multiplicative -lump {{0 1}}
2447: test:
2448: requires: p4est
2449: suffix: 0_p4est
2450: nsize: {{1 2}}
2451: args: -dm_refine 1 -dm_plex_convert_type p4est -lump 0
2453: test:
2454: suffix: 0_periodic
2455: nsize: {{1 2}}
2456: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -lump {{0 1}}
2458: test:
2459: requires: p4est
2460: suffix: 0_p4est_periodic
2461: nsize: {{1 2}}
2462: args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -dm_plex_convert_type p4est -lump 0
2464: test:
2465: requires: p4est
2466: suffix: 0_p4est_mg
2467: nsize: {{1 2}}
2468: 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
2470: testset:
2471: requires: hdf5
2472: 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
2474: test:
2475: requires: !single
2476: suffix: restart
2477: nsize: {{1 2}separate output}
2478: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 0
2480: test:
2481: requires: triangle
2482: suffix: restart_simplex
2483: nsize: {{1 2}separate output}
2484: args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 1
2486: test:
2487: requires: !single
2488: suffix: restart_refonly
2489: nsize: {{1 2}separate output}
2490: args: -dm_refine 1 -dm_plex_simplex 0
2492: test:
2493: requires: triangle
2494: suffix: restart_simplex_refonly
2495: nsize: {{1 2}separate output}
2496: args: -dm_refine 1 -dm_plex_simplex 1
2498: test:
2499: suffix: annulus
2500: requires: exodusii
2501: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2
2503: test:
2504: suffix: hdf5_diagnostic
2505: requires: hdf5 exodusii
2506: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2 -monitor_hdf5 diagnostic.h5 -ic_num 3
2508: test:
2509: suffix: vtk_diagnostic
2510: requires: exodusii
2511: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2 -monitor_vtk 'diagnostic-%03d.vtu' -ic_num 3
2513: test:
2514: suffix: full_cdisc
2515: args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition selfp -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_mat_schur_complement_ainv_type blockdiag -fieldsplit_potential_ksp_type preonly -fieldsplit_potential_pc_type svd
2517: test:
2518: suffix: full_cdisc_split
2519: args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_pc_type gamg -split -monitor_norms
2521: test:
2522: suffix: full_cdisc_minres
2523: args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type diag -pc_fieldsplit_schur_precondition selfp -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_mat_schur_complement_ainv_type blockdiag -fieldsplit_potential_ksp_type preonly -fieldsplit_potential_pc_type svd -ksp_type minres
2525: TEST*/