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