Actual source code: ex18.c

      1: static char help[] = "Hybrid Finite Element-Finite Volume Example.\n";
      2: /*F
      3:   Here we are advecting a passive tracer in a harmonic velocity field, defined by
      4: a forcing function $f$:
      5: \begin{align}
      6:   -\Delta \mathbf{u} + f &= 0 \\
      7:   \frac{\partial\phi}{\partial t} + \nabla\cdot \phi \mathbf{u} &= 0
      8: \end{align}
      9: F*/

     11: #include <petscdmplex.h>
     12: #include <petscds.h>
     13: #include <petscts.h>

     15: #include <petsc/private/dmpleximpl.h>

     17: typedef enum {
     18:   VEL_ZERO,
     19:   VEL_CONSTANT,
     20:   VEL_HARMONIC,
     21:   VEL_SHEAR
     22: } VelocityDistribution;

     24: typedef enum {
     25:   ZERO,
     26:   CONSTANT,
     27:   GAUSSIAN,
     28:   TILTED,
     29:   DELTA
     30: } PorosityDistribution;

     32: static PetscErrorCode constant_u_2d(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);

     34: /*
     35:   FunctionalFn - Calculates the value of a functional of the solution at a point

     37:   Input Parameters:
     38: + dm   - The DM
     39: . time - The TS time
     40: . x    - The coordinates of the evaluation point
     41: . u    - The field values at point x
     42: - ctx  - A user context, or NULL

     44:   Output Parameter:
     45: . f    - The value of the functional at point x

     47: */
     48: typedef PetscErrorCode (*FunctionalFn)(DM, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);

     50: typedef struct _n_Functional *Functional;
     51: struct _n_Functional {
     52:   char        *name;
     53:   FunctionalFn func;
     54:   void        *ctx;
     55:   PetscInt     offset;
     56:   Functional   next;
     57: };

     59: typedef struct {
     60:   /* Problem definition */
     61:   PetscBool useFV; /* Use a finite volume scheme for advection */
     62:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
     63:   VelocityDistribution velocityDist;
     64:   PorosityDistribution porosityDist;
     65:   PetscReal            inflowState;
     66:   PetscReal            source[3];
     67:   /* Monitoring */
     68:   PetscInt    numMonitorFuncs, maxMonitorFunc;
     69:   Functional *monitorFuncs;
     70:   PetscInt    errorFunctional;
     71:   Functional  functionalRegistry;
     72: } AppCtx;

     74: static AppCtx *globalUser;

     76: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
     77: {
     78:   const char *velocityDist[4] = {"zero", "constant", "harmonic", "shear"};
     79:   const char *porosityDist[5] = {"zero", "constant", "gaussian", "tilted", "delta"};
     80:   PetscInt    vd, pd, d;
     81:   PetscBool   flg;

     83:   PetscFunctionBeginUser;
     84:   options->useFV           = PETSC_FALSE;
     85:   options->velocityDist    = VEL_HARMONIC;
     86:   options->porosityDist    = ZERO;
     87:   options->inflowState     = -2.0;
     88:   options->numMonitorFuncs = 0;
     89:   options->source[0]       = 0.5;
     90:   options->source[1]       = 0.5;
     91:   options->source[2]       = 0.5;

     93:   PetscOptionsBegin(comm, "", "Magma Dynamics Options", "DMPLEX");
     94:   PetscCall(PetscOptionsBool("-use_fv", "Use the finite volume method for advection", "ex18.c", options->useFV, &options->useFV, NULL));
     95:   vd = options->velocityDist;
     96:   PetscCall(PetscOptionsEList("-velocity_dist", "Velocity distribution type", "ex18.c", velocityDist, 4, velocityDist[options->velocityDist], &vd, NULL));
     97:   options->velocityDist = (VelocityDistribution)vd;
     98:   pd                    = options->porosityDist;
     99:   PetscCall(PetscOptionsEList("-porosity_dist", "Initial porosity distribution type", "ex18.c", porosityDist, 5, porosityDist[options->porosityDist], &pd, NULL));
    100:   options->porosityDist = (PorosityDistribution)pd;
    101:   PetscCall(PetscOptionsReal("-inflow_state", "The inflow state", "ex18.c", options->inflowState, &options->inflowState, NULL));
    102:   d = 2;
    103:   PetscCall(PetscOptionsRealArray("-source_loc", "The source location", "ex18.c", options->source, &d, &flg));
    104:   PetscCheck(!flg || d == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must give dim coordinates for the source location, not %" PetscInt_FMT, d);
    105:   PetscOptionsEnd();
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: static PetscErrorCode ProcessMonitorOptions(MPI_Comm comm, AppCtx *options)
    110: {
    111:   Functional func;
    112:   char      *names[256];
    113:   PetscInt   f;

    115:   PetscFunctionBeginUser;
    116:   PetscOptionsBegin(comm, "", "Simulation Monitor Options", "DMPLEX");
    117:   options->numMonitorFuncs = PETSC_STATIC_ARRAY_LENGTH(names);
    118:   PetscCall(PetscOptionsStringArray("-monitor", "List of functionals to monitor", "", names, &options->numMonitorFuncs, NULL));
    119:   PetscCall(PetscMalloc1(options->numMonitorFuncs, &options->monitorFuncs));
    120:   for (f = 0; f < options->numMonitorFuncs; ++f) {
    121:     for (func = options->functionalRegistry; func; func = func->next) {
    122:       PetscBool match;

    124:       PetscCall(PetscStrcasecmp(names[f], func->name, &match));
    125:       if (match) break;
    126:     }
    127:     PetscCheck(func, comm, PETSC_ERR_USER, "No known functional '%s'", names[f]);
    128:     options->monitorFuncs[f] = func;
    129:     /* Jed inserts a de-duplication of functionals here */
    130:     PetscCall(PetscFree(names[f]));
    131:   }
    132:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
    133:   options->maxMonitorFunc = -1;
    134:   for (func = options->functionalRegistry; func; func = func->next) {
    135:     for (f = 0; f < options->numMonitorFuncs; ++f) {
    136:       Functional call = options->monitorFuncs[f];

    138:       if (func->func == call->func && func->ctx == call->ctx) options->maxMonitorFunc = PetscMax(options->maxMonitorFunc, func->offset);
    139:     }
    140:   }
    141:   PetscOptionsEnd();
    142:   PetscFunctionReturn(PETSC_SUCCESS);
    143: }

    145: static PetscErrorCode FunctionalRegister(Functional *functionalRegistry, const char name[], PetscInt *offset, FunctionalFn func, void *ctx)
    146: {
    147:   Functional *ptr, f;
    148:   PetscInt    lastoffset = -1;

    150:   PetscFunctionBeginUser;
    151:   for (ptr = functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
    152:   PetscCall(PetscNew(&f));
    153:   PetscCall(PetscStrallocpy(name, &f->name));
    154:   f->offset = lastoffset + 1;
    155:   f->func   = func;
    156:   f->ctx    = ctx;
    157:   f->next   = NULL;
    158:   *ptr      = f;
    159:   *offset   = f->offset;
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: static PetscErrorCode FunctionalDestroy(Functional *link)
    164: {
    165:   Functional next, l;

    167:   PetscFunctionBeginUser;
    168:   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
    169:   l     = *link;
    170:   *link = NULL;
    171:   for (; l; l = next) {
    172:     next = l->next;
    173:     PetscCall(PetscFree(l->name));
    174:     PetscCall(PetscFree(l));
    175:   }
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: static void f0_zero_u(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[])
    180: {
    181:   PetscInt comp;
    182:   for (comp = 0; comp < dim; ++comp) f0[comp] = u[comp];
    183: }

    185: static void f0_constant_u(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[])
    186: {
    187:   PetscScalar wind[3] = {0.0, 0.0, 0.0};
    188:   PetscInt    comp;

    190:   PetscCallAbort(PETSC_COMM_SELF, constant_u_2d(dim, t, x, Nf, wind, NULL));
    191:   for (comp = 0; comp < dim && comp < 3; ++comp) f0[comp] = u[comp] - wind[comp];
    192: }

    194: static void f1_constant_u(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[])
    195: {
    196:   PetscInt comp;
    197:   for (comp = 0; comp < dim * dim; ++comp) f1[comp] = 0.0;
    198: }

    200: static void g0_constant_uu(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 g0[])
    201: {
    202:   PetscInt d;
    203:   for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
    204: }

    206: static void g0_constant_pp(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 g0[])
    207: {
    208:   g0[0] = 1.0;
    209: }

    211: static void f0_lap_u(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[])
    212: {
    213:   PetscInt comp;
    214:   for (comp = 0; comp < dim; ++comp) f0[comp] = 4.0;
    215: }

    217: static void f1_lap_u(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[])
    218: {
    219:   PetscInt comp, d;
    220:   for (comp = 0; comp < dim; ++comp) {
    221:     for (d = 0; d < dim; ++d) f1[comp * dim + d] = u_x[comp * dim + d];
    222:   }
    223: }

    225: static void f0_lap_periodic_u(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[])
    226: {
    227:   f0[0] = -PetscSinReal(2.0 * PETSC_PI * x[0]);
    228:   f0[1] = 2.0 * PETSC_PI * x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]);
    229: }

    231: static void f0_lap_doubly_periodic_u(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[])
    232: {
    233:   f0[0] = -2.0 * PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]);
    234:   f0[1] = 2.0 * PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]);
    235: }

    237: void g3_uu(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 g3[])
    238: {
    239:   const PetscInt Ncomp = dim;
    240:   PetscInt       compI, d;

    242:   for (compI = 0; compI < Ncomp; ++compI) {
    243:     for (d = 0; d < dim; ++d) g3[((compI * Ncomp + compI) * dim + d) * dim + d] = 1.0;
    244:   }
    245: }

    247: /* \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \mathbf{u} + \phi \nabla \cdot \mathbf{u} = 0 */
    248: static void f0_advection(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[])
    249: {
    250:   PetscInt d;
    251:   f0[0] = u_t[dim];
    252:   for (d = 0; d < dim; ++d) f0[0] += u[dim] * u_x[d * dim + d] + u_x[dim * dim + d] * u[d];
    253: }

    255: static void f1_advection(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[])
    256: {
    257:   PetscInt d;
    258:   for (d = 0; d < dim; ++d) f1[0] = 0.0;
    259: }

    261: void g0_adv_pp(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 g0[])
    262: {
    263:   PetscInt d;
    264:   g0[0] = u_tShift;
    265:   for (d = 0; d < dim; ++d) g0[0] += u_x[d * dim + d];
    266: }

    268: void g1_adv_pp(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 g1[])
    269: {
    270:   PetscInt d;
    271:   for (d = 0; d < dim; ++d) g1[d] = u[d];
    272: }

    274: void g0_adv_pu(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 g0[])
    275: {
    276:   PetscInt d;
    277:   for (d = 0; d < dim; ++d) g0[0] += u_x[dim * dim + d];
    278: }

    280: void g1_adv_pu(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 g1[])
    281: {
    282:   PetscInt d;
    283:   for (d = 0; d < dim; ++d) g1[d * dim + d] = u[dim];
    284: }

    286: static void riemann_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx)
    287: {
    288:   PetscReal wind[3] = {0.0, 1.0, 0.0};
    289:   PetscReal wn      = DMPlex_DotRealD_Internal(PetscMin(dim, 3), wind, n);

    291:   flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
    292: }

    294: static void riemann_coupled_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx)
    295: {
    296:   PetscReal wn = DMPlex_DotD_Internal(dim, uL, n);

    298: #if 1
    299:   flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
    300: #else
    301:   /* if (fabs(uL[0] - wind[0]) > 1.0e-7 || fabs(uL[1] - wind[1]) > 1.0e-7) PetscPrintf(PETSC_COMM_SELF, "wind (%g, %g) uL (%g, %g) uR (%g, %g)\n", wind[0], wind[1], uL[0], uL[1], uR[0], uR[1]); */
    302:   /* Smear it out */
    303:   flux[0] = 0.5 * ((uL[dim] + uR[dim]) + (uL[dim] - uR[dim]) * tanh(1.0e5 * wn)) * wn;
    304: #endif
    305: }

    307: static PetscErrorCode zero_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    308: {
    309:   u[0] = 0.0;
    310:   u[1] = 0.0;
    311:   return PETSC_SUCCESS;
    312: }

    314: static PetscErrorCode constant_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    315: {
    316:   u[0] = 0.0;
    317:   u[1] = 1.0;
    318:   return PETSC_SUCCESS;
    319: }

    321: /* Coordinates of the point which was at x at t = 0 */
    322: static PetscErrorCode constant_x_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    323: {
    324:   const PetscReal t = *((PetscReal *)ctx);
    325:   u[0]              = x[0];
    326:   u[1]              = x[1] + t;
    327: #if 0
    328:   PetscCall(DMLocalizeCoordinate(globalUser->dm, u, PETSC_FALSE, u));
    329: #else
    330:   u[1] = u[1] - (int)PetscRealPart(u[1]);
    331: #endif
    332:   return PETSC_SUCCESS;
    333: }

    335: /*
    336:   In 2D we use the exact solution:

    338:     u   = x^2 + y^2
    339:     v   = 2 x^2 - 2xy
    340:     phi = h(x + y + (u + v) t)
    341:     f_x = f_y = 4

    343:   so that

    345:     -\Delta u + f = <-4, -4> + <4, 4> = 0
    346:     {\partial\phi}{\partial t} - \nabla\cdot \phi u = 0
    347:     h_t(x + y + (u + v) t) - u . grad phi - phi div u
    348:   = u h' + v h'              - u h_x - v h_y
    349:   = 0

    351: We will conserve phi since

    353:     \nabla \cdot u = 2x - 2x = 0

    355: Also try h((x + ut)^2 + (y + vt)^2), so that

    357:     h_t((x + ut)^2 + (y + vt)^2) - u . grad phi - phi div u
    358:   = 2 h' (u (x + ut) + v (y + vt)) - u h_x - v h_y
    359:   = 2 h' (u (x + ut) + v (y + vt)) - u h' 2 (x + u t) - v h' 2 (y + vt)
    360:   = 2 h' (u (x + ut) + v (y + vt)  - u (x + u t) - v (y + vt))
    361:   = 0

    363: */
    364: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    365: {
    366:   u[0] = x[0] * x[0] + x[1] * x[1];
    367:   u[1] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[1];
    368:   return PETSC_SUCCESS;
    369: }

    371: /*
    372:   In 2D we use the exact, periodic solution:

    374:     u   =  sin(2 pi x)/4 pi^2
    375:     v   = -y cos(2 pi x)/2 pi
    376:     phi = h(x + y + (u + v) t)
    377:     f_x = -sin(2 pi x)
    378:     f_y = 2 pi y cos(2 pi x)

    380:   so that

    382:     -\Delta u + f = <sin(2pi x),  -2pi y cos(2pi x)> + <-sin(2pi x), 2pi y cos(2pi x)> = 0

    384: We will conserve phi since

    386:     \nabla \cdot u = cos(2pi x)/2pi - cos(2pi x)/2pi = 0
    387: */
    388: static PetscErrorCode periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    389: {
    390:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
    391:   u[1] = -x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]) / (2.0 * PETSC_PI);
    392:   return PETSC_SUCCESS;
    393: }

    395: /*
    396:   In 2D we use the exact, doubly periodic solution:

    398:     u   =  sin(2 pi x) cos(2 pi y)/4 pi^2
    399:     v   = -sin(2 pi y) cos(2 pi x)/4 pi^2
    400:     phi = h(x + y + (u + v) t)
    401:     f_x = -2sin(2 pi x) cos(2 pi y)
    402:     f_y =  2sin(2 pi y) cos(2 pi x)

    404:   so that

    406:     -\Delta u + f = <2 sin(2pi x) cos(2pi y),  -2 sin(2pi y) cos(2pi x)> + <-2 sin(2pi x) cos(2pi y), 2 sin(2pi y) cos(2pi x)> = 0

    408: We will conserve phi since

    410:     \nabla \cdot u = cos(2pi x) cos(2pi y)/2pi - cos(2pi y) cos(2pi x)/2pi = 0
    411: */
    412: static PetscErrorCode doubly_periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    413: {
    414:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]) / PetscSqr(2.0 * PETSC_PI);
    415:   u[1] = -PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
    416:   return PETSC_SUCCESS;
    417: }

    419: static PetscErrorCode shear_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    420: {
    421:   u[0] = x[1] - 0.5;
    422:   u[1] = 0.0;
    423:   return PETSC_SUCCESS;
    424: }

    426: static PetscErrorCode initialVelocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    427: {
    428:   PetscInt d;
    429:   for (d = 0; d < dim; ++d) u[d] = 0.0;
    430:   return PETSC_SUCCESS;
    431: }

    433: static PetscErrorCode zero_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    434: {
    435:   u[0] = 0.0;
    436:   return PETSC_SUCCESS;
    437: }

    439: static PetscErrorCode constant_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    440: {
    441:   u[0] = 1.0;
    442:   return PETSC_SUCCESS;
    443: }

    445: static PetscErrorCode delta_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    446: {
    447:   PetscReal   x0[2];
    448:   PetscScalar xn[2];

    450:   x0[0] = globalUser->source[0];
    451:   x0[1] = globalUser->source[1];
    452:   PetscCall(constant_x_2d(dim, time, x0, Nf, xn, ctx));
    453:   {
    454:     const PetscReal xi  = x[0] - PetscRealPart(xn[0]);
    455:     const PetscReal eta = x[1] - PetscRealPart(xn[1]);
    456:     const PetscReal r2  = xi * xi + eta * eta;

    458:     u[0] = r2 < 1.0e-7 ? 1.0 : 0.0;
    459:   }
    460:   return PETSC_SUCCESS;
    461: }

    463: /*
    464:   Gaussian blob, initially centered on (0.5, 0.5)

    466:   xi = x(t) - x0, eta = y(t) - y0

    468: where x(t), y(t) are the integral curves of v(t),

    470:   dx/dt . grad f = v . f

    472: Check: constant v(t) = {v0, w0}, x(t) = {x0 + v0 t, y0 + w0 t}

    474:   v0 f_x + w0 f_y = v . f
    475: */
    476: static PetscErrorCode gaussian_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    477: {
    478:   const PetscReal x0[2] = {0.5, 0.5};
    479:   const PetscReal sigma = 1.0 / 6.0;
    480:   PetscScalar     xn[2];

    482:   PetscCall(constant_x_2d(dim, time, x0, Nf, xn, ctx));
    483:   {
    484:     /* const PetscReal xi  = x[0] + (sin(2.0*PETSC_PI*x[0])/(4.0*PETSC_PI*PETSC_PI))*t - x0[0]; */
    485:     /* const PetscReal eta = x[1] + (-x[1]*cos(2.0*PETSC_PI*x[0])/(2.0*PETSC_PI))*t - x0[1]; */
    486:     const PetscReal xi  = x[0] - PetscRealPart(xn[0]);
    487:     const PetscReal eta = x[1] - PetscRealPart(xn[1]);
    488:     const PetscReal r2  = xi * xi + eta * eta;

    490:     u[0] = PetscExpReal(-r2 / (2.0 * sigma * sigma)) / (sigma * PetscSqrtReal(2.0 * PETSC_PI));
    491:   }
    492:   return PETSC_SUCCESS;
    493: }

    495: static PetscErrorCode tilted_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    496: {
    497:   PetscReal       x0[3];
    498:   const PetscReal wind[3] = {0.0, 1.0, 0.0};
    499:   const PetscReal t       = *((PetscReal *)ctx);

    501:   DMPlex_WaxpyD_Internal(2, -t, wind, x, x0);
    502:   if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
    503:   else u[0] = -2.0; /* Inflow state */
    504:   return PETSC_SUCCESS;
    505: }

    507: static PetscErrorCode tilted_phi_coupled_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    508: {
    509:   PetscReal       ur[3];
    510:   PetscReal       x0[3];
    511:   const PetscReal t = *((PetscReal *)ctx);

    513:   ur[0] = PetscRealPart(u[0]);
    514:   ur[1] = PetscRealPart(u[1]);
    515:   ur[2] = PetscRealPart(u[2]);
    516:   DMPlex_WaxpyD_Internal(2, -t, ur, x, x0);
    517:   if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
    518:   else u[0] = -2.0; /* Inflow state */
    519:   return PETSC_SUCCESS;
    520: }

    522: static PetscErrorCode advect_inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    523: {
    524:   AppCtx *user = (AppCtx *)ctx;

    526:   PetscFunctionBeginUser;
    527:   xG[0] = user->inflowState;
    528:   PetscFunctionReturn(PETSC_SUCCESS);
    529: }

    531: static PetscErrorCode advect_outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    532: {
    533:   PetscFunctionBeginUser;
    534:   //xG[0] = xI[dim];
    535:   xG[0] = xI[2];
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

    539: static PetscErrorCode ExactSolution(DM dm, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
    540: {
    541:   AppCtx  *user = (AppCtx *)ctx;
    542:   PetscInt dim;

    544:   PetscFunctionBeginUser;
    545:   PetscCall(DMGetDimension(dm, &dim));
    546:   switch (user->porosityDist) {
    547:   case TILTED:
    548:     if (user->velocityDist == VEL_ZERO) PetscCall(tilted_phi_2d(dim, time, x, 2, u, (void *)&time));
    549:     else PetscCall(tilted_phi_coupled_2d(dim, time, x, 2, u, (void *)&time));
    550:     break;
    551:   case GAUSSIAN:
    552:     PetscCall(gaussian_phi_2d(dim, time, x, 2, u, (void *)&time));
    553:     break;
    554:   case DELTA:
    555:     PetscCall(delta_phi_2d(dim, time, x, 2, u, (void *)&time));
    556:     break;
    557:   default:
    558:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
    559:   }
    560:   PetscFunctionReturn(PETSC_SUCCESS);
    561: }

    563: static PetscErrorCode Functional_Error(DM dm, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
    564: {
    565:   AppCtx     *user      = (AppCtx *)ctx;
    566:   PetscScalar yexact[3] = {0, 0, 0};

    568:   PetscFunctionBeginUser;
    569:   PetscCall(ExactSolution(dm, time, x, yexact, ctx));
    570:   f[user->errorFunctional] = PetscAbsScalar(y[0] - yexact[0]);
    571:   PetscFunctionReturn(PETSC_SUCCESS);
    572: }

    574: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    575: {
    576:   PetscFunctionBeginUser;
    577:   PetscCall(DMCreate(comm, dm));
    578:   PetscCall(DMSetType(*dm, DMPLEX));
    579:   PetscCall(DMSetFromOptions(*dm));
    580:   PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
    581:   PetscFunctionReturn(PETSC_SUCCESS);
    582: }

    584: static PetscErrorCode SetupBC(DM dm, AppCtx *user)
    585: {
    586:   PetscErrorCode (*exactFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    587:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    588:   PetscDS        prob;
    589:   DMLabel        label;
    590:   PetscBool      check;
    591:   PetscInt       dim, n = 3;
    592:   const char    *prefix;

    594:   PetscFunctionBeginUser;
    595:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    596:   PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL));
    597:   PetscCall(DMGetDimension(dm, &dim));
    598:   /* Set initial guesses and exact solutions */
    599:   switch (dim) {
    600:   case 2:
    601:     user->initialGuess[0] = initialVelocity;
    602:     switch (user->porosityDist) {
    603:     case ZERO:
    604:       user->initialGuess[1] = zero_phi;
    605:       break;
    606:     case CONSTANT:
    607:       user->initialGuess[1] = constant_phi;
    608:       break;
    609:     case GAUSSIAN:
    610:       user->initialGuess[1] = gaussian_phi_2d;
    611:       break;
    612:     case DELTA:
    613:       user->initialGuess[1] = delta_phi_2d;
    614:       break;
    615:     case TILTED:
    616:       if (user->velocityDist == VEL_ZERO) user->initialGuess[1] = tilted_phi_2d;
    617:       else user->initialGuess[1] = tilted_phi_coupled_2d;
    618:       break;
    619:     }
    620:     break;
    621:   default:
    622:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
    623:   }
    624:   exactFuncs[0] = user->initialGuess[0];
    625:   exactFuncs[1] = user->initialGuess[1];
    626:   switch (dim) {
    627:   case 2:
    628:     switch (user->velocityDist) {
    629:     case VEL_ZERO:
    630:       exactFuncs[0] = zero_u_2d;
    631:       break;
    632:     case VEL_CONSTANT:
    633:       exactFuncs[0] = constant_u_2d;
    634:       break;
    635:     case VEL_HARMONIC:
    636:       switch (bdt[0]) {
    637:       case DM_BOUNDARY_PERIODIC:
    638:         switch (bdt[1]) {
    639:         case DM_BOUNDARY_PERIODIC:
    640:           exactFuncs[0] = doubly_periodic_u_2d;
    641:           break;
    642:         default:
    643:           exactFuncs[0] = periodic_u_2d;
    644:           break;
    645:         }
    646:         break;
    647:       default:
    648:         exactFuncs[0] = quadratic_u_2d;
    649:         break;
    650:       }
    651:       break;
    652:     case VEL_SHEAR:
    653:       exactFuncs[0] = shear_bc;
    654:       break;
    655:     default:
    656:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %" PetscInt_FMT, dim);
    657:     }
    658:     break;
    659:   default:
    660:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
    661:   }
    662:   {
    663:     PetscBool isImplicit = PETSC_FALSE;

    665:     PetscCall(PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit));
    666:     if (user->velocityDist == VEL_CONSTANT && !isImplicit) user->initialGuess[0] = exactFuncs[0];
    667:   }
    668:   PetscCall(PetscOptionsHasName(NULL, NULL, "-dmts_check", &check));
    669:   if (check) {
    670:     user->initialGuess[0] = exactFuncs[0];
    671:     user->initialGuess[1] = exactFuncs[1];
    672:   }
    673:   /* Set BC */
    674:   PetscCall(DMGetDS(dm, &prob));
    675:   PetscCall(DMGetLabel(dm, "marker", &label));
    676:   PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], user));
    677:   PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], user));
    678:   if (label) {
    679:     const PetscInt id = 1;

    681:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, user, NULL));
    682:   }
    683:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
    684:   if (label && user->useFV) {
    685:     const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};

    687:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 1, 0, NULL, (PetscVoidFn *)advect_inflow, NULL, user, NULL));
    688:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 1, 0, NULL, (PetscVoidFn *)advect_outflow, NULL, user, NULL));
    689:   }
    690:   PetscFunctionReturn(PETSC_SUCCESS);
    691: }

    693: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
    694: {
    695:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    696:   PetscDS        prob;
    697:   PetscInt       n = 3;
    698:   const char    *prefix;

    700:   PetscFunctionBeginUser;
    701:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    702:   PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL));
    703:   PetscCall(DMGetDS(dm, &prob));
    704:   switch (user->velocityDist) {
    705:   case VEL_ZERO:
    706:     PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_constant_u));
    707:     break;
    708:   case VEL_CONSTANT:
    709:     PetscCall(PetscDSSetResidual(prob, 0, f0_constant_u, f1_constant_u));
    710:     PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_constant_uu, NULL, NULL, NULL));
    711:     PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_constant_pp, NULL, NULL, NULL));
    712:     break;
    713:   case VEL_HARMONIC:
    714:     switch (bdt[0]) {
    715:     case DM_BOUNDARY_PERIODIC:
    716:       switch (bdt[1]) {
    717:       case DM_BOUNDARY_PERIODIC:
    718:         PetscCall(PetscDSSetResidual(prob, 0, f0_lap_doubly_periodic_u, f1_lap_u));
    719:         break;
    720:       default:
    721:         PetscCall(PetscDSSetResidual(prob, 0, f0_lap_periodic_u, f1_lap_u));
    722:         break;
    723:       }
    724:       break;
    725:     default:
    726:       PetscCall(PetscDSSetResidual(prob, 0, f0_lap_u, f1_lap_u));
    727:       break;
    728:     }
    729:     PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
    730:     break;
    731:   case VEL_SHEAR:
    732:     PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_lap_u));
    733:     PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
    734:     break;
    735:   }
    736:   PetscCall(PetscDSSetResidual(prob, 1, f0_advection, f1_advection));
    737:   PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_adv_pp, g1_adv_pp, NULL, NULL));
    738:   PetscCall(PetscDSSetJacobian(prob, 1, 0, g0_adv_pu, g1_adv_pu, NULL, NULL));
    739:   if (user->velocityDist == VEL_ZERO) PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_advection));
    740:   else PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_coupled_advection));

    742:   PetscCall(FunctionalRegister(&user->functionalRegistry, "Error", &user->errorFunctional, Functional_Error, user));
    743:   PetscFunctionReturn(PETSC_SUCCESS);
    744: }

    746: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
    747: {
    748:   DM              cdm = dm;
    749:   PetscQuadrature q;
    750:   PetscFE         fe[2];
    751:   PetscFV         fv;
    752:   MPI_Comm        comm;
    753:   PetscInt        dim;

    755:   PetscFunctionBeginUser;
    756:   /* Create finite element */
    757:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    758:   PetscCall(DMGetDimension(dm, &dim));
    759:   PetscCall(PetscFECreateDefault(comm, dim, dim, PETSC_FALSE, "velocity_", PETSC_DEFAULT, &fe[0]));
    760:   PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
    761:   PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_FALSE, "porosity_", PETSC_DEFAULT, &fe[1]));
    762:   PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
    763:   PetscCall(PetscObjectSetName((PetscObject)fe[1], "porosity"));

    765:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv));
    766:   PetscCall(PetscObjectSetName((PetscObject)fv, "porosity"));
    767:   PetscCall(PetscFVSetFromOptions(fv));
    768:   PetscCall(PetscFVSetNumComponents(fv, 1));
    769:   PetscCall(PetscFVSetSpatialDimension(fv, dim));
    770:   PetscCall(PetscFEGetQuadrature(fe[0], &q));
    771:   PetscCall(PetscFVSetQuadrature(fv, q));

    773:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
    774:   if (user->useFV) PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fv));
    775:   else PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
    776:   PetscCall(DMCreateDS(dm));
    777:   PetscCall(SetupProblem(dm, user));

    779:   /* Set discretization and boundary conditions for each mesh */
    780:   while (cdm) {
    781:     PetscCall(DMCopyDisc(dm, cdm));
    782:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    783:     /* Coordinates were never localized for coarse meshes */
    784:     if (cdm) PetscCall(DMLocalizeCoordinates(cdm));
    785:   }
    786:   PetscCall(PetscFEDestroy(&fe[0]));
    787:   PetscCall(PetscFEDestroy(&fe[1]));
    788:   PetscCall(PetscFVDestroy(&fv));
    789:   PetscFunctionReturn(PETSC_SUCCESS);
    790: }

    792: static PetscErrorCode CreateDM(MPI_Comm comm, AppCtx *user, DM *dm)
    793: {
    794:   PetscFunctionBeginUser;
    795:   PetscCall(CreateMesh(comm, user, dm));
    796:   /* Handle refinement, etc. */
    797:   PetscCall(DMSetFromOptions(*dm));
    798:   /* Construct ghost cells */
    799:   if (user->useFV) {
    800:     DM gdm;

    802:     PetscCall(DMPlexConstructGhostCells(*dm, NULL, NULL, &gdm));
    803:     PetscCall(DMDestroy(dm));
    804:     *dm = gdm;
    805:   }
    806:   /* Localize coordinates */
    807:   PetscCall(DMLocalizeCoordinates(*dm));
    808:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
    809:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    810:   /* Setup problem */
    811:   PetscCall(SetupDiscretization(*dm, user));
    812:   /* Setup BC */
    813:   PetscCall(SetupBC(*dm, user));
    814:   PetscFunctionReturn(PETSC_SUCCESS);
    815: }

    817: static PetscErrorCode SetInitialConditionFVM(DM dm, Vec X, PetscInt field, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx)
    818: {
    819:   PetscDS            prob;
    820:   DM                 dmCell;
    821:   Vec                cellgeom;
    822:   const PetscScalar *cgeom;
    823:   PetscScalar       *x;
    824:   PetscInt           dim, Nf, cStart, cEnd, c;

    826:   PetscFunctionBeginUser;
    827:   PetscCall(DMGetDS(dm, &prob));
    828:   PetscCall(DMGetDimension(dm, &dim));
    829:   PetscCall(PetscDSGetNumFields(prob, &Nf));
    830:   PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL));
    831:   PetscCall(VecGetDM(cellgeom, &dmCell));
    832:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    833:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    834:   PetscCall(VecGetArray(X, &x));
    835:   for (c = cStart; c < cEnd; ++c) {
    836:     PetscFVCellGeom *cg;
    837:     PetscScalar     *xc;

    839:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
    840:     PetscCall(DMPlexPointGlobalFieldRef(dm, c, field, x, &xc));
    841:     if (xc) PetscCall((*func)(dim, 0.0, cg->centroid, Nf, xc, ctx));
    842:   }
    843:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    844:   PetscCall(VecRestoreArray(X, &x));
    845:   PetscFunctionReturn(PETSC_SUCCESS);
    846: }

    848: static PetscErrorCode MonitorFunctionals(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
    849: {
    850:   AppCtx            *user   = (AppCtx *)ctx;
    851:   char              *ftable = NULL;
    852:   DM                 dm;
    853:   PetscSection       s;
    854:   Vec                cellgeom;
    855:   const PetscScalar *x;
    856:   PetscScalar       *a;
    857:   PetscReal         *xnorms;
    858:   PetscInt           pStart, pEnd, p, Nf, f;

    860:   PetscFunctionBeginUser;
    861:   PetscCall(VecViewFromOptions(X, (PetscObject)ts, "-view_solution"));
    862:   PetscCall(VecGetDM(X, &dm));
    863:   PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL));
    864:   PetscCall(DMGetLocalSection(dm, &s));
    865:   PetscCall(PetscSectionGetNumFields(s, &Nf));
    866:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    867:   PetscCall(PetscCalloc1(Nf * 2, &xnorms));
    868:   PetscCall(VecGetArrayRead(X, &x));
    869:   for (p = pStart; p < pEnd; ++p) {
    870:     for (f = 0; f < Nf; ++f) {
    871:       PetscInt dof, cdof, d;

    873:       PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
    874:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
    875:       PetscCall(DMPlexPointGlobalFieldRead(dm, p, f, x, &a));
    876:       /* TODO Use constrained indices here */
    877:       for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 0] = PetscMax(xnorms[f * 2 + 0], PetscAbsScalar(a[d]));
    878:       for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 1] += PetscAbsScalar(a[d]);
    879:     }
    880:   }
    881:   PetscCall(VecRestoreArrayRead(X, &x));
    882:   if (stepnum >= 0) { /* No summary for final time */
    883:     DM                 dmCell, *fdm;
    884:     Vec               *fv;
    885:     const PetscScalar *cgeom;
    886:     PetscScalar      **fx;
    887:     PetscReal         *fmin, *fmax, *fint, *ftmp, t;
    888:     PetscInt           cStart, cEnd, c, fcount, f, num;

    890:     size_t ftableused, ftablealloc;

    892:     /* Functionals have indices after registering, this is an upper bound */
    893:     fcount = user->numMonitorFuncs;
    894:     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fint, fcount, &ftmp));
    895:     PetscCall(PetscMalloc3(fcount, &fdm, fcount, &fv, fcount, &fx));
    896:     for (f = 0; f < fcount; ++f) {
    897:       PetscSection fs;
    898:       const char  *name = user->monitorFuncs[f]->name;

    900:       fmin[f] = PETSC_MAX_REAL;
    901:       fmax[f] = PETSC_MIN_REAL;
    902:       fint[f] = 0;
    903:       /* Make monitor vecs */
    904:       PetscCall(DMClone(dm, &fdm[f]));
    905:       PetscCall(DMGetOutputSequenceNumber(dm, &num, &t));
    906:       PetscCall(DMSetOutputSequenceNumber(fdm[f], num, t));
    907:       PetscCall(PetscSectionClone(s, &fs));
    908:       PetscCall(PetscSectionSetFieldName(fs, 0, NULL));
    909:       PetscCall(PetscSectionSetFieldName(fs, 1, name));
    910:       PetscCall(DMSetLocalSection(fdm[f], fs));
    911:       PetscCall(PetscSectionDestroy(&fs));
    912:       PetscCall(DMGetGlobalVector(fdm[f], &fv[f]));
    913:       PetscCall(PetscObjectSetName((PetscObject)fv[f], name));
    914:       PetscCall(VecGetArray(fv[f], &fx[f]));
    915:     }
    916:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    917:     PetscCall(VecGetDM(cellgeom, &dmCell));
    918:     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    919:     PetscCall(VecGetArrayRead(X, &x));
    920:     for (c = cStart; c < cEnd; ++c) {
    921:       PetscFVCellGeom *cg;
    922:       PetscScalar     *cx;

    924:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
    925:       PetscCall(DMPlexPointGlobalFieldRead(dm, c, 1, x, &cx));
    926:       if (!cx) continue; /* not a global cell */
    927:       for (f = 0; f < user->numMonitorFuncs; ++f) {
    928:         Functional   func = user->monitorFuncs[f];
    929:         PetscScalar *fxc;

    931:         PetscCall(DMPlexPointGlobalFieldRef(dm, c, 1, fx[f], &fxc));
    932:         /* I need to make it easier to get interpolated values here */
    933:         PetscCall((*func->func)(dm, time, cg->centroid, cx, ftmp, func->ctx));
    934:         fxc[0] = ftmp[user->monitorFuncs[f]->offset];
    935:       }
    936:       for (f = 0; f < fcount; ++f) {
    937:         fmin[f] = PetscMin(fmin[f], ftmp[f]);
    938:         fmax[f] = PetscMax(fmax[f], ftmp[f]);
    939:         fint[f] += cg->volume * ftmp[f];
    940:       }
    941:     }
    942:     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    943:     PetscCall(VecRestoreArrayRead(X, &x));
    944:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
    945:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
    946:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fint, (PetscMPIInt)fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
    947:     /* Output functional data */
    948:     ftablealloc = fcount * 100;
    949:     ftableused  = 0;
    950:     PetscCall(PetscCalloc1(ftablealloc, &ftable));
    951:     for (f = 0; f < user->numMonitorFuncs; ++f) {
    952:       Functional func      = user->monitorFuncs[f];
    953:       PetscInt   id        = func->offset;
    954:       char       newline[] = "\n";
    955:       char       buffer[256], *p, *prefix;
    956:       size_t     countused, len;

    958:       /* Create string with functional outputs */
    959:       if (f % 3) {
    960:         PetscCall(PetscArraycpy(buffer, "  ", 2));
    961:         p = buffer + 2;
    962:       } else if (f) {
    963:         PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
    964:         p = buffer + sizeof(newline) - 1;
    965:       } else {
    966:         p = buffer;
    967:       }
    968:       PetscCall(PetscSNPrintfCount(p, sizeof buffer - (p - buffer), "%12s [%12.6g,%12.6g] int %12.6g", &countused, func->name, (double)fmin[id], (double)fmax[id], (double)fint[id]));
    969:       countused += p - buffer;
    970:       /* reallocate */
    971:       if (countused > ftablealloc - ftableused - 1) {
    972:         char *ftablenew;

    974:         ftablealloc = 2 * ftablealloc + countused;
    975:         PetscCall(PetscMalloc1(ftablealloc, &ftablenew));
    976:         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
    977:         PetscCall(PetscFree(ftable));
    978:         ftable = ftablenew;
    979:       }
    980:       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
    981:       ftableused += countused;
    982:       ftable[ftableused] = 0;
    983:       /* Output vecs */
    984:       PetscCall(VecRestoreArray(fv[f], &fx[f]));
    985:       PetscCall(PetscStrlen(func->name, &len));
    986:       PetscCall(PetscMalloc1(len + 2, &prefix));
    987:       PetscCall(PetscStrncpy(prefix, func->name, len + 2));
    988:       PetscCall(PetscStrlcat(prefix, "_", len + 2));
    989:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fv[f], prefix));
    990:       PetscCall(VecViewFromOptions(fv[f], NULL, "-vec_view"));
    991:       PetscCall(PetscFree(prefix));
    992:       PetscCall(DMRestoreGlobalVector(fdm[f], &fv[f]));
    993:       PetscCall(DMDestroy(&fdm[f]));
    994:     }
    995:     PetscCall(PetscFree4(fmin, fmax, fint, ftmp));
    996:     PetscCall(PetscFree3(fdm, fv, fx));
    997:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| (", stepnum, (double)time));
    998:     for (f = 0; f < Nf; ++f) {
    999:       if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
    1000:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 0]));
    1001:     }
    1002:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ") |x|_1 ("));
    1003:     for (f = 0; f < Nf; ++f) {
    1004:       if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
    1005:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 1]));
    1006:     }
    1007:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ")  %s\n", ftable ? ftable : ""));
    1008:     PetscCall(PetscFree(ftable));
    1009:   }
    1010:   PetscCall(PetscFree(xnorms));
    1011:   PetscFunctionReturn(PETSC_SUCCESS);
    1012: }

    1014: int main(int argc, char **argv)
    1015: {
    1016:   MPI_Comm  comm;
    1017:   TS        ts;
    1018:   DM        dm;
    1019:   Vec       u;
    1020:   AppCtx    user;
    1021:   PetscReal t0, t = 0.0;
    1022:   void     *ctxs[2] = {&t, &t};

    1024:   PetscFunctionBeginUser;
    1025:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    1026:   comm                    = PETSC_COMM_WORLD;
    1027:   user.functionalRegistry = NULL;
    1028:   globalUser              = &user;
    1029:   PetscCall(ProcessOptions(comm, &user));
    1030:   PetscCall(TSCreate(comm, &ts));
    1031:   PetscCall(TSSetType(ts, TSBEULER));
    1032:   PetscCall(CreateDM(comm, &user, &dm));
    1033:   PetscCall(TSSetDM(ts, dm));
    1034:   PetscCall(ProcessMonitorOptions(comm, &user));

    1036:   PetscCall(DMCreateGlobalVector(dm, &u));
    1037:   PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
    1038:   if (user.useFV) {
    1039:     PetscBool isImplicit = PETSC_FALSE;

    1041:     PetscCall(PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit));
    1042:     if (isImplicit) {
    1043:       PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
    1044:       PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
    1045:     }
    1046:     PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
    1047:     PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &user));
    1048:   } else {
    1049:     PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
    1050:     PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
    1051:     PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
    1052:   }
    1053:   if (user.useFV) PetscCall(TSMonitorSet(ts, MonitorFunctionals, &user, NULL));
    1054:   PetscCall(TSSetMaxSteps(ts, 1));
    1055:   PetscCall(TSSetMaxTime(ts, 2.0));
    1056:   PetscCall(TSSetTimeStep(ts, 0.01));
    1057:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1058:   PetscCall(TSSetFromOptions(ts));

    1060:   PetscCall(DMProjectFunction(dm, 0.0, user.initialGuess, ctxs, INSERT_VALUES, u));
    1061:   if (user.useFV) PetscCall(SetInitialConditionFVM(dm, u, 1, user.initialGuess[1], ctxs[1]));
    1062:   PetscCall(VecViewFromOptions(u, NULL, "-init_vec_view"));
    1063:   PetscCall(TSGetTime(ts, &t));
    1064:   t0 = t;
    1065:   PetscCall(DMTSCheckFromOptions(ts, u));
    1066:   PetscCall(TSSolve(ts, u));
    1067:   PetscCall(TSGetTime(ts, &t));
    1068:   if (t > t0) PetscCall(DMTSCheckFromOptions(ts, u));
    1069:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
    1070:   {
    1071:     PetscReal         ftime;
    1072:     PetscInt          nsteps;
    1073:     TSConvergedReason reason;

    1075:     PetscCall(TSGetSolveTime(ts, &ftime));
    1076:     PetscCall(TSGetStepNumber(ts, &nsteps));
    1077:     PetscCall(TSGetConvergedReason(ts, &reason));
    1078:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
    1079:   }

    1081:   PetscCall(VecDestroy(&u));
    1082:   PetscCall(DMDestroy(&dm));
    1083:   PetscCall(TSDestroy(&ts));
    1084:   PetscCall(PetscFree(user.monitorFuncs));
    1085:   PetscCall(FunctionalDestroy(&user.functionalRegistry));
    1086:   PetscCall(PetscFinalize());
    1087:   return 0;
    1088: }

    1090: /*TEST

    1092:   testset:
    1093:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3

    1095:     # 2D harmonic velocity, no porosity
    1096:     test:
    1097:       suffix: p1p1
    1098:       requires: !complex !single
    1099:       args: -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type ilu -pc_factor_shift_type nonzero -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1100:     test:
    1101:       suffix: p1p1_xper
    1102:       requires: !complex !single
    1103:       args: -dm_refine 1 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1104:     test:
    1105:       suffix: p1p1_xper_ref
    1106:       requires: !complex !single
    1107:       args: -dm_refine 2 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1108:     test:
    1109:       suffix: p1p1_xyper
    1110:       requires: !complex !single
    1111:       args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1112:     test:
    1113:       suffix: p1p1_xyper_ref
    1114:       requires: !complex !single
    1115:       args: -dm_refine 2 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1116:     test:
    1117:       suffix: p2p1
    1118:       requires: !complex !single
    1119:       args: -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
    1120:     test:
    1121:       suffix: p2p1_xyper
    1122:       requires: !complex !single
    1123:       args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check

    1125:     test:
    1126:       suffix: adv_1
    1127:       requires: !complex !single
    1128:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -ts_view -dm_view

    1130:     test:
    1131:       suffix: adv_2
    1132:       requires: !complex
    1133:       TODO: broken memory corruption
    1134:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3,4 -bc_outflow 1,2 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason

    1136:     test:
    1137:       suffix: adv_3
    1138:       requires: !complex
    1139:       TODO: broken memory corruption
    1140:       args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason

    1142:     test:
    1143:       suffix: adv_3_ex
    1144:       requires: !complex
    1145:       args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.1 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -ksp_max_it 100 -ts_view -dm_view

    1147:     test:
    1148:       suffix: adv_4
    1149:       requires: !complex
    1150:       TODO: broken memory corruption
    1151:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason

    1153:     # 2D Advection, box, delta
    1154:     test:
    1155:       suffix: adv_delta_yper_0
    1156:       requires: !complex
    1157:       TODO: broken
    1158:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error

    1160:     test:
    1161:       suffix: adv_delta_yper_1
    1162:       requires: !complex
    1163:       TODO: broken
    1164:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 1 -source_loc 0.416666,0.416666

    1166:     test:
    1167:       suffix: adv_delta_yper_2
    1168:       requires: !complex
    1169:       TODO: broken
    1170:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 2 -source_loc 0.458333,0.458333

    1172:     test:
    1173:       suffix: adv_delta_yper_fim_0
    1174:       requires: !complex
    1175:       TODO: broken
    1176:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

    1178:     test:
    1179:       suffix: adv_delta_yper_fim_1
    1180:       requires: !complex
    1181:       TODO: broken
    1182:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic

    1184:     test:
    1185:       suffix: adv_delta_yper_fim_2
    1186:       requires: !complex
    1187:       TODO: broken
    1188:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic

    1190:     test:
    1191:       suffix: adv_delta_yper_im_0
    1192:       requires: !complex
    1193:       TODO: broken
    1194:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

    1196:     test:
    1197:       suffix: adv_delta_yper_im_1
    1198:       requires: !complex
    1199:       TODO: broken
    1200:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666

    1202:     test:
    1203:       suffix: adv_delta_yper_im_2
    1204:       requires: !complex
    1205:       TODO: broken
    1206:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333

    1208:     test:
    1209:       suffix: adv_delta_yper_im_3
    1210:       requires: !complex
    1211:       TODO: broken
    1212:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

    1214:     #    I believe the nullspace is sin(pi y)
    1215:     test:
    1216:       suffix: adv_delta_yper_im_4
    1217:       requires: !complex
    1218:       TODO: broken
    1219:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666

    1221:     test:
    1222:       suffix: adv_delta_yper_im_5
    1223:       requires: !complex
    1224:       TODO: broken
    1225:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333

    1227:     test:
    1228:       suffix: adv_delta_yper_im_6
    1229:       requires: !complex
    1230:       TODO: broken
    1231:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type svd -snes_converged_reason
    1232:     # 2D Advection, magma benchmark 1

    1234:     test:
    1235:       suffix: adv_delta_shear_im_0
    1236:       requires: !complex
    1237:       TODO: broken
    1238:       args: -dm_plex_box_bd periodic,none -dm_refine 2 -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist shear -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 1,3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -source_loc 0.458333,0.708333
    1239:     # 2D Advection, box, gaussian

    1241:     test:
    1242:       suffix: adv_gauss
    1243:       requires: !complex
    1244:       TODO: broken
    1245:       args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view

    1247:     test:
    1248:       suffix: adv_gauss_im
    1249:       requires: !complex
    1250:       TODO: broken
    1251:       args: -use_fv -use_implicit -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

    1253:     test:
    1254:       suffix: adv_gauss_im_1
    1255:       requires: !complex
    1256:       TODO: broken
    1257:       args: -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

    1259:     test:
    1260:       suffix: adv_gauss_im_2
    1261:       requires: !complex
    1262:       TODO: broken
    1263:       args: -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

    1265:     test:
    1266:       suffix: adv_gauss_corner
    1267:       requires: !complex
    1268:       TODO: broken
    1269:       args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 2 -ts_view -dm_view

    1271:     # 2D Advection+Harmonic 12-
    1272:     test:
    1273:       suffix: adv_harm_0
    1274:       requires: !complex
    1275:       TODO: broken memory corruption
    1276:       args: -velocity_petscspace_degree 2 -use_fv -velocity_dist harmonic -porosity_dist gaussian -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -use_implicit -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason -snes_monitor -dmts_check

    1278:   #   Must check that FV BCs propagate to coarse meshes
    1279:   #   Must check that FV BC ids propagate to coarse meshes
    1280:   #   Must check that FE+FV BCs work at the same time
    1281:   # 2D Advection, matching wind in ex11 8-11
    1282:   #   NOTE implicit solves are limited by accuracy of FD Jacobian
    1283:   test:
    1284:     suffix: adv_0
    1285:     requires: !complex !single exodusii
    1286:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view

    1288:   test:
    1289:     suffix: adv_0_im
    1290:     requires: !complex exodusii
    1291:     TODO: broken  memory corruption
    1292:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu

    1294:   test:
    1295:     suffix: adv_0_im_2
    1296:     requires: !complex exodusii
    1297:     TODO: broken
    1298:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -snes_rtol 1.0e-7

    1300:   test:
    1301:     suffix: adv_0_im_3
    1302:     requires: !complex exodusii
    1303:     TODO: broken
    1304:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7

    1306:   test:
    1307:     suffix: adv_0_im_4
    1308:     requires: !complex exodusii
    1309:     TODO: broken
    1310:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7
    1311:   # 2D Advection, misc

    1313: TEST*/