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), &section));
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(&section));

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*/