Actual source code: ex30.c

      1: static char help[] = "Biological network from https://link.springer.com/article/10.1007/s42967-023-00297-3\n\n\n";

      3: #include <petscts.h>
      4: #include <petscsf.h>
      5: #include <petscdmplex.h>
      6: #include <petscdmplextransform.h>
      7: #include <petscdmforest.h>
      8: #include <petscviewerhdf5.h>
      9: #include <petscds.h>

     11: /*
     12:     Here we solve the system of PDEs on \Omega \in R^2:

     14:     * dC/dt - D^2 \Delta C - c^2 \nabla p \cross \nabla p + \alpha sqrt(||C||^2_F + eps)^(\gamma-2) C = 0
     15:     * - \nabla \cdot ((r + C) \nabla p) = S

     17:     where:
     18:       C = symmetric 2x2 conductivity tensor
     19:       p = potential
     20:       S = source

     22:     with natural boundary conditions on \partial\Omega:
     23:       \nabla C \cdot n  = 0
     24:       \nabla ((r + C)\nabla p) \cdot n  = 0

     26:     Parameters:
     27:       D = diffusion constant
     28:       c = activation parameter
     29:       \alpha = metabolic coefficient
     30:       \gamma = metabolic exponent
     31:       r, eps are regularization parameters

     33:     We use Lagrange elements for C_ij and P.
     34: */

     36: typedef enum _fieldidx {
     37:   C_FIELD_ID = 0,
     38:   P_FIELD_ID,
     39:   NUM_FIELDS
     40: } FieldIdx;

     42: typedef enum _constantidx {
     43:   R_ID = 0,
     44:   EPS_ID,
     45:   ALPHA_ID,
     46:   GAMMA_ID,
     47:   D_ID,
     48:   C2_ID,
     49:   FIXC_ID,
     50:   NUM_CONSTANTS
     51: } ConstantIdx;

     53: PetscLogStage SetupStage, SolveStage;

     55: #define NORM2C(c00, c01, c11) PetscSqr(c00) + 2 * PetscSqr(c01) + PetscSqr(c11)

     57: /* residual for C when tested against basis functions */
     58: static void C_0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
     59: {
     60:   const PetscReal   c2       = PetscRealPart(constants[C2_ID]);
     61:   const PetscReal   alpha    = PetscRealPart(constants[ALPHA_ID]);
     62:   const PetscReal   gamma    = PetscRealPart(constants[GAMMA_ID]);
     63:   const PetscReal   eps      = PetscRealPart(constants[EPS_ID]);
     64:   const PetscScalar gradp[]  = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
     65:   const PetscScalar crossp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[1] * gradp[1]};
     66:   const PetscScalar C00      = u[uOff[C_FIELD_ID]];
     67:   const PetscScalar C01      = u[uOff[C_FIELD_ID] + 1];
     68:   const PetscScalar C11      = u[uOff[C_FIELD_ID] + 2];
     69:   const PetscScalar norm     = NORM2C(C00, C01, C11) + eps;
     70:   const PetscScalar nexp     = (gamma - 2.0) / 2.0;
     71:   const PetscScalar fnorm    = PetscPowScalar(norm, nexp);

     73:   for (PetscInt k = 0; k < 3; k++) f0[k] = u_t[uOff[C_FIELD_ID] + k] - c2 * crossp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k];
     74: }

     76: /* Jacobian for C against C basis functions */
     77: static void JC_0_c0c0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
     78: {
     79:   const PetscReal   alpha  = PetscRealPart(constants[ALPHA_ID]);
     80:   const PetscReal   gamma  = PetscRealPart(constants[GAMMA_ID]);
     81:   const PetscReal   eps    = PetscRealPart(constants[EPS_ID]);
     82:   const PetscScalar C00    = u[uOff[C_FIELD_ID]];
     83:   const PetscScalar C01    = u[uOff[C_FIELD_ID] + 1];
     84:   const PetscScalar C11    = u[uOff[C_FIELD_ID] + 2];
     85:   const PetscScalar norm   = NORM2C(C00, C01, C11) + eps;
     86:   const PetscScalar nexp   = (gamma - 2.0) / 2.0;
     87:   const PetscScalar fnorm  = PetscPowScalar(norm, nexp);
     88:   const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0);
     89:   const PetscScalar dC[]   = {2 * C00, 4 * C01, 2 * C11};

     91:   for (PetscInt k = 0; k < 3; k++) {
     92:     for (PetscInt j = 0; j < 3; j++) J[k * 3 + j] = alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k];
     93:     J[k * 3 + k] += alpha * fnorm + u_tShift;
     94:   }
     95: }

     97: /* Jacobian for C against C basis functions and gradients of P basis functions */
     98: static void JC_0_c0p1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
     99: {
    100:   const PetscReal   c2      = PetscRealPart(constants[C2_ID]);
    101:   const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};

    103:   J[0] = -c2 * 2 * gradp[0];
    104:   J[1] = 0.0;
    105:   J[2] = -c2 * gradp[1];
    106:   J[3] = -c2 * gradp[0];
    107:   J[4] = 0.0;
    108:   J[5] = -c2 * 2 * gradp[1];
    109: }

    111: /* residual for C when tested against gradients of basis functions */
    112: static void C_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    113: {
    114:   const PetscReal D = PetscRealPart(constants[D_ID]);
    115:   for (PetscInt k = 0; k < 3; k++)
    116:     for (PetscInt d = 0; d < 2; d++) f1[k * 2 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 2 + d];
    117: }

    119: /* Jacobian for C against gradients of C basis functions */
    120: static void JC_1_c1c1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
    121: {
    122:   const PetscReal D = PetscRealPart(constants[D_ID]);
    123:   for (PetscInt k = 0; k < 3; k++)
    124:     for (PetscInt d = 0; d < 2; d++) J[k * (3 + 1) * 2 * 2 + d * 2 + d] = PetscSqr(D);
    125: }

    127: /* residual for P when tested against basis functions.
    128:    The source term always comes from the auxiliary vec because it needs to have zero mean */
    129: static void P_0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    130: {
    131:   PetscScalar S = a[aOff[P_FIELD_ID]];

    133:   f0[0] = -S;
    134: }

    136: static inline void QuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal x[2]);

    138: /* compute shift to make C positive definite */
    139: static inline PetscReal FIX_C(PetscScalar C00, PetscScalar C01, PetscScalar C11)
    140: {
    141: #if !PetscDefined(USE_COMPLEX)
    142:   PetscReal eigs[2], s = 0.0;

    144:   QuadraticRoots(1, -(C00 + C11), C00 * C11 - PetscSqr(C01), eigs);
    145:   if (eigs[0] < 0 || eigs[1] < 0) s = -PetscMin(eigs[0], eigs[1]) + PETSC_SMALL;
    146:   return s;
    147: #else
    148:   return 0.0;
    149: #endif
    150: }

    152: /* residual for P when tested against gradients of basis functions */
    153: static void P_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    154: {
    155:   const PetscReal   r       = PetscRealPart(constants[R_ID]);
    156:   const PetscScalar C00     = u[uOff[C_FIELD_ID]] + r;
    157:   const PetscScalar C01     = u[uOff[C_FIELD_ID] + 1];
    158:   const PetscScalar C10     = C01;
    159:   const PetscScalar C11     = u[uOff[C_FIELD_ID] + 2] + r;
    160:   const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
    161:   const PetscBool   fix_c   = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
    162:   const PetscScalar s       = fix_c ? FIX_C(C00, C01, C11) : 0.0;

    164:   f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1];
    165:   f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1];
    166: }

    168: /* Same as above for the P-only subproblem for initial conditions: the conductivity values come from the auxiliary vec */
    169: static void P_1_aux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    170: {
    171:   const PetscReal   r       = PetscRealPart(constants[R_ID]);
    172:   const PetscScalar C00     = a[aOff[C_FIELD_ID]] + r;
    173:   const PetscScalar C01     = a[aOff[C_FIELD_ID] + 1];
    174:   const PetscScalar C10     = C01;
    175:   const PetscScalar C11     = a[aOff[C_FIELD_ID] + 2] + r;
    176:   const PetscScalar gradp[] = {u_x[uOff_x[0]], u_x[uOff_x[0] + 1]};
    177:   const PetscBool   fix_c   = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0);
    178:   const PetscScalar s       = fix_c ? FIX_C(C00, C01, C11) : 0.0;

    180:   f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1];
    181:   f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1];
    182: }

    184: /* Jacobian for P against gradients of P basis functions */
    185: static void JP_1_p1p1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
    186: {
    187:   const PetscReal   r     = PetscRealPart(constants[R_ID]);
    188:   const PetscScalar C00   = u[uOff[C_FIELD_ID]] + r;
    189:   const PetscScalar C01   = u[uOff[C_FIELD_ID] + 1];
    190:   const PetscScalar C10   = C01;
    191:   const PetscScalar C11   = u[uOff[C_FIELD_ID] + 2] + r;
    192:   const PetscBool   fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
    193:   const PetscScalar s     = fix_c ? FIX_C(C00, C01, C11) : 0.0;

    195:   J[0] = C00 + s;
    196:   J[1] = C01;
    197:   J[2] = C10;
    198:   J[3] = C11 + s;
    199: }

    201: /* Same as above for the P-only subproblem for initial conditions */
    202: static void JP_1_p1p1_aux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
    203: {
    204:   const PetscReal   r     = PetscRealPart(constants[R_ID]);
    205:   const PetscScalar C00   = a[aOff[C_FIELD_ID]] + r;
    206:   const PetscScalar C01   = a[aOff[C_FIELD_ID] + 1];
    207:   const PetscScalar C10   = C01;
    208:   const PetscScalar C11   = a[aOff[C_FIELD_ID] + 2] + r;
    209:   const PetscBool   fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0);
    210:   const PetscScalar s     = fix_c ? FIX_C(C00, C01, C11) : 0.0;

    212:   J[0] = C00 + s;
    213:   J[1] = C01;
    214:   J[2] = C10;
    215:   J[3] = C11 + s;
    216: }

    218: /* Jacobian for P against gradients of P basis functions and C basis functions */
    219: static void JP_1_p1c0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[])
    220: {
    221:   const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};

    223:   J[0] = gradp[0];
    224:   J[1] = 0;
    225:   J[2] = gradp[1];
    226:   J[3] = gradp[0];
    227:   J[4] = 0;
    228:   J[5] = gradp[1];
    229: }

    231: /* the source term S(x) = exp(-500*||x - x0||^2) */
    232: static PetscErrorCode source_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    233: {
    234:   PetscReal *x0 = (PetscReal *)ctx;
    235:   PetscReal  n  = 0;

    237:   for (PetscInt d = 0; d < dim; ++d) n += (x[d] - x0[d]) * (x[d] - x0[d]);
    238:   u[0] = PetscExpReal(-500 * n);
    239:   return PETSC_SUCCESS;
    240: }

    242: static PetscErrorCode source_1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    243: {
    244:   PetscScalar     ut[1];
    245:   const PetscReal x0[] = {0.25, 0.25};
    246:   const PetscReal x1[] = {0.75, 0.75};

    248:   PetscCall(source_0(dim, time, x, Nf, ut, (void *)x0));
    249:   PetscCall(source_0(dim, time, x, Nf, u, (void *)x1));
    250:   u[0] += ut[0];
    251:   return PETSC_SUCCESS;
    252: }

    254: /* functionals to be integrated: average -> \int_\Omega u dx */
    255: static void average(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
    256: {
    257:   obj[0] = u[uOff[P_FIELD_ID]];
    258: }

    260: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
    261: static inline void QuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal x[2])
    262: {
    263:   PetscReal delta = PetscMax(b * b - 4 * a * c, 0); /* eigenvalues symmetric matrix */
    264:   PetscReal temp  = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(delta));

    266:   x[0] = temp / a;
    267:   x[1] = c / temp;
    268: }

    270: /* functionals to be integrated: energy -> D^2/2 * ||\nabla C||^2 + c^2\nabla p * (r + C) * \nabla p + \alpha/ \gamma * ||C||^\gamma */
    271: static void energy(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
    272: {
    273:   const PetscReal   D         = PetscRealPart(constants[D_ID]);
    274:   const PetscReal   c2        = PetscRealPart(constants[C2_ID]);
    275:   const PetscReal   r         = PetscRealPart(constants[R_ID]);
    276:   const PetscReal   alpha     = PetscRealPart(constants[ALPHA_ID]);
    277:   const PetscReal   gamma     = PetscRealPart(constants[GAMMA_ID]);
    278:   const PetscScalar C00       = u[uOff[C_FIELD_ID]];
    279:   const PetscScalar C01       = u[uOff[C_FIELD_ID] + 1];
    280:   const PetscScalar C10       = C01;
    281:   const PetscScalar C11       = u[uOff[C_FIELD_ID] + 2];
    282:   const PetscScalar gradp[]   = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};
    283:   const PetscScalar gradC00[] = {u_x[uOff_x[C_FIELD_ID] + 0], u_x[uOff_x[C_FIELD_ID] + 1]};
    284:   const PetscScalar gradC01[] = {u_x[uOff_x[C_FIELD_ID] + 2], u_x[uOff_x[C_FIELD_ID] + 3]};
    285:   const PetscScalar gradC11[] = {u_x[uOff_x[C_FIELD_ID] + 4], u_x[uOff_x[C_FIELD_ID] + 5]};
    286:   const PetscScalar normC     = NORM2C(C00, C01, C11);
    287:   const PetscScalar normgradC = NORM2C(gradC00[0], gradC01[0], gradC11[0]) + NORM2C(gradC00[1], gradC01[1], gradC11[1]);
    288:   const PetscScalar nexp      = gamma / 2.0;

    290:   const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC;
    291:   const PetscScalar t1 = c2 * (gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1]));
    292:   const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp);

    294:   obj[0] = t0 + t1 + t2;
    295: }

    297: /* functionals to be integrated: ellipticity_fail -> 0 means C+r is elliptic at quadrature point, otherwise it returns the absolute value of the most negative eigenvalue */
    298: static void ellipticity_fail(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
    299: {
    300:   const PetscReal r   = PetscRealPart(constants[R_ID]);
    301:   const PetscReal C00 = PetscRealPart(u[uOff[C_FIELD_ID]] + r);
    302:   const PetscReal C01 = PetscRealPart(u[uOff[C_FIELD_ID] + 1]);
    303:   const PetscReal C11 = PetscRealPart(u[uOff[C_FIELD_ID] + 2] + r);

    305:   PetscReal eigs[2];
    306:   QuadraticRoots(1, -(C00 + C11), C00 * C11 - PetscSqr(C01), eigs);
    307:   if (eigs[0] < 0 || eigs[1] < 0) obj[0] = -PetscMin(eigs[0], eigs[1]);
    308:   else obj[0] = 0.0;
    309: }

    311: /* initial conditions for C: eq. 16 */
    312: static PetscErrorCode initial_conditions_C_0(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
    313: {
    314:   u[0] = 1;
    315:   u[1] = 0;
    316:   u[2] = 1;
    317:   return PETSC_SUCCESS;
    318: }

    320: /* initial conditions for C: eq. 17 */
    321: static PetscErrorCode initial_conditions_C_1(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
    322: {
    323:   const PetscReal x = xx[0];
    324:   const PetscReal y = xx[1];

    326:   u[0] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
    327:   u[1] = 0;
    328:   u[2] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y));
    329:   return PETSC_SUCCESS;
    330: }

    332: /* initial conditions for C: eq. 18 */
    333: static PetscErrorCode initial_conditions_C_2(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
    334: {
    335:   u[0] = 0;
    336:   u[1] = 0;
    337:   u[2] = 0;
    338:   return PETSC_SUCCESS;
    339: }

    341: /* functionals to be sampled: C * \grad p */
    342: static void flux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
    343: {
    344:   const PetscScalar C00     = u[uOff[C_FIELD_ID]];
    345:   const PetscScalar C01     = u[uOff[C_FIELD_ID] + 1];
    346:   const PetscScalar C10     = C01;
    347:   const PetscScalar C11     = u[uOff[C_FIELD_ID] + 2];
    348:   const PetscScalar gradp[] = {u_x[uOff_x[P_FIELD_ID]], u_x[uOff_x[P_FIELD_ID] + 1]};

    350:   f[0] = C00 * gradp[0] + C01 * gradp[1];
    351:   f[1] = C10 * gradp[0] + C11 * gradp[1];
    352: }

    354: /* functionals to be sampled: ||C|| */
    355: static void normc(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
    356: {
    357:   const PetscScalar C00 = u[uOff[C_FIELD_ID]];
    358:   const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1];
    359:   const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2];

    361:   f[0] = PetscSqrtScalar(NORM2C(C00, C01, C11));
    362: }

    364: /* functionals to be sampled: zero */
    365: static void zero(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
    366: {
    367:   f[0] = 0.0;
    368: }

    370: /* functions to be sampled: zero function */
    371: static PetscErrorCode zerof(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
    372: {
    373:   for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
    374:   return PETSC_SUCCESS;
    375: }

    377: /* functions to be sampled: constant function */
    378: static PetscErrorCode constantf(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
    379: {
    380:   PetscInt d;
    381:   for (d = 0; d < Nc; ++d) u[d] = 1.0;
    382:   return PETSC_SUCCESS;
    383: }

    385: /* application context: customizable parameters */
    386: typedef struct {
    387:   PetscReal r;
    388:   PetscReal eps;
    389:   PetscReal alpha;
    390:   PetscReal gamma;
    391:   PetscReal D;
    392:   PetscReal c;
    393:   PetscInt  ic_num;
    394:   PetscInt  source_num;
    395:   PetscReal x0[2];
    396:   PetscBool lump;
    397:   PetscBool amr;
    398:   PetscBool load;
    399:   char      load_filename[PETSC_MAX_PATH_LEN];
    400:   PetscBool save;
    401:   char      save_filename[PETSC_MAX_PATH_LEN];
    402:   PetscInt  save_every;
    403:   PetscBool test_restart;
    404:   PetscBool ellipticity;
    405:   PetscInt  fix_c;
    406: } AppCtx;

    408: /* process command line options */
    409: static PetscErrorCode ProcessOptions(AppCtx *options)
    410: {
    411:   PetscInt dim = PETSC_STATIC_ARRAY_LENGTH(options->x0);

    413:   PetscFunctionBeginUser;
    414:   options->r            = 1.e-1;
    415:   options->eps          = 1.e-3;
    416:   options->alpha        = 0.75;
    417:   options->gamma        = 0.75;
    418:   options->c            = 5;
    419:   options->D            = 1.e-2;
    420:   options->ic_num       = 0;
    421:   options->source_num   = 0;
    422:   options->x0[0]        = 0.25;
    423:   options->x0[1]        = 0.25;
    424:   options->lump         = PETSC_FALSE;
    425:   options->amr          = PETSC_FALSE;
    426:   options->load         = PETSC_FALSE;
    427:   options->save         = PETSC_FALSE;
    428:   options->save_every   = -1;
    429:   options->test_restart = PETSC_FALSE;
    430:   options->ellipticity  = PETSC_FALSE;
    431:   options->fix_c        = 1; /* 1 means only Jac, 2 means function and Jac */

    433:   PetscOptionsBegin(PETSC_COMM_WORLD, "", __FILE__, "DMPLEX");
    434:   PetscCall(PetscOptionsReal("-alpha", "alpha", __FILE__, options->alpha, &options->alpha, NULL));
    435:   PetscCall(PetscOptionsReal("-gamma", "gamma", __FILE__, options->gamma, &options->gamma, NULL));
    436:   PetscCall(PetscOptionsReal("-c", "c", __FILE__, options->c, &options->c, NULL));
    437:   PetscCall(PetscOptionsReal("-d", "D", __FILE__, options->D, &options->D, NULL));
    438:   PetscCall(PetscOptionsReal("-eps", "eps", __FILE__, options->eps, &options->eps, NULL));
    439:   PetscCall(PetscOptionsReal("-r", "r", __FILE__, options->r, &options->r, NULL));
    440:   PetscCall(PetscOptionsRealArray("-x0", "x0", __FILE__, options->x0, &dim, NULL));
    441:   PetscCall(PetscOptionsInt("-ic_num", "ic_num", __FILE__, options->ic_num, &options->ic_num, NULL));
    442:   PetscCall(PetscOptionsInt("-source_num", "source_num", __FILE__, options->source_num, &options->source_num, NULL));
    443:   PetscCall(PetscOptionsBool("-lump", "use mass lumping", __FILE__, options->lump, &options->lump, NULL));
    444:   PetscCall(PetscOptionsInt("-fix_c", "shift conductivity to always be positive semi-definite", __FILE__, options->fix_c, &options->fix_c, NULL));
    445:   PetscCall(PetscOptionsBool("-amr", "use adaptive mesh refinement", __FILE__, options->amr, &options->amr, NULL));
    446:   PetscCall(PetscOptionsBool("-test_restart", "test restarting files", __FILE__, options->test_restart, &options->test_restart, NULL));
    447:   if (!options->test_restart) {
    448:     PetscCall(PetscOptionsString("-load", "filename with data to be loaded for restarting", __FILE__, options->load_filename, options->load_filename, PETSC_MAX_PATH_LEN, &options->load));
    449:     PetscCall(PetscOptionsString("-save", "filename with data to be saved for restarting", __FILE__, options->save_filename, options->save_filename, PETSC_MAX_PATH_LEN, &options->save));
    450:     if (options->save) PetscCall(PetscOptionsInt("-save_every", "save every n timestep (-1 saves only the last)", __FILE__, options->save_every, &options->save_every, NULL));
    451:   }
    452:   PetscCall(PetscOptionsBool("-monitor_ellipticity", "Dump locations of ellipticity violation", __FILE__, options->ellipticity, &options->ellipticity, NULL));
    453:   PetscOptionsEnd();
    454:   PetscFunctionReturn(PETSC_SUCCESS);
    455: }

    457: static PetscErrorCode SaveToFile(DM dm, Vec u, const char *filename)
    458: {
    459: #if defined(PETSC_HAVE_HDF5)
    460:   PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
    461:   PetscViewer       viewer;
    462:   DM                cdm       = dm;
    463:   PetscInt          numlevels = 0;

    465:   PetscFunctionBeginUser;
    466:   while (cdm) {
    467:     numlevels++;
    468:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    469:   }
    470:   /* Cannot be set programmatically */
    471:   PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_view_hdf5_storage_version 3.0.0"));
    472:   PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), filename, FILE_MODE_WRITE, &viewer));
    473:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numlevels", PETSC_INT, &numlevels));
    474:   PetscCall(PetscViewerPushFormat(viewer, format));
    475:   for (PetscInt level = numlevels - 1; level >= 0; level--) {
    476:     PetscInt    cc, rr;
    477:     PetscBool   isRegular, isUniform;
    478:     const char *dmname;
    479:     char        groupname[PETSC_MAX_PATH_LEN];

    481:     PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
    482:     PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
    483:     PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
    484:     PetscCall(DMGetCoarsenLevel(dm, &cc));
    485:     PetscCall(DMGetRefineLevel(dm, &rr));
    486:     PetscCall(DMPlexGetRegularRefinement(dm, &isRegular));
    487:     PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
    488:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "meshname", PETSC_STRING, dmname));
    489:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refinelevel", PETSC_INT, &rr));
    490:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, &cc));
    491:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refRegular", PETSC_BOOL, &isRegular));
    492:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refUniform", PETSC_BOOL, &isUniform));
    493:     PetscCall(DMPlexTopologyView(dm, viewer));
    494:     PetscCall(DMPlexLabelsView(dm, viewer));
    495:     PetscCall(DMPlexCoordinatesView(dm, viewer));
    496:     PetscCall(DMPlexSectionView(dm, viewer, NULL));
    497:     if (level == numlevels - 1) {
    498:       PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
    499:       PetscCall(DMPlexGlobalVectorView(dm, viewer, NULL, u));
    500:     }
    501:     if (level) {
    502:       PetscInt        cStart, cEnd, ccStart, ccEnd, cpStart;
    503:       DMPolytopeType  ct;
    504:       DMPlexTransform tr;
    505:       DM              sdm;
    506:       PetscScalar    *array;
    507:       PetscSection    section;
    508:       Vec             map;
    509:       IS              gis;
    510:       const PetscInt *gidx;

    512:       PetscCall(DMGetCoarseDM(dm, &cdm));
    513:       PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    514:       PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section));
    515:       PetscCall(PetscSectionSetChart(section, cStart, cEnd));
    516:       for (PetscInt c = cStart; c < cEnd; c++) PetscCall(PetscSectionSetDof(section, c, 1));
    517:       PetscCall(PetscSectionSetUp(section));

    519:       PetscCall(DMClone(dm, &sdm));
    520:       PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
    521:       PetscCall(PetscObjectSetName((PetscObject)section, "pdm_section"));
    522:       PetscCall(DMSetLocalSection(sdm, section));
    523:       PetscCall(PetscSectionDestroy(&section));

    525:       PetscCall(DMGetLocalVector(sdm, &map));
    526:       PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
    527:       PetscCall(VecGetArray(map, &array));
    528:       PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
    529:       PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
    530:       PetscCall(DMPlexTransformSetDM(tr, cdm));
    531:       PetscCall(DMPlexTransformSetFromOptions(tr));
    532:       PetscCall(DMPlexTransformSetUp(tr));
    533:       PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
    534:       PetscCall(DMPlexGetChart(cdm, &cpStart, NULL));
    535:       PetscCall(DMPlexCreatePointNumbering(cdm, &gis));
    536:       PetscCall(ISGetIndices(gis, &gidx));
    537:       for (PetscInt c = ccStart; c < ccEnd; c++) {
    538:         PetscInt       *rsize, *rcone, *rornt, Nt;
    539:         DMPolytopeType *rct;
    540:         PetscInt        gnum = gidx[c - cpStart] >= 0 ? gidx[c - cpStart] : -(gidx[c - cpStart] + 1);

    542:         PetscCall(DMPlexGetCellType(cdm, c, &ct));
    543:         PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nt, &rct, &rsize, &rcone, &rornt));
    544:         for (PetscInt r = 0; r < rsize[Nt - 1]; ++r) {
    545:           PetscInt pNew;

    547:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[Nt - 1], c, r, &pNew));
    548:           array[pNew - cStart] = gnum;
    549:         }
    550:       }
    551:       PetscCall(ISRestoreIndices(gis, &gidx));
    552:       PetscCall(ISDestroy(&gis));
    553:       PetscCall(VecRestoreArray(map, &array));
    554:       PetscCall(DMPlexTransformDestroy(&tr));
    555:       PetscCall(DMPlexSectionView(dm, viewer, sdm));
    556:       PetscCall(DMPlexLocalVectorView(dm, viewer, sdm, map));
    557:       PetscCall(DMRestoreLocalVector(sdm, &map));
    558:       PetscCall(DMDestroy(&sdm));
    559:     }
    560:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    561:     PetscCall(DMGetCoarseDM(dm, &dm));
    562:   }
    563:   PetscCall(PetscViewerPopFormat(viewer));
    564:   PetscCall(PetscViewerDestroy(&viewer));
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: #else
    567:   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
    568: #endif
    569: }

    571: static PetscErrorCode LoadFromFile(MPI_Comm comm, const char *filename, DM *odm)
    572: {
    573: #if defined(PETSC_HAVE_HDF5)
    574:   PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
    575:   PetscViewer       viewer;
    576:   DM                dm, cdm = NULL;
    577:   PetscSF           sfXC      = NULL;
    578:   PetscInt          numlevels = -1;

    580:   PetscFunctionBeginUser;
    581:   PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer));
    582:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numlevels", PETSC_INT, NULL, &numlevels));
    583:   PetscCall(PetscViewerPushFormat(viewer, format));
    584:   for (PetscInt level = 0; level < numlevels; level++) {
    585:     char             groupname[PETSC_MAX_PATH_LEN], *dmname;
    586:     PetscSF          sfXB, sfBC, sfG;
    587:     PetscPartitioner part;
    588:     PetscInt         rr, cc;
    589:     PetscBool        isRegular, isUniform;

    591:     PetscCall(DMCreate(comm, &dm));
    592:     PetscCall(DMSetType(dm, DMPLEX));
    593:     PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level));
    594:     PetscCall(PetscViewerHDF5PushGroup(viewer, groupname));
    595:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "meshname", PETSC_STRING, NULL, &dmname));
    596:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refinelevel", PETSC_INT, NULL, &rr));
    597:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, NULL, &cc));
    598:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refRegular", PETSC_BOOL, NULL, &isRegular));
    599:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refUniform", PETSC_BOOL, NULL, &isUniform));
    600:     PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
    601:     PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
    602:     PetscCall(DMPlexLabelsLoad(dm, viewer, sfXB));
    603:     PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXB));
    604:     PetscCall(DMPlexGetPartitioner(dm, &part));
    605:     if (!level) { /* partition the coarse level only */
    606:       PetscCall(PetscPartitionerSetFromOptions(part));
    607:     } else { /* propagate partitioning information from coarser to finer level */
    608:       DM           sdm;
    609:       Vec          map;
    610:       PetscSF      sf;
    611:       PetscLayout  clayout;
    612:       PetscScalar *array;
    613:       PetscInt    *cranks_leaf, *cranks_root, *npoints, *points, *ranks, *starts, *gidxs;
    614:       PetscInt     nparts, cStart, cEnd, nr, ccStart, ccEnd, cpStart, cpEnd;
    615:       PetscMPIInt  size, rank;

    617:       PetscCall(DMClone(dm, &sdm));
    618:       PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm"));
    619:       PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXB, NULL, &sf));
    620:       PetscCall(DMGetLocalVector(sdm, &map));
    621:       PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map"));
    622:       PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, sf, map));

    624:       PetscCallMPI(MPI_Comm_size(comm, &size));
    625:       PetscCallMPI(MPI_Comm_rank(comm, &rank));
    626:       nparts = size;
    627:       PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    628:       PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd));
    629:       PetscCall(DMPlexGetChart(cdm, &cpStart, &cpEnd));
    630:       PetscCall(PetscCalloc1(nparts, &npoints));
    631:       PetscCall(PetscMalloc4(cEnd - cStart, &points, cEnd - cStart, &ranks, nparts + 1, &starts, cEnd - cStart, &gidxs));
    632:       PetscCall(PetscSFGetGraph(sfXC, &nr, NULL, NULL, NULL));
    633:       PetscCall(PetscMalloc2(cpEnd - cpStart, &cranks_leaf, nr, &cranks_root));
    634:       for (PetscInt c = 0; c < cpEnd - cpStart; c++) cranks_leaf[c] = rank;
    635:       PetscCall(PetscSFReduceBegin(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));
    636:       PetscCall(PetscSFReduceEnd(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE));

    638:       PetscCall(VecGetArray(map, &array));
    639:       for (PetscInt c = 0; c < cEnd - cStart; c++) gidxs[c] = (PetscInt)PetscRealPart(array[c]);
    640:       PetscCall(VecRestoreArray(map, &array));

    642:       PetscCall(PetscLayoutCreate(comm, &clayout));
    643:       PetscCall(PetscLayoutSetLocalSize(clayout, nr));
    644:       PetscCall(PetscSFSetGraphLayout(sf, clayout, cEnd - cStart, NULL, PETSC_OWN_POINTER, gidxs));
    645:       PetscCall(PetscLayoutDestroy(&clayout));

    647:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
    648:       PetscCall(PetscSFBcastEnd(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE));
    649:       PetscCall(PetscSFDestroy(&sf));
    650:       PetscCall(PetscFree2(cranks_leaf, cranks_root));
    651:       for (PetscInt c = 0; c < cEnd - cStart; c++) npoints[ranks[c]]++;

    653:       starts[0] = 0;
    654:       for (PetscInt c = 0; c < nparts; c++) starts[c + 1] = starts[c] + npoints[c];
    655:       for (PetscInt c = 0; c < cEnd - cStart; c++) points[starts[ranks[c]]++] = c;
    656:       PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
    657:       PetscCall(PetscPartitionerShellSetPartition(part, nparts, npoints, points));
    658:       PetscCall(PetscFree(npoints));
    659:       PetscCall(PetscFree4(points, ranks, starts, gidxs));
    660:       PetscCall(DMRestoreLocalVector(sdm, &map));
    661:       PetscCall(DMDestroy(&sdm));
    662:     }
    663:     PetscCall(PetscSFDestroy(&sfXC));
    664:     PetscCall(DMPlexDistribute(dm, 0, &sfBC, odm));
    665:     if (*odm) {
    666:       PetscCall(DMDestroy(&dm));
    667:       dm   = *odm;
    668:       *odm = NULL;
    669:       PetscCall(PetscObjectSetName((PetscObject)dm, dmname));
    670:     }
    671:     if (sfBC) PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
    672:     else {
    673:       PetscCall(PetscObjectReference((PetscObject)sfXB));
    674:       sfXC = sfXB;
    675:     }
    676:     PetscCall(PetscSFDestroy(&sfXB));
    677:     PetscCall(PetscSFDestroy(&sfBC));
    678:     PetscCall(DMSetCoarsenLevel(dm, cc));
    679:     PetscCall(DMSetRefineLevel(dm, rr));
    680:     PetscCall(DMPlexSetRegularRefinement(dm, isRegular));
    681:     PetscCall(DMPlexSetRefinementUniform(dm, isUniform));
    682:     PetscCall(DMPlexSectionLoad(dm, viewer, NULL, sfXC, &sfG, NULL));
    683:     if (level == numlevels - 1) {
    684:       Vec u;

    686:       PetscCall(DMGetNamedGlobalVector(dm, "solution_", &u));
    687:       PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
    688:       PetscCall(DMPlexGlobalVectorLoad(dm, viewer, NULL, sfG, u));
    689:       PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &u));
    690:     }
    691:     PetscCall(PetscFree(dmname));
    692:     PetscCall(PetscSFDestroy(&sfG));
    693:     PetscCall(DMSetCoarseDM(dm, cdm));
    694:     PetscCall(DMDestroy(&cdm));
    695:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    696:     cdm = dm;
    697:   }
    698:   *odm = dm;
    699:   PetscCall(PetscViewerPopFormat(viewer));
    700:   PetscCall(PetscViewerDestroy(&viewer));
    701:   PetscCall(PetscSFDestroy(&sfXC));
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: #else
    704:   SETERRQ(comm, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5");
    705: #endif
    706: }

    708: /* Project source function and make it zero-mean */
    709: static PetscErrorCode ProjectSource(DM dm, PetscReal time, AppCtx *ctx)
    710: {
    711:   PetscInt    id = C_FIELD_ID;
    712:   DM          dmAux;
    713:   Vec         u, lu;
    714:   IS          is;
    715:   void       *ctxs[NUM_FIELDS];
    716:   PetscScalar vals[NUM_FIELDS];
    717:   PetscDS     ds;
    718:   PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);

    720:   PetscFunctionBeginUser;
    721:   switch (ctx->source_num) {
    722:   case 0:
    723:     funcs[P_FIELD_ID] = source_0;
    724:     ctxs[P_FIELD_ID]  = ctx->x0;
    725:     break;
    726:   case 1:
    727:     funcs[P_FIELD_ID] = source_1;
    728:     ctxs[P_FIELD_ID]  = NULL;
    729:     break;
    730:   default:
    731:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown source");
    732:   }
    733:   funcs[C_FIELD_ID] = zerof;
    734:   ctxs[C_FIELD_ID]  = NULL;
    735:   PetscCall(DMGetDS(dm, &ds));
    736:   PetscCall(DMGetGlobalVector(dm, &u));
    737:   PetscCall(DMProjectFunction(dm, time, funcs, ctxs, INSERT_ALL_VALUES, u));
    738:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
    739:   PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
    740:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
    741:   PetscCall(VecShift(u, -vals[P_FIELD_ID]));
    742:   PetscCall(DMCreateSubDM(dm, 1, &id, &is, NULL));
    743:   PetscCall(VecISSet(u, is, 0));
    744:   PetscCall(ISDestroy(&is));

    746:   /* Attach source vector as auxiliary vector:
    747:      Use a different DM to break ref cycles */
    748:   PetscCall(DMClone(dm, &dmAux));
    749:   PetscCall(DMCopyDisc(dm, dmAux));
    750:   PetscCall(DMCreateLocalVector(dmAux, &lu));
    751:   PetscCall(DMDestroy(&dmAux));
    752:   PetscCall(DMGlobalToLocal(dm, u, INSERT_VALUES, lu));
    753:   PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, lu));
    754:   PetscCall(VecViewFromOptions(lu, NULL, "-aux_view"));
    755:   PetscCall(VecDestroy(&lu));
    756:   PetscCall(DMRestoreGlobalVector(dm, &u));
    757:   PetscFunctionReturn(PETSC_SUCCESS);
    758: }

    760: /* callback for the creation of the potential null space */
    761: static PetscErrorCode CreatePotentialNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
    762: {
    763:   Vec vec;
    764:   PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zerof};

    766:   PetscFunctionBeginUser;
    767:   funcs[nfield] = constantf;
    768:   PetscCall(DMCreateGlobalVector(dm, &vec));
    769:   PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
    770:   PetscCall(VecNormalize(vec, NULL));
    771:   PetscCall(PetscObjectSetName((PetscObject)vec, "Potential Null Space"));
    772:   PetscCall(VecViewFromOptions(vec, NULL, "-potential_nullspace_view"));
    773:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
    774:   /* break ref cycles */
    775:   PetscCall(VecSetDM(vec, NULL));
    776:   PetscCall(VecDestroy(&vec));
    777:   PetscFunctionReturn(PETSC_SUCCESS);
    778: }

    780: static PetscErrorCode DMGetLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
    781: {
    782:   PetscBool has;

    784:   PetscFunctionBeginUser;
    785:   if (local) {
    786:     PetscCall(DMHasNamedLocalVector(dm, "lumped_mass", &has));
    787:     PetscCall(DMGetNamedLocalVector(dm, "lumped_mass", lumped_mass));
    788:   } else {
    789:     PetscCall(DMHasNamedGlobalVector(dm, "lumped_mass", &has));
    790:     PetscCall(DMGetNamedGlobalVector(dm, "lumped_mass", lumped_mass));
    791:   }
    792:   if (!has) {
    793:     Vec w;
    794:     IS  is;

    796:     PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
    797:     if (!is) {
    798:       PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};

    800:       PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &is, NULL));
    801:       PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)is));
    802:       PetscCall(PetscObjectDereference((PetscObject)is));
    803:     }
    804:     if (local) {
    805:       Vec w2, wg;

    807:       PetscCall(DMCreateMassMatrixLumped(dm, &w, NULL));
    808:       PetscCall(DMGetGlobalVector(dm, &wg));
    809:       PetscCall(DMGetLocalVector(dm, &w2));
    810:       PetscCall(VecSet(w2, 0.0));
    811:       PetscCall(VecSet(wg, 1.0));
    812:       PetscCall(VecISSet(wg, is, 0.0));
    813:       PetscCall(DMGlobalToLocal(dm, wg, INSERT_VALUES, w2));
    814:       PetscCall(VecPointwiseMult(w, w, w2));
    815:       PetscCall(DMRestoreGlobalVector(dm, &wg));
    816:       PetscCall(DMRestoreLocalVector(dm, &w2));
    817:     } else {
    818:       PetscCall(DMCreateMassMatrixLumped(dm, NULL, &w));
    819:       PetscCall(VecISSet(w, is, 0.0));
    820:     }
    821:     PetscCall(VecCopy(w, *lumped_mass));
    822:     PetscCall(VecDestroy(&w));
    823:   }
    824:   PetscFunctionReturn(PETSC_SUCCESS);
    825: }

    827: static PetscErrorCode DMRestoreLumpedMass(DM dm, PetscBool local, Vec *lumped_mass)
    828: {
    829:   PetscFunctionBeginUser;
    830:   if (local) PetscCall(DMRestoreNamedLocalVector(dm, "lumped_mass", lumped_mass));
    831:   else PetscCall(DMRestoreNamedGlobalVector(dm, "lumped_mass", lumped_mass));
    832:   PetscFunctionReturn(PETSC_SUCCESS);
    833: }

    835: /* callbacks for lumped mass matrix residual and Jacobian */
    836: static PetscErrorCode DMPlexTSComputeIFunctionFEM_Lumped(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
    837: {
    838:   Vec work, local_lumped_mass;

    840:   PetscFunctionBeginUser;
    841:   PetscCall(DMGetLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
    842:   PetscCall(DMGetLocalVector(dm, &work));
    843:   PetscCall(VecSet(work, 0.0));
    844:   PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, work, locF, user));
    845:   PetscCall(VecPointwiseMult(work, locX_t, local_lumped_mass));
    846:   PetscCall(VecAXPY(locF, 1.0, work));
    847:   PetscCall(DMRestoreLocalVector(dm, &work));
    848:   PetscCall(DMRestoreLumpedMass(dm, PETSC_TRUE, &local_lumped_mass));
    849:   PetscFunctionReturn(PETSC_SUCCESS);
    850: }

    852: static PetscErrorCode DMPlexTSComputeIJacobianFEM_Lumped(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
    853: {
    854:   Vec lumped_mass, work;

    856:   PetscFunctionBeginUser;
    857:   // XXX CHECK DIRK
    858:   PetscCall(DMGetLumpedMass(dm, PETSC_FALSE, &lumped_mass));
    859:   PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, 0.0, Jac, JacP, user));
    860:   PetscCall(DMGetGlobalVector(dm, &work));
    861:   PetscCall(VecAXPBY(work, X_tShift, 0.0, lumped_mass));
    862:   PetscCall(MatDiagonalSet(JacP, work, ADD_VALUES));
    863:   PetscCall(DMRestoreGlobalVector(dm, &work));
    864:   PetscCall(DMRestoreLumpedMass(dm, PETSC_FALSE, &lumped_mass));
    865:   PetscFunctionReturn(PETSC_SUCCESS);
    866: }

    868: /* customize residuals and Jacobians */
    869: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
    870: {
    871:   PetscDS     ds;
    872:   PetscInt    cdim, dim;
    873:   PetscScalar constants[NUM_CONSTANTS];

    875:   PetscFunctionBeginUser;
    876:   constants[R_ID]     = ctx->r;
    877:   constants[EPS_ID]   = ctx->eps;
    878:   constants[ALPHA_ID] = ctx->alpha;
    879:   constants[GAMMA_ID] = ctx->gamma;
    880:   constants[D_ID]     = ctx->D;
    881:   constants[C2_ID]    = ctx->c * ctx->c;
    882:   constants[FIXC_ID]  = ctx->fix_c;

    884:   PetscCall(DMGetDimension(dm, &dim));
    885:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    886:   PetscCheck(dim == 2 && cdim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2D meshes");
    887:   PetscCall(DMGetDS(dm, &ds));
    888:   PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
    889:   PetscCall(PetscDSSetImplicit(ds, C_FIELD_ID, PETSC_TRUE));
    890:   PetscCall(PetscDSSetImplicit(ds, P_FIELD_ID, PETSC_TRUE));
    891:   PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
    892:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
    893:   PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0, C_1));
    894:   PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1));
    895:   PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0, NULL, NULL, JC_1_c1c1));
    896:   PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1, NULL, NULL));
    897:   PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0, NULL));
    898:   PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1));

    900:   /* Attach potential nullspace */
    901:   PetscCall(DMSetNullSpaceConstructor(dm, P_FIELD_ID, CreatePotentialNullSpace));

    903:   /* Attach source function as auxiliary vector */
    904:   PetscCall(ProjectSource(dm, 0, ctx));

    906:   /* Add callbacks */
    907:   if (ctx->lump) {
    908:     PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM_Lumped, NULL));
    909:     PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM_Lumped, NULL));
    910:   } else {
    911:     PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, NULL));
    912:     PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, NULL));
    913:   }
    914:   /* This is not really needed because we use Neumann boundaries */
    915:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL));
    916:   PetscFunctionReturn(PETSC_SUCCESS);
    917: }

    919: /* setup discrete spaces and residuals */
    920: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
    921: {
    922:   DM           plex, cdm = dm;
    923:   PetscFE      feC, feP;
    924:   PetscBool    simplex;
    925:   PetscInt     dim;
    926:   MPI_Comm     comm = PetscObjectComm((PetscObject)dm);
    927:   MatNullSpace nsp;

    929:   PetscFunctionBeginUser;
    930:   PetscCall(DMGetDimension(dm, &dim));

    932:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    933:   PetscCall(DMPlexIsSimplex(plex, &simplex));
    934:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &simplex, 1, MPIU_BOOL, MPI_LOR, comm));
    935:   PetscCall(DMDestroy(&plex));

    937:   /* We model Cij with Cij = Cji -> dim*(dim+1)/2 components */
    938:   PetscCall(PetscFECreateDefault(comm, dim, (dim * (dim + 1)) / 2, simplex, "c_", -1, &feC));
    939:   PetscCall(PetscObjectSetName((PetscObject)feC, "conductivity"));
    940:   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "p_", -1, &feP));
    941:   PetscCall(PetscObjectSetName((PetscObject)feP, "potential"));
    942:   PetscCall(MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, &nsp));
    943:   PetscCall(PetscObjectCompose((PetscObject)feP, "nullspace", (PetscObject)nsp));
    944:   PetscCall(MatNullSpaceDestroy(&nsp));
    945:   PetscCall(PetscFECopyQuadrature(feP, feC));

    947:   PetscCall(DMSetNumFields(dm, 2));
    948:   PetscCall(DMSetField(dm, C_FIELD_ID, NULL, (PetscObject)feC));
    949:   PetscCall(DMSetField(dm, P_FIELD_ID, NULL, (PetscObject)feP));
    950:   PetscCall(PetscFEDestroy(&feC));
    951:   PetscCall(PetscFEDestroy(&feP));
    952:   PetscCall(DMCreateDS(dm));
    953:   while (cdm) {
    954:     PetscCall(DMCopyDisc(dm, cdm));
    955:     PetscCall(SetupProblem(cdm, ctx));
    956:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    957:   }
    958:   PetscFunctionReturn(PETSC_SUCCESS);
    959: }

    961: /* Create mesh by command line options */
    962: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
    963: {
    964:   PetscFunctionBeginUser;
    965:   if (ctx->load) {
    966:     PetscInt  refine = 0;
    967:     PetscBool isHierarchy;
    968:     DM       *dms;
    969:     char      typeName[256];
    970:     PetscBool flg;

    972:     PetscCall(LoadFromFile(comm, ctx->load_filename, dm));
    973:     PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
    974:     PetscCall(PetscOptionsFList("-dm_mat_type", "Matrix type used for created matrices", "DMSetMatType", MatList, MATAIJ, typeName, sizeof(typeName), &flg));
    975:     if (flg) PetscCall(DMSetMatType(*dm, typeName));
    976:     PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0));
    977:     PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0));
    978:     PetscOptionsEnd();
    979:     if (refine) {
    980:       PetscCall(SetupDiscretization(*dm, ctx));
    981:       PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
    982:     }
    983:     PetscCall(PetscCalloc1(refine, &dms));
    984:     if (isHierarchy) PetscCall(DMRefineHierarchy(*dm, refine, dms));
    985:     for (PetscInt r = 0; r < refine; r++) {
    986:       Mat M;
    987:       DM  dmr = dms[r];
    988:       Vec u, ur;

    990:       if (!isHierarchy) {
    991:         PetscCall(DMRefine(*dm, PetscObjectComm((PetscObject)*dm), &dmr));
    992:         PetscCall(DMSetCoarseDM(dmr, *dm));
    993:       }
    994:       PetscCall(DMCreateInterpolation(*dm, dmr, &M, NULL));
    995:       PetscCall(DMGetNamedGlobalVector(*dm, "solution_", &u));
    996:       PetscCall(DMGetNamedGlobalVector(dmr, "solution_", &ur));
    997:       PetscCall(MatInterpolate(M, u, ur));
    998:       PetscCall(DMRestoreNamedGlobalVector(*dm, "solution_", &u));
    999:       PetscCall(DMRestoreNamedGlobalVector(dmr, "solution_", &ur));
    1000:       PetscCall(MatDestroy(&M));
    1001:       if (!isHierarchy) PetscCall(DMSetCoarseDM(dmr, NULL));
    1002:       PetscCall(DMDestroy(dm));
    1003:       *dm = dmr;
    1004:     }
    1005:     if (refine && !isHierarchy) PetscCall(DMSetRefineLevel(*dm, 0));
    1006:     PetscCall(PetscFree(dms));
    1007:   } else {
    1008:     PetscCall(DMCreate(comm, dm));
    1009:     PetscCall(DMSetType(*dm, DMPLEX));
    1010:     PetscCall(DMSetFromOptions(*dm));
    1011:     PetscCall(DMLocalizeCoordinates(*dm));
    1012:     {
    1013:       char      convType[256];
    1014:       PetscBool flg;
    1015:       PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX");
    1016:       PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg));
    1017:       PetscOptionsEnd();
    1018:       if (flg) {
    1019:         DM dmConv;
    1020:         PetscCall(DMConvert(*dm, convType, &dmConv));
    1021:         if (dmConv) {
    1022:           PetscCall(DMDestroy(dm));
    1023:           *dm = dmConv;
    1024:           PetscCall(DMSetFromOptions(*dm));
    1025:           PetscCall(DMSetUp(*dm));
    1026:         }
    1027:       }
    1028:     }
    1029:   }
    1030:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    1031:   PetscFunctionReturn(PETSC_SUCCESS);
    1032: }

    1034: /* Make potential field zero mean */
    1035: static PetscErrorCode ZeroMeanPotential(DM dm, Vec u)
    1036: {
    1037:   PetscScalar vals[NUM_FIELDS];
    1038:   PetscDS     ds;
    1039:   IS          is;

    1041:   PetscFunctionBeginUser;
    1042:   PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is));
    1043:   PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS");
    1044:   PetscCall(DMGetDS(dm, &ds));
    1045:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
    1046:   PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
    1047:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));
    1048:   PetscCall(VecISShift(u, is, -vals[P_FIELD_ID]));
    1049:   PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
    1050:   PetscFunctionReturn(PETSC_SUCCESS);
    1051: }

    1053: /* Compute initial conditions and exclude potential from local truncation error
    1054:    Since we are solving a DAE, once the initial conditions for the differential
    1055:    variables are set, we need to compute the corresponding value for the
    1056:    algebraic variables. We do so by creating a subDM for the potential only
    1057:    and solve a static problem with SNES */
    1058: static PetscErrorCode SetInitialConditionsAndTolerances(TS ts, PetscInt nv, Vec vecs[], PetscBool valid)
    1059: {
    1060:   DM         dm;
    1061:   Vec        tu, u, p, lsource, subaux, vatol, vrtol;
    1062:   PetscReal  t, atol, rtol;
    1063:   PetscInt   fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID};
    1064:   IS         isp;
    1065:   DM         dmp;
    1066:   VecScatter sctp = NULL;
    1067:   PetscDS    ds;
    1068:   SNES       snes;
    1069:   KSP        ksp;
    1070:   PC         pc;
    1071:   AppCtx    *ctx;

    1073:   PetscFunctionBeginUser;
    1074:   PetscCall(TSGetDM(ts, &dm));
    1075:   PetscCall(TSGetApplicationContext(ts, &ctx));
    1076:   if (valid) {
    1077:     PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, NULL));
    1078:     PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
    1079:     PetscCall(DMCreateGlobalVector(dm, &vatol));
    1080:     PetscCall(DMCreateGlobalVector(dm, &vrtol));
    1081:     PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
    1082:     PetscCall(VecSet(vatol, atol));
    1083:     PetscCall(VecISSet(vatol, isp, -1));
    1084:     PetscCall(VecSet(vrtol, rtol));
    1085:     PetscCall(VecISSet(vrtol, isp, -1));
    1086:     PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
    1087:     PetscCall(VecDestroy(&vatol));
    1088:     PetscCall(VecDestroy(&vrtol));
    1089:     PetscCall(ISDestroy(&isp));
    1090:     for (PetscInt i = 0; i < nv; i++) { PetscCall(ZeroMeanPotential(dm, vecs[i])); }
    1091:     PetscFunctionReturn(PETSC_SUCCESS);
    1092:   }
    1093:   PetscCall(DMCreateSubDM(dm, NUM_FIELDS - 1, fields + 1, &isp, &dmp));
    1094:   PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)isp));
    1095:   PetscCall(DMSetMatType(dmp, MATAIJ));
    1096:   PetscCall(DMGetDS(dmp, &ds));
    1097:   //PetscCall(PetscDSSetResidual(ds, 0, P_0, P_1_aux));
    1098:   PetscCall(PetscDSSetResidual(ds, 0, 0, P_1_aux));
    1099:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux));
    1100:   PetscCall(DMPlexSetSNESLocalFEM(dmp, PETSC_FALSE, NULL));

    1102:   PetscCall(DMCreateGlobalVector(dmp, &p));

    1104:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)dmp), &snes));
    1105:   PetscCall(SNESSetDM(snes, dmp));
    1106:   PetscCall(SNESSetOptionsPrefix(snes, "initial_"));
    1107:   PetscCall(SNESSetErrorIfNotConverged(snes, PETSC_TRUE));
    1108:   PetscCall(SNESGetKSP(snes, &ksp));
    1109:   PetscCall(KSPSetType(ksp, KSPFGMRES));
    1110:   PetscCall(KSPGetPC(ksp, &pc));
    1111:   PetscCall(PCSetType(pc, PCGAMG));
    1112:   PetscCall(SNESSetFromOptions(snes));
    1113:   PetscCall(SNESSetUp(snes));

    1115:   /* Loop over input vectors and compute corresponding potential */
    1116:   for (PetscInt i = 0; i < nv; i++) {
    1117:     PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);

    1119:     u = vecs[i];
    1120:     if (!valid) { /* Assumes entries in u are not valid */
    1121:       PetscCall(TSGetTime(ts, &t));
    1122:       switch (ctx->ic_num) {
    1123:       case 0:
    1124:         funcs[C_FIELD_ID] = initial_conditions_C_0;
    1125:         break;
    1126:       case 1:
    1127:         funcs[C_FIELD_ID] = initial_conditions_C_1;
    1128:         break;
    1129:       case 2:
    1130:         funcs[C_FIELD_ID] = initial_conditions_C_2;
    1131:         break;
    1132:       default:
    1133:         SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unknown IC");
    1134:       }
    1135:       funcs[P_FIELD_ID] = zerof;
    1136:       PetscCall(DMProjectFunction(dm, t, funcs, NULL, INSERT_ALL_VALUES, u));
    1137:     }

    1139:     /* pass conductivity and source information via auxiliary data */
    1140:     PetscCall(DMGetGlobalVector(dm, &tu));
    1141:     PetscCall(VecCopy(u, tu));
    1142:     PetscCall(VecISSet(tu, isp, 0.0));
    1143:     PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &lsource));
    1144:     PetscCall(DMCreateLocalVector(dm, &subaux));
    1145:     PetscCall(DMGlobalToLocal(dm, tu, INSERT_VALUES, subaux));
    1146:     PetscCall(DMRestoreGlobalVector(dm, &tu));
    1147:     PetscCall(VecAXPY(subaux, 1.0, lsource));
    1148:     PetscCall(VecViewFromOptions(subaux, NULL, "-initial_aux_view"));
    1149:     PetscCall(DMSetAuxiliaryVec(dmp, NULL, 0, 0, subaux));
    1150:     PetscCall(VecDestroy(&subaux));

    1152:     /* solve the subproblem */
    1153:     if (!sctp) PetscCall(VecScatterCreate(u, isp, p, NULL, &sctp));
    1154:     PetscCall(VecScatterBegin(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
    1155:     PetscCall(VecScatterEnd(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD));
    1156:     PetscCall(SNESSolve(snes, NULL, p));

    1158:     /* scatter from potential only to full space */
    1159:     PetscCall(VecScatterBegin(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
    1160:     PetscCall(VecScatterEnd(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE));
    1161:     PetscCall(ZeroMeanPotential(dm, u));
    1162:   }
    1163:   PetscCall(VecDestroy(&p));
    1164:   PetscCall(DMDestroy(&dmp));
    1165:   PetscCall(SNESDestroy(&snes));
    1166:   PetscCall(VecScatterDestroy(&sctp));

    1168:   /* exclude potential from computation of the LTE */
    1169:   PetscCall(DMCreateGlobalVector(dm, &vatol));
    1170:   PetscCall(DMCreateGlobalVector(dm, &vrtol));
    1171:   PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL));
    1172:   PetscCall(VecSet(vatol, atol));
    1173:   PetscCall(VecISSet(vatol, isp, -1));
    1174:   PetscCall(VecSet(vrtol, rtol));
    1175:   PetscCall(VecISSet(vrtol, isp, -1));
    1176:   PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol));
    1177:   PetscCall(VecDestroy(&vatol));
    1178:   PetscCall(VecDestroy(&vrtol));
    1179:   PetscCall(ISDestroy(&isp));
    1180:   PetscFunctionReturn(PETSC_SUCCESS);
    1181: }

    1183: /* mesh adaption context */
    1184: typedef struct {
    1185:   VecTagger refineTag;
    1186:   DMLabel   adaptLabel;
    1187:   PetscInt  cnt;
    1188: } AdaptCtx;

    1190: static PetscErrorCode ResizeSetUp(TS ts, PetscInt nstep, PetscReal time, Vec u, PetscBool *resize, void *vctx)
    1191: {
    1192:   AdaptCtx *ctx = (AdaptCtx *)vctx;
    1193:   Vec       ellVecCells, ellVecCellsF;
    1194:   DM        dm, plex;
    1195:   PetscDS   ds;
    1196:   PetscReal norm;
    1197:   PetscInt  cStart, cEnd;

    1199:   PetscFunctionBeginUser;
    1200:   PetscCall(TSGetDM(ts, &dm));
    1201:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1202:   PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
    1203:   PetscCall(DMDestroy(&plex));
    1204:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), NUM_FIELDS * (cEnd - cStart), PETSC_DECIDE, &ellVecCellsF));
    1205:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), cEnd - cStart, PETSC_DECIDE, &ellVecCells));
    1206:   PetscCall(VecSetBlockSize(ellVecCellsF, NUM_FIELDS));
    1207:   PetscCall(DMGetDS(dm, &ds));
    1208:   PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
    1209:   PetscCall(DMPlexComputeCellwiseIntegralFEM(dm, u, ellVecCellsF, NULL));
    1210:   PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));
    1211:   PetscCall(VecStrideGather(ellVecCellsF, C_FIELD_ID, ellVecCells, INSERT_VALUES));
    1212:   PetscCall(VecDestroy(&ellVecCellsF));
    1213:   PetscCall(VecNorm(ellVecCells, NORM_1, &norm));
    1214:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "STEP %d norm %g\n", (int)nstep, (double)norm));
    1215:   if (norm && !ctx->cnt) {
    1216:     IS refineIS;

    1218:     *resize = PETSC_TRUE;
    1219:     if (!ctx->refineTag) {
    1220:       VecTaggerBox refineBox;
    1221:       refineBox.min = PETSC_MACHINE_EPSILON;
    1222:       refineBox.max = PETSC_MAX_REAL;

    1224:       PetscCall(VecTaggerCreate(PETSC_COMM_SELF, &ctx->refineTag));
    1225:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->refineTag, "refine_"));
    1226:       PetscCall(VecTaggerSetType(ctx->refineTag, VECTAGGERABSOLUTE));
    1227:       PetscCall(VecTaggerAbsoluteSetBox(ctx->refineTag, &refineBox));
    1228:       PetscCall(VecTaggerSetFromOptions(ctx->refineTag));
    1229:       PetscCall(VecTaggerSetUp(ctx->refineTag));
    1230:       PetscCall(PetscObjectViewFromOptions((PetscObject)ctx->refineTag, NULL, "-tag_view"));
    1231:     }
    1232:     PetscCall(DMLabelDestroy(&ctx->adaptLabel));
    1233:     PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)ts), "adapt", &ctx->adaptLabel));
    1234:     PetscCall(VecTaggerComputeIS(ctx->refineTag, ellVecCells, &refineIS, NULL));
    1235:     PetscCall(DMLabelSetStratumIS(ctx->adaptLabel, DM_ADAPT_REFINE, refineIS));
    1236:     PetscCall(ISDestroy(&refineIS));
    1237: #if 0
    1238:     void (*funcs[NUM_FIELDS])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
    1239:     Vec ellVec;

    1241:     funcs[P_FIELD_ID] = ellipticity_fail;
    1242:     funcs[C_FIELD_ID] = NULL;

    1244:     PetscCall(DMGetGlobalVector(dm, &ellVec));
    1245:     PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
    1246:     PetscCall(VecViewFromOptions(ellVec,NULL,"-view_amr_ell"));
    1247:     PetscCall(DMRestoreGlobalVector(dm, &ellVec));
    1248: #endif
    1249:     ctx->cnt++;
    1250:   } else {
    1251:     ctx->cnt = 0;
    1252:   }
    1253:   PetscCall(VecDestroy(&ellVecCells));
    1254:   PetscFunctionReturn(PETSC_SUCCESS);
    1255: }

    1257: static PetscErrorCode ResizeTransfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *vctx)
    1258: {
    1259:   AdaptCtx *actx = (AdaptCtx *)vctx;
    1260:   AppCtx   *ctx;
    1261:   DM        dm, adm;
    1262:   PetscReal time;

    1264:   PetscFunctionBeginUser;
    1265:   PetscCall(TSGetDM(ts, &dm));
    1266:   PetscCheck(actx->adaptLabel, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptLabel");
    1267:   PetscCall(DMAdaptLabel(dm, actx->adaptLabel, &adm));
    1268:   PetscCheck(adm, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adapted DM");
    1269:   PetscCall(TSGetTime(ts, &time));
    1270:   PetscCall(DMLabelDestroy(&actx->adaptLabel));
    1271:   for (PetscInt i = 0; i < nv; i++) {
    1272:     PetscCall(DMCreateGlobalVector(adm, &vecsout[i]));
    1273:     PetscCall(DMForestTransferVec(dm, vecsin[i], adm, vecsout[i], PETSC_TRUE, time));
    1274:   }
    1275:   PetscCall(DMForestSetAdaptivityForest(adm, NULL));
    1276:   PetscCall(DMSetCoarseDM(adm, NULL));
    1277:   PetscCall(DMSetLocalSection(adm, NULL));
    1278:   PetscCall(TSSetDM(ts, adm));
    1279:   PetscCall(TSGetTime(ts, &time));
    1280:   PetscCall(TSGetApplicationContext(ts, &ctx));
    1281:   PetscCall(DMSetNullSpaceConstructor(adm, P_FIELD_ID, CreatePotentialNullSpace));
    1282:   PetscCall(ProjectSource(adm, time, ctx));
    1283:   PetscCall(SetInitialConditionsAndTolerances(ts, nv, vecsout, PETSC_TRUE));
    1284:   PetscCall(DMDestroy(&adm));
    1285:   PetscFunctionReturn(PETSC_SUCCESS);
    1286: }

    1288: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
    1289: {
    1290:   PetscFunctionBeginUser;
    1291:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
    1292:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
    1293:   PetscCall(PetscViewerFileSetName(*viewer, filename));
    1294:   PetscFunctionReturn(PETSC_SUCCESS);
    1295: }

    1297: /* Monitor relevant functionals */
    1298: static PetscErrorCode Monitor(TS ts, PetscInt stepnum, PetscReal time, Vec u, void *vctx)
    1299: {
    1300:   PetscScalar vals[2 * NUM_FIELDS];
    1301:   DM          dm;
    1302:   PetscDS     ds;
    1303:   AppCtx     *ctx;

    1305:   PetscFunctionBeginUser;
    1306:   PetscCall(TSGetDM(ts, &dm));
    1307:   PetscCall(TSGetApplicationContext(ts, &ctx));
    1308:   PetscCall(DMGetDS(dm, &ds));

    1310:   /* monitor energy and potential average */
    1311:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average));
    1312:   PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL));
    1313:   PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero));

    1315:   /* monitor ellipticity_fail */
    1316:   PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail));
    1317:   PetscCall(DMPlexComputeIntegralFEM(dm, u, vals + NUM_FIELDS, NULL));
    1318:   if (ctx->ellipticity) {
    1319:     void (*funcs[NUM_FIELDS])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
    1320:     Vec         ellVec;
    1321:     PetscViewer viewer;
    1322:     char        filename[PETSC_MAX_PATH_LEN];

    1324:     funcs[P_FIELD_ID] = ellipticity_fail;
    1325:     funcs[C_FIELD_ID] = NULL;

    1327:     PetscCall(DMGetGlobalVector(dm, &ellVec));
    1328:     PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec));
    1329:     PetscCall(PetscSNPrintf(filename, sizeof filename, "ellipticity_fail-%03" PetscInt_FMT ".vtu", stepnum));
    1330:     PetscCall(OutputVTK(dm, filename, &viewer));
    1331:     PetscCall(VecView(ellVec, viewer));
    1332:     PetscCall(PetscViewerDestroy(&viewer));
    1333:     PetscCall(DMRestoreGlobalVector(dm, &ellVec));
    1334:   }
    1335:   PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy));

    1337:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%4" PetscInt_FMT " TS: time %g, energy %g, intp %g, ell %g\n", stepnum, (double)time, (double)PetscRealPart(vals[C_FIELD_ID]), (double)PetscRealPart(vals[P_FIELD_ID]), (double)PetscRealPart(vals[NUM_FIELDS + C_FIELD_ID])));
    1338:   PetscFunctionReturn(PETSC_SUCCESS);
    1339: }

    1341: /* Save restart information */
    1342: static PetscErrorCode MonitorSave(TS ts, PetscInt steps, PetscReal time, Vec u, void *vctx)
    1343: {
    1344:   DM                dm;
    1345:   AppCtx           *ctx        = (AppCtx *)vctx;
    1346:   PetscInt          save_every = ctx->save_every;
    1347:   TSConvergedReason reason;

    1349:   PetscFunctionBeginUser;
    1350:   if (!ctx->save) PetscFunctionReturn(PETSC_SUCCESS);
    1351:   PetscCall(TSGetDM(ts, &dm));
    1352:   PetscCall(TSGetConvergedReason(ts, &reason));
    1353:   if ((save_every > 0 && steps % save_every == 0) || (save_every == -1 && reason) || save_every < -1) PetscCall(SaveToFile(dm, u, ctx->save_filename));
    1354:   PetscFunctionReturn(PETSC_SUCCESS);
    1355: }

    1357: /* Make potential zero mean after SNES solve */
    1358: static PetscErrorCode PostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
    1359: {
    1360:   DM       dm;
    1361:   Vec      u = Y[stageindex];
    1362:   SNES     snes;
    1363:   PetscInt nits, lits, stepnum;
    1364:   AppCtx  *ctx;

    1366:   PetscFunctionBeginUser;
    1367:   PetscCall(TSGetDM(ts, &dm));
    1368:   PetscCall(ZeroMeanPotential(dm, u));

    1370:   PetscCall(TSGetApplicationContext(ts, &ctx));
    1371:   if (ctx->test_restart) PetscFunctionReturn(PETSC_SUCCESS);

    1373:   /* monitor linear and nonlinear iterations */
    1374:   PetscCall(TSGetStepNumber(ts, &stepnum));
    1375:   PetscCall(TSGetSNES(ts, &snes));
    1376:   PetscCall(SNESGetIterationNumber(snes, &nits));
    1377:   PetscCall(SNESGetLinearSolveIterations(snes, &lits));

    1379:   /* if function evals in TSDIRK are zero in the first stage, it is FSAL */
    1380:   if (stageindex == 0) {
    1381:     PetscBool dirk;
    1382:     PetscInt  nf;

    1384:     PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
    1385:     PetscCall(SNESGetNumberFunctionEvals(snes, &nf));
    1386:     if (dirk && nf == 0) nits = 0;
    1387:   }
    1388:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "         step %" PetscInt_FMT " stage %" PetscInt_FMT " nonlinear its %" PetscInt_FMT ", linear its %" PetscInt_FMT "\n", stepnum, stageindex, nits, lits));
    1389:   PetscFunctionReturn(PETSC_SUCCESS);
    1390: }

    1392: static PetscErrorCode VecViewFlux(Vec u, const char *opts)
    1393: {
    1394:   Vec       fluxVec;
    1395:   DM        dmFlux, dm, plex;
    1396:   PetscInt  dim;
    1397:   PetscFE   feC, feFluxC, feNormC;
    1398:   PetscBool simplex, has;

    1400:   void (*funcs[])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) = {normc, flux};

    1402:   PetscFunctionBeginUser;
    1403:   PetscCall(PetscOptionsHasName(NULL, NULL, opts, &has));
    1404:   if (!has) PetscFunctionReturn(PETSC_SUCCESS);
    1405:   PetscCall(VecGetDM(u, &dm));
    1406:   PetscCall(DMGetDimension(dm, &dim));
    1407:   PetscCall(DMGetField(dm, C_FIELD_ID, NULL, (PetscObject *)&feC));
    1408:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1409:   PetscCall(DMPlexIsSimplex(plex, &simplex));
    1410:   PetscCall(DMDestroy(&plex));
    1411:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, simplex, "flux_", -1, &feFluxC));
    1412:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, "normc_", -1, &feNormC));
    1413:   PetscCall(PetscFECopyQuadrature(feC, feFluxC));
    1414:   PetscCall(PetscFECopyQuadrature(feC, feNormC));
    1415:   PetscCall(DMClone(dm, &dmFlux));
    1416:   PetscCall(DMSetNumFields(dmFlux, 1));
    1417:   PetscCall(DMSetField(dmFlux, 0, NULL, (PetscObject)feNormC));
    1418:   /* paraview segfaults! */
    1419:   //PetscCall(DMSetField(dmFlux, 1, NULL, (PetscObject)feFluxC));
    1420:   PetscCall(DMCreateDS(dmFlux));
    1421:   PetscCall(PetscFEDestroy(&feFluxC));
    1422:   PetscCall(PetscFEDestroy(&feNormC));

    1424:   PetscCall(DMGetGlobalVector(dmFlux, &fluxVec));
    1425:   PetscCall(DMProjectField(dmFlux, 0.0, u, funcs, INSERT_VALUES, fluxVec));
    1426:   PetscCall(VecViewFromOptions(fluxVec, NULL, opts));
    1427:   PetscCall(DMRestoreGlobalVector(dmFlux, &fluxVec));
    1428:   PetscCall(DMDestroy(&dmFlux));
    1429:   PetscFunctionReturn(PETSC_SUCCESS);
    1430: }

    1432: static PetscErrorCode Run(MPI_Comm comm, AppCtx *ctx)
    1433: {
    1434:   DM        dm;
    1435:   TS        ts;
    1436:   Vec       u;
    1437:   AdaptCtx *actx;
    1438:   PetscBool flg;

    1440:   PetscFunctionBeginUser;
    1441:   if (ctx->test_restart) {
    1442:     PetscViewer subviewer;
    1443:     PetscMPIInt rank;

    1445:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    1446:     PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
    1447:     if (ctx->load) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d loading from %s\n", rank, ctx->load_filename));
    1448:     if (ctx->save) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d saving to %s\n", rank, ctx->save_filename));
    1449:     PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
    1450:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    1451:   } else {
    1452:     PetscCall(PetscPrintf(comm, "----------------------------\n"));
    1453:     PetscCall(PetscPrintf(comm, "Simulation parameters:\n"));
    1454:     PetscCall(PetscPrintf(comm, "  r    : %g\n", (double)ctx->r));
    1455:     PetscCall(PetscPrintf(comm, "  eps  : %g\n", (double)ctx->eps));
    1456:     PetscCall(PetscPrintf(comm, "  alpha: %g\n", (double)ctx->alpha));
    1457:     PetscCall(PetscPrintf(comm, "  gamma: %g\n", (double)ctx->gamma));
    1458:     PetscCall(PetscPrintf(comm, "  D    : %g\n", (double)ctx->D));
    1459:     PetscCall(PetscPrintf(comm, "  c    : %g\n", (double)ctx->c));
    1460:     if (ctx->load) PetscCall(PetscPrintf(comm, "  load : %s\n", ctx->load_filename));
    1461:     else PetscCall(PetscPrintf(comm, "  IC   : %" PetscInt_FMT "\n", ctx->ic_num));
    1462:     PetscCall(PetscPrintf(comm, "  S    : %" PetscInt_FMT "\n", ctx->source_num));
    1463:     PetscCall(PetscPrintf(comm, "  x0   : (%g,%g)\n", (double)ctx->x0[0], (double)ctx->x0[1]));
    1464:     PetscCall(PetscPrintf(comm, "----------------------------\n"));
    1465:   }

    1467:   if (!ctx->test_restart) PetscCall(PetscLogStagePush(SetupStage));
    1468:   PetscCall(CreateMesh(comm, &dm, ctx));
    1469:   PetscCall(SetupDiscretization(dm, ctx));

    1471:   PetscCall(TSCreate(comm, &ts));
    1472:   PetscCall(TSSetApplicationContext(ts, ctx));

    1474:   PetscCall(TSSetDM(ts, dm));
    1475:   if (ctx->test_restart) {
    1476:     PetscViewer subviewer;
    1477:     PetscMPIInt rank;
    1478:     PetscInt    level;

    1480:     PetscCall(DMGetRefineLevel(dm, &level));
    1481:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    1482:     PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
    1483:     PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d DM refinement level %" PetscInt_FMT "\n", rank, level));
    1484:     PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer));
    1485:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    1486:   }

    1488:   if (ctx->test_restart) PetscCall(TSSetMaxSteps(ts, 1));
    1489:   PetscCall(TSSetMaxTime(ts, 10.0));
    1490:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1491:   if (!ctx->test_restart) PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
    1492:   PetscCall(TSMonitorSet(ts, MonitorSave, ctx, NULL));
    1493:   PetscCall(PetscNew(&actx));
    1494:   if (ctx->amr) PetscCall(TSSetResize(ts, PETSC_TRUE, ResizeSetUp, ResizeTransfer, actx));
    1495:   PetscCall(TSSetPostStage(ts, PostStage));
    1496:   PetscCall(TSSetMaxSNESFailures(ts, -1));
    1497:   PetscCall(TSSetFromOptions(ts));

    1499:   PetscCall(DMCreateGlobalVector(dm, &u));
    1500:   PetscCall(PetscObjectSetName((PetscObject)u, "solution_"));
    1501:   PetscCall(DMHasNamedGlobalVector(dm, "solution_", &flg));
    1502:   if (flg) {
    1503:     Vec ru;

    1505:     PetscCall(DMGetNamedGlobalVector(dm, "solution_", &ru));
    1506:     PetscCall(VecCopy(ru, u));
    1507:     PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &ru));
    1508:   }
    1509:   PetscCall(SetInitialConditionsAndTolerances(ts, 1, &u, PETSC_FALSE));
    1510:   PetscCall(TSSetSolution(ts, u));
    1511:   PetscCall(VecDestroy(&u));
    1512:   PetscCall(DMDestroy(&dm));
    1513:   if (!ctx->test_restart) PetscCall(PetscLogStagePop());

    1515:   if (!ctx->test_restart) PetscCall(PetscLogStagePush(SolveStage));
    1516:   PetscCall(TSSolve(ts, NULL));
    1517:   if (!ctx->test_restart) PetscCall(PetscLogStagePop());

    1519:   PetscCall(TSGetSolution(ts, &u));
    1520:   PetscCall(VecViewFromOptions(u, NULL, "-final_view"));
    1521:   PetscCall(VecViewFlux(u, "-final_flux_view"));

    1523:   PetscCall(TSDestroy(&ts));
    1524:   PetscCall(VecTaggerDestroy(&actx->refineTag));
    1525:   PetscCall(PetscFree(actx));
    1526:   PetscFunctionReturn(PETSC_SUCCESS);
    1527: }

    1529: int main(int argc, char **argv)
    1530: {
    1531:   AppCtx ctx;

    1533:   PetscFunctionBeginUser;
    1534:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    1535:   PetscCall(ProcessOptions(&ctx));
    1536:   PetscCall(PetscLogStageRegister("Setup", &SetupStage));
    1537:   PetscCall(PetscLogStageRegister("Solve", &SolveStage));
    1538:   if (ctx.test_restart) { /* Test sequences of save and loads */
    1539:     PetscMPIInt rank;

    1541:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

    1543:     /* test saving */
    1544:     ctx.load = PETSC_FALSE;
    1545:     ctx.save = PETSC_TRUE;
    1546:     /* sequential save */
    1547:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential save\n"));
    1548:     PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_seq_%d.h5", rank));
    1549:     PetscCall(Run(PETSC_COMM_SELF, &ctx));
    1550:     /* parallel save */
    1551:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel save\n"));
    1552:     PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_par.h5"));
    1553:     PetscCall(Run(PETSC_COMM_WORLD, &ctx));

    1555:     /* test loading */
    1556:     ctx.load = PETSC_TRUE;
    1557:     ctx.save = PETSC_FALSE;
    1558:     /* sequential and parallel runs from sequential save */
    1559:     PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_seq_0.h5"));
    1560:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from sequential save\n"));
    1561:     PetscCall(Run(PETSC_COMM_SELF, &ctx));
    1562:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from sequential save\n"));
    1563:     PetscCall(Run(PETSC_COMM_WORLD, &ctx));
    1564:     /* sequential and parallel runs from parallel save */
    1565:     PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_par.h5"));
    1566:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from parallel save\n"));
    1567:     PetscCall(Run(PETSC_COMM_SELF, &ctx));
    1568:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from parallel save\n"));
    1569:     PetscCall(Run(PETSC_COMM_WORLD, &ctx));
    1570:   } else { /* Run the simulation */
    1571:     PetscCall(Run(PETSC_COMM_WORLD, &ctx));
    1572:   }
    1573:   PetscCall(PetscFinalize());
    1574:   return 0;
    1575: }

    1577: /*TEST

    1579:   testset:
    1580:     args: -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 1 -initial_snes_test_jacobian -snes_test_jacobian -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none

    1582:     test:
    1583:       suffix: 0
    1584:       nsize: {{1 2}}
    1585:       args: -dm_refine 1 -lump {{0 1}}

    1587:     test:
    1588:       suffix: 0_dirk
    1589:       nsize: {{1 2}}
    1590:       args: -dm_refine 1 -ts_type dirk

    1592:     test:
    1593:       suffix: 0_dirk_mg
    1594:       nsize: {{1 2}}
    1595:       args: -dm_refine_hierarchy 1 -ts_type dirk -pc_type mg -mg_levels_pc_type bjacobi -mg_levels_sub_pc_factor_levels 2 -mg_levels_sub_pc_factor_mat_ordering_type rcm -mg_levels_sub_pc_factor_reuse_ordering -mg_coarse_pc_type svd -lump {{0 1}}

    1597:     test:
    1598:       suffix: 0_dirk_fieldsplit
    1599:       nsize: {{1 2}}
    1600:       args: -dm_refine 1 -ts_type dirk -pc_type fieldsplit -pc_fieldsplit_type multiplicative -lump {{0 1}}

    1602:     test:
    1603:       requires: p4est
    1604:       suffix: 0_p4est
    1605:       nsize: {{1 2}}
    1606:       args: -dm_refine 1 -dm_plex_convert_type p4est -lump 0

    1608:     test:
    1609:       suffix: 0_periodic
    1610:       nsize: {{1 2}}
    1611:       args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -lump {{0 1}}

    1613:     test:
    1614:       requires: p4est
    1615:       suffix: 0_p4est_periodic
    1616:       nsize: {{1 2}}
    1617:       args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -dm_plex_convert_type p4est -lump 0

    1619:     test:
    1620:       requires: p4est
    1621:       suffix: 0_p4est_mg
    1622:       nsize: {{1 2}}
    1623:       args: -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_plex_convert_type p4est -pc_type mg -mg_coarse_pc_type svd -mg_levels_pc_type svd -lump 0

    1625:   testset:
    1626:     requires: hdf5
    1627:     args: -test_restart -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type mg -mg_levels_pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -petscpartitioner_type simple -test_restart

    1629:     test:
    1630:       requires: !single
    1631:       suffix: restart
    1632:       nsize: {{1 2}separate output}
    1633:       args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 0

    1635:     test:
    1636:       requires: triangle
    1637:       suffix: restart_simplex
    1638:       nsize: {{1 2}separate output}
    1639:       args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 1

    1641:     test:
    1642:       requires: !single
    1643:       suffix: restart_refonly
    1644:       nsize: {{1 2}separate output}
    1645:       args: -dm_refine 1 -dm_plex_simplex 0

    1647:     test:
    1648:       requires: triangle
    1649:       suffix: restart_simplex_refonly
    1650:       nsize: {{1 2}separate output}
    1651:       args: -dm_refine 1 -dm_plex_simplex 1

    1653: TEST*/