Actual source code: ex11_sa.c

      1: static char help[] = "Second Order TVD Finite Volume Example.\n";
      2: /*F

      4: We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
      5: over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
      6: \begin{equation}
      7:   f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
      8: \end{equation}
      9: and then update the cell values given the cell volume.
     10: \begin{eqnarray}
     11:     f^L_i &-=& \frac{f_i}{vol^L} \\
     12:     f^R_i &+=& \frac{f_i}{vol^R}
     13: \end{eqnarray}

     15: As an example, we can consider the shallow water wave equation,
     16: \begin{eqnarray}
     17:      h_t + \nabla\cdot \left( uh                              \right) &=& 0 \\
     18:   (uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
     19: \end{eqnarray}
     20: where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.

     22: A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
     23: \begin{eqnarray}
     24:   f^{L,R}_h    &=& uh^{L,R} \cdot \hat n \\
     25:   f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
     26:   c^{L,R}      &=& \sqrt{g h^{L,R}} \\
     27:   s            &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\
     28:   f_i          &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
     29: \end{eqnarray}
     30: where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.

     32: The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
     33: over a neighborhood of the given element.

     35: The mesh is read in from an ExodusII file, usually generated by Cubit.
     36: F*/
     37: #include <petscts.h>
     38: #include <petscfv.h>
     39: #include <petscdmplex.h>
     40: #include <petscsf.h>
     41: #include <petscblaslapack.h>

     43: #define DIM 2 /* Geometric dimension */

     45: static PetscFunctionList PhysicsList;

     47: /* Represents continuum physical equations. */
     48: typedef struct _n_Physics *Physics;

     50: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
     51:  * discretization-independent, but its members depend on the scenario being solved. */
     52: typedef struct _n_Model *Model;

     54: /* 'User' implements a discretization of a continuous model. */
     55: typedef struct _n_User *User;

     57: typedef PetscErrorCode (*RiemannFunction)(const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscScalar *, void *);
     58: typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
     59: typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
     60: typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
     61: static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
     62: static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
     63: static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);

     65: struct FieldDescription {
     66:   const char *name;
     67:   PetscInt    dof;
     68: };

     70: typedef struct _n_FunctionalLink *FunctionalLink;
     71: struct _n_FunctionalLink {
     72:   char              *name;
     73:   FunctionalFunction func;
     74:   void              *ctx;
     75:   PetscInt           offset;
     76:   FunctionalLink     next;
     77: };

     79: struct _n_Physics {
     80:   RiemannFunction                riemann;
     81:   PetscInt                       dof;      /* number of degrees of freedom per cell */
     82:   PetscReal                      maxspeed; /* kludge to pick initial time step, need to add monitoring and step control */
     83:   void                          *data;
     84:   PetscInt                       nfields;
     85:   const struct FieldDescription *field_desc;
     86: };

     88: struct _n_Model {
     89:   MPI_Comm         comm; /* Does not do collective communication, but some error conditions can be collective */
     90:   Physics          physics;
     91:   FunctionalLink   functionalRegistry;
     92:   PetscInt         maxComputed;
     93:   PetscInt         numMonitored;
     94:   FunctionalLink  *functionalMonitored;
     95:   PetscInt         numCall;
     96:   FunctionalLink  *functionalCall;
     97:   SolutionFunction solution;
     98:   void            *solutionctx;
     99:   PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
    100: };

    102: struct _n_User {
    103:   PetscInt numSplitFaces;
    104:   PetscInt vtkInterval; /* For monitor */
    105:   Model    model;
    106: };

    108: static inline PetscScalar DotDIM(const PetscScalar *x, const PetscScalar *y)
    109: {
    110:   PetscInt    i;
    111:   PetscScalar prod = 0.0;

    113:   for (i = 0; i < DIM; i++) prod += x[i] * y[i];
    114:   return prod;
    115: }
    116: static inline PetscReal NormDIM(const PetscScalar *x)
    117: {
    118:   return PetscSqrtReal(PetscAbsScalar(DotDIM(x, x)));
    119: }
    120: static inline void axDIM(const PetscScalar a, PetscScalar *x)
    121: {
    122:   PetscInt i;
    123:   for (i = 0; i < DIM; i++) x[i] *= a;
    124: }
    125: static inline void waxDIM(const PetscScalar a, const PetscScalar *x, PetscScalar *w)
    126: {
    127:   PetscInt i;
    128:   for (i = 0; i < DIM; i++) w[i] = x[i] * a;
    129: }
    130: static inline void NormalSplitDIM(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
    131: { /* Split x into normal and tangential components */
    132:   PetscInt    i;
    133:   PetscScalar c;
    134:   c = DotDIM(x, n) / DotDIM(n, n);
    135:   for (i = 0; i < DIM; i++) {
    136:     xn[i] = c * n[i];
    137:     xt[i] = x[i] - xn[i];
    138:   }
    139: }

    141: static inline PetscScalar Dot2(const PetscScalar *x, const PetscScalar *y)
    142: {
    143:   return x[0] * y[0] + x[1] * y[1];
    144: }
    145: static inline PetscReal Norm2(const PetscScalar *x)
    146: {
    147:   return PetscSqrtReal(PetscAbsScalar(Dot2(x, x)));
    148: }
    149: static inline void Normalize2(PetscScalar *x)
    150: {
    151:   PetscReal a = 1. / Norm2(x);
    152:   x[0] *= a;
    153:   x[1] *= a;
    154: }
    155: static inline void Waxpy2(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
    156: {
    157:   w[0] = a * x[0] + y[0];
    158:   w[1] = a * x[1] + y[1];
    159: }
    160: static inline void Scale2(PetscScalar a, const PetscScalar *x, PetscScalar *y)
    161: {
    162:   y[0] = a * x[0];
    163:   y[1] = a * x[1];
    164: }

    166: static inline void WaxpyD(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
    167: {
    168:   PetscInt d;
    169:   for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
    170: }
    171: static inline PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
    172: {
    173:   PetscScalar sum = 0.0;
    174:   PetscInt    d;
    175:   for (d = 0; d < dim; ++d) sum += x[d] * y[d];
    176:   return sum;
    177: }
    178: static inline PetscReal NormD(PetscInt dim, const PetscScalar *x)
    179: {
    180:   return PetscSqrtReal(PetscAbsScalar(DotD(dim, x, x)));
    181: }

    183: static inline void NormalSplit(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
    184: { /* Split x into normal and tangential components */
    185:   Scale2(Dot2(x, n) / Dot2(n, n), n, xn);
    186:   Waxpy2(-1, xn, x, xt);
    187: }

    189: /******************* Advect ********************/
    190: typedef enum {
    191:   ADVECT_SOL_TILTED,
    192:   ADVECT_SOL_BUMP
    193: } AdvectSolType;
    194: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "AdvectSolType", "ADVECT_SOL_", 0};
    195: typedef enum {
    196:   ADVECT_SOL_BUMP_CONE,
    197:   ADVECT_SOL_BUMP_COS
    198: } AdvectSolBumpType;
    199: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};

    201: typedef struct {
    202:   PetscReal wind[DIM];
    203: } Physics_Advect_Tilted;
    204: typedef struct {
    205:   PetscReal         center[DIM];
    206:   PetscReal         radius;
    207:   AdvectSolBumpType type;
    208: } Physics_Advect_Bump;

    210: typedef struct {
    211:   PetscReal     inflowState;
    212:   AdvectSolType soltype;
    213:   union
    214:   {
    215:     Physics_Advect_Tilted tilted;
    216:     Physics_Advect_Bump   bump;
    217:   } sol;
    218:   struct {
    219:     PetscInt Error;
    220:   } functional;
    221: } Physics_Advect;

    223: static const struct FieldDescription PhysicsFields_Advect[] = {
    224:   {"U",  1},
    225:   {NULL, 0}
    226: };

    228: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    229: {
    230:   Physics         phys   = (Physics)ctx;
    231:   Physics_Advect *advect = (Physics_Advect *)phys->data;

    233:   PetscFunctionBeginUser;
    234:   xG[0] = advect->inflowState;
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }

    238: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    239: {
    240:   PetscFunctionBeginUser;
    241:   xG[0] = xI[0];
    242:   PetscFunctionReturn(PETSC_SUCCESS);
    243: }

    245: static PetscErrorCode PhysicsRiemann_Advect(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
    246: {
    247:   Physics_Advect *advect = (Physics_Advect *)phys->data;
    248:   PetscReal       wind[DIM], wn;

    250:   PetscFunctionBeginUser;
    251:   switch (advect->soltype) {
    252:   case ADVECT_SOL_TILTED: {
    253:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
    254:     wind[0]                       = tilted->wind[0];
    255:     wind[1]                       = tilted->wind[1];
    256:   } break;
    257:   case ADVECT_SOL_BUMP:
    258:     wind[0] = -qp[1];
    259:     wind[1] = qp[0];
    260:     break;
    261:   default:
    262:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for solution type %s", AdvectSolBumpTypes[advect->soltype]);
    263:   }
    264:   wn      = Dot2(wind, n);
    265:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
    266:   PetscFunctionReturn(PETSC_SUCCESS);
    267: }

    269: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
    270: {
    271:   Physics         phys   = (Physics)ctx;
    272:   Physics_Advect *advect = (Physics_Advect *)phys->data;

    274:   PetscFunctionBeginUser;
    275:   switch (advect->soltype) {
    276:   case ADVECT_SOL_TILTED: {
    277:     PetscReal              x0[DIM];
    278:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
    279:     Waxpy2(-time, tilted->wind, x, x0);
    280:     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
    281:     else u[0] = advect->inflowState;
    282:   } break;
    283:   case ADVECT_SOL_BUMP: {
    284:     Physics_Advect_Bump *bump = &advect->sol.bump;
    285:     PetscReal            x0[DIM], v[DIM], r, cost, sint;
    286:     cost  = PetscCosReal(time);
    287:     sint  = PetscSinReal(time);
    288:     x0[0] = cost * x[0] + sint * x[1];
    289:     x0[1] = -sint * x[0] + cost * x[1];
    290:     Waxpy2(-1, bump->center, x0, v);
    291:     r = Norm2(v);
    292:     switch (bump->type) {
    293:     case ADVECT_SOL_BUMP_CONE:
    294:       u[0] = PetscMax(1 - r / bump->radius, 0);
    295:       break;
    296:     case ADVECT_SOL_BUMP_COS:
    297:       u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
    298:       break;
    299:     }
    300:   } break;
    301:   default:
    302:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
    303:   }
    304:   PetscFunctionReturn(PETSC_SUCCESS);
    305: }

    307: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscScalar *x, const PetscScalar *y, PetscReal *f, void *ctx)
    308: {
    309:   Physics         phys   = (Physics)ctx;
    310:   Physics_Advect *advect = (Physics_Advect *)phys->data;
    311:   PetscScalar     yexact[1];

    313:   PetscFunctionBeginUser;
    314:   PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
    315:   f[advect->functional.Error] = PetscAbsScalar(y[0] - yexact[0]);
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: static PetscErrorCode PhysicsCreate_Advect(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
    320: {
    321:   Physics_Advect *advect;

    323:   PetscFunctionBeginUser;
    324:   phys->field_desc = PhysicsFields_Advect;
    325:   phys->riemann    = (RiemannFunction)PhysicsRiemann_Advect;
    326:   PetscCall(PetscNew(&advect));
    327:   phys->data = advect;
    328:   PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
    329:   {
    330:     PetscInt two = 2, dof = 1;
    331:     advect->soltype = ADVECT_SOL_TILTED;
    332:     PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
    333:     switch (advect->soltype) {
    334:     case ADVECT_SOL_TILTED: {
    335:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
    336:       two                           = 2;
    337:       tilted->wind[0]               = 0.0;
    338:       tilted->wind[1]               = 1.0;
    339:       PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
    340:       advect->inflowState = -2.0;
    341:       PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
    342:       phys->maxspeed = Norm2(tilted->wind);
    343:     } break;
    344:     case ADVECT_SOL_BUMP: {
    345:       Physics_Advect_Bump *bump = &advect->sol.bump;
    346:       two                       = 2;
    347:       bump->center[0]           = 2.;
    348:       bump->center[1]           = 0.;
    349:       PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
    350:       bump->radius = 0.9;
    351:       PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
    352:       bump->type = ADVECT_SOL_BUMP_CONE;
    353:       PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
    354:       phys->maxspeed = 3.; /* radius of mesh, kludge */
    355:     } break;
    356:     }
    357:   }
    358:   PetscOptionsHeadEnd();
    359:   {
    360:     const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
    361:     DMLabel        label;

    363:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
    364:     /* Register "canned" boundary conditions and defaults for where to apply. */
    365:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
    366:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Outflow, NULL, phys, NULL));
    367:     /* Initial/transient solution with default boundary conditions */
    368:     PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys));
    369:     /* Register "canned" functionals */
    370:     PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys));
    371:   }
    372:   PetscFunctionReturn(PETSC_SUCCESS);
    373: }

    375: /******************* Shallow Water ********************/
    376: typedef struct {
    377:   PetscReal gravity;
    378:   PetscReal boundaryHeight;
    379:   struct {
    380:     PetscInt Height;
    381:     PetscInt Speed;
    382:     PetscInt Energy;
    383:   } functional;
    384: } Physics_SW;
    385: typedef struct {
    386:   PetscScalar vals[0];
    387:   PetscScalar h;
    388:   PetscScalar uh[DIM];
    389: } SWNode;

    391: static const struct FieldDescription PhysicsFields_SW[] = {
    392:   {"Height",   1  },
    393:   {"Momentum", DIM},
    394:   {NULL,       0  }
    395: };

    397: /*
    398:  * h_t + div(uh) = 0
    399:  * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
    400:  *
    401:  * */
    402: static PetscErrorCode SWFlux(Physics phys, const PetscReal *n, const SWNode *x, SWNode *f)
    403: {
    404:   Physics_SW *sw = (Physics_SW *)phys->data;
    405:   PetscScalar uhn, u[DIM];
    406:   PetscInt    i;

    408:   PetscFunctionBeginUser;
    409:   Scale2(1. / x->h, x->uh, u);
    410:   uhn  = Dot2(x->uh, n);
    411:   f->h = uhn;
    412:   for (i = 0; i < DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
    413:   PetscFunctionReturn(PETSC_SUCCESS);
    414: }

    416: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    417: {
    418:   PetscFunctionBeginUser;
    419:   xG[0] = xI[0];
    420:   xG[1] = -xI[1];
    421:   xG[2] = -xI[2];
    422:   PetscFunctionReturn(PETSC_SUCCESS);
    423: }

    425: static PetscErrorCode PhysicsRiemann_SW(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
    426: {
    427:   Physics_SW   *sw = (Physics_SW *)phys->data;
    428:   PetscReal     cL, cR, speed, nn[DIM];
    429:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
    430:   SWNode        fL, fR;
    431:   PetscInt      i;

    433:   PetscFunctionBeginUser;
    434:   PetscCheck(uL->h >= 0 && uR->h >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
    435:   nn[0] = n[0];
    436:   nn[1] = n[1];
    437:   Normalize2(nn);
    438:   SWFlux(phys, nn, uL, &fL);
    439:   SWFlux(phys, nn, uR, &fR);
    440:   cL    = PetscSqrtReal(sw->gravity * PetscRealPart(uL->h));
    441:   cR    = PetscSqrtReal(sw->gravity * PetscRealPart(uR->h)); /* gravity wave speed */
    442:   speed = PetscMax(PetscAbsScalar(Dot2(uL->uh, nn) / uL->h) + cL, PetscAbsScalar(Dot2(uR->uh, nn) / uR->h) + cR);
    443:   for (i = 0; i < 1 + DIM; i++) flux[i] = (0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i])) * Norm2(n);
    444:   PetscFunctionReturn(PETSC_SUCCESS);
    445: }

    447: static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
    448: {
    449:   PetscReal dx[2], r, sigma;

    451:   PetscFunctionBeginUser;
    452:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
    453:   dx[0] = x[0] - 1.5;
    454:   dx[1] = x[1] - 1.0;
    455:   r     = Norm2(dx);
    456:   sigma = 0.5;
    457:   u[0]  = 1 + 2 * PetscExpScalar(-PetscSqr(r) / (2 * PetscSqr(sigma)));
    458:   u[1]  = 0.0;
    459:   u[2]  = 0.0;
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }

    463: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
    464: {
    465:   Physics       phys = (Physics)ctx;
    466:   Physics_SW   *sw   = (Physics_SW *)phys->data;
    467:   const SWNode *x    = (const SWNode *)xx;
    468:   PetscScalar   u[2];
    469:   PetscReal     h;

    471:   PetscFunctionBeginUser;
    472:   h = PetscRealPart(x->h);
    473:   Scale2(1. / x->h, x->uh, u);
    474:   f[sw->functional.Height] = h;
    475:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity * h);
    476:   f[sw->functional.Energy] = 0.5 * (Dot2(x->uh, u) + sw->gravity * PetscSqr(h));
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

    480: static PetscErrorCode PhysicsCreate_SW(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
    481: {
    482:   Physics_SW *sw;

    484:   PetscFunctionBeginUser;
    485:   phys->field_desc = PhysicsFields_SW;
    486:   phys->riemann    = (RiemannFunction)PhysicsRiemann_SW;
    487:   PetscCall(PetscNew(&sw));
    488:   phys->data = sw;
    489:   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
    490:   {
    491:     sw->gravity = 1.0;
    492:     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
    493:   }
    494:   PetscOptionsHeadEnd();
    495:   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */

    497:   {
    498:     const PetscInt wallids[] = {100, 101, 200, 300};
    499:     DMLabel        label;

    501:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
    502:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_SW_Wall, NULL, phys, NULL));
    503:     PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
    504:     PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
    505:     PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
    506:     PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
    507:   }
    508:   PetscFunctionReturn(PETSC_SUCCESS);
    509: }

    511: /******************* Euler ********************/
    512: typedef struct {
    513:   PetscScalar vals[0];
    514:   PetscScalar r;
    515:   PetscScalar ru[DIM];
    516:   PetscScalar e;
    517: } EulerNode;
    518: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode *, PetscScalar *);
    519: typedef struct {
    520:   PetscInt        npars;
    521:   PetscReal       pars[DIM];
    522:   EquationOfState pressure;
    523:   EquationOfState sound;
    524:   struct {
    525:     PetscInt Density;
    526:     PetscInt Momentum;
    527:     PetscInt Energy;
    528:     PetscInt Pressure;
    529:     PetscInt Speed;
    530:   } monitor;
    531: } Physics_Euler;

    533: static const struct FieldDescription PhysicsFields_Euler[] = {
    534:   {"Density",  1  },
    535:   {"Momentum", DIM},
    536:   {"Energy",   1  },
    537:   {NULL,       0  }
    538: };

    540: static PetscErrorCode Pressure_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *p)
    541: {
    542:   PetscScalar ru2;

    544:   PetscFunctionBeginUser;
    545:   ru2 = DotDIM(x->ru, x->ru);
    546:   ru2 /= x->r;
    547:   /* kinematic dof = params[0] */
    548:   (*p) = 2.0 * (x->e - 0.5 * ru2) / pars[0];
    549:   PetscFunctionReturn(PETSC_SUCCESS);
    550: }

    552: static PetscErrorCode SpeedOfSound_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *c)
    553: {
    554:   PetscScalar p;

    556:   PetscFunctionBeginUser;
    557:   /* TODO remove direct usage of Pressure_PG */
    558:   Pressure_PG(pars, x, &p);
    559:   /* TODO check the sign of p */
    560:   /* pars[1] = heat capacity ratio */
    561:   (*c) = PetscSqrtScalar(pars[1] * p / x->r);
    562:   PetscFunctionReturn(PETSC_SUCCESS);
    563: }

    565: /*
    566:  * x = (rho,rho*(u_1),...,rho*e)^T
    567:  * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
    568:  *
    569:  * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
    570:  *
    571:  * */
    572: static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
    573: {
    574:   Physics_Euler *eu = (Physics_Euler *)phys->data;
    575:   PetscScalar    u, nu, p;
    576:   PetscInt       i;

    578:   PetscFunctionBeginUser;
    579:   u = DotDIM(x->ru, x->ru);
    580:   u /= (x->r * x->r);
    581:   nu = DotDIM(x->ru, n);
    582:   /* TODO check the sign of p */
    583:   eu->pressure(eu->pars, x, &p);
    584:   f->r = nu * x->r;
    585:   for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p;
    586:   f->e = nu * (x->e + p);
    587:   PetscFunctionReturn(PETSC_SUCCESS);
    588: }

    590: /* PetscReal* => EulerNode* conversion */
    591: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
    592: {
    593:   PetscInt    i;
    594:   PetscScalar xn[DIM], xt[DIM];

    596:   PetscFunctionBeginUser;
    597:   xG[0] = xI[0];
    598:   NormalSplitDIM(n, xI + 1, xn, xt);
    599:   for (i = 0; i < DIM; i++) xG[i + 1] = -xn[i] + xt[i];
    600:   xG[DIM + 1] = xI[DIM + 1];
    601:   PetscFunctionReturn(PETSC_SUCCESS);
    602: }

    604: /* PetscReal* => EulerNode* conversion */
    605: static PetscErrorCode PhysicsRiemann_Euler_Rusanov(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
    606: {
    607:   Physics_Euler   *eu = (Physics_Euler *)phys->data;
    608:   PetscScalar      cL, cR, speed;
    609:   const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
    610:   EulerNode        fL, fR;
    611:   PetscInt         i;

    613:   PetscFunctionBeginUser;
    614:   PetscCheck(uL->r >= 0 && uR->r >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
    615:   EulerFlux(phys, n, uL, &fL);
    616:   EulerFlux(phys, n, uR, &fR);
    617:   eu->sound(eu->pars, uL, &cL);
    618:   eu->sound(eu->pars, uR, &cR);
    619:   speed = PetscMax(cL, cR) + PetscMax(PetscAbsScalar(DotDIM(uL->ru, n) / NormDIM(n)), PetscAbsScalar(DotDIM(uR->ru, n) / NormDIM(n)));
    620:   for (i = 0; i < 2 + DIM; i++) flux[i] = 0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i]);
    621:   PetscFunctionReturn(PETSC_SUCCESS);
    622: }

    624: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
    625: {
    626:   PetscInt i;

    628:   PetscFunctionBeginUser;
    629:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
    630:   u[0]       = 1.0;
    631:   u[DIM + 1] = 1.0 + PetscAbsReal(x[0]);
    632:   for (i = 1; i < DIM + 1; i++) u[i] = 0.0;
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
    637: {
    638:   Physics          phys = (Physics)ctx;
    639:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
    640:   const EulerNode *x    = (const EulerNode *)xx;
    641:   PetscScalar      p;

    643:   PetscFunctionBeginUser;
    644:   f[eu->monitor.Density]  = x->r;
    645:   f[eu->monitor.Momentum] = NormDIM(x->ru);
    646:   f[eu->monitor.Energy]   = x->e;
    647:   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
    648:   eu->pressure(eu->pars, x, &p);
    649:   f[eu->monitor.Pressure] = p;
    650:   PetscFunctionReturn(PETSC_SUCCESS);
    651: }

    653: static PetscErrorCode PhysicsCreate_Euler(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
    654: {
    655:   Physics_Euler *eu;

    657:   PetscFunctionBeginUser;
    658:   phys->field_desc = PhysicsFields_Euler;
    659:   phys->riemann    = (RiemannFunction)PhysicsRiemann_Euler_Rusanov;
    660:   PetscCall(PetscNew(&eu));
    661:   phys->data = eu;
    662:   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
    663:   {
    664:     eu->pars[0] = 3.0;
    665:     eu->pars[1] = 1.67;
    666:     PetscCall(PetscOptionsReal("-eu_f", "Degrees of freedom", "", eu->pars[0], &eu->pars[0], NULL));
    667:     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[1], &eu->pars[1], NULL));
    668:   }
    669:   PetscOptionsHeadEnd();
    670:   eu->pressure   = Pressure_PG;
    671:   eu->sound      = SpeedOfSound_PG;
    672:   phys->maxspeed = 1.0;
    673:   {
    674:     const PetscInt wallids[] = {100, 101, 200, 300};
    675:     DMLabel        label;

    677:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
    678:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
    679:     PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
    680:     PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
    681:     PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
    682:     PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
    683:     PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
    684:     PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
    685:   }
    686:   PetscFunctionReturn(PETSC_SUCCESS);
    687: }

    689: PetscErrorCode ConstructCellBoundary(DM dm, User user)
    690: {
    691:   const char     *name   = "Cell Sets";
    692:   const char     *bdname = "split faces";
    693:   IS              regionIS, innerIS;
    694:   const PetscInt *regions, *cells;
    695:   PetscInt        numRegions, innerRegion, numCells, c;
    696:   PetscInt        cStart, cEnd, cEndInterior, fStart, fEnd;
    697:   PetscBool       hasLabel;

    699:   PetscFunctionBeginUser;
    700:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    701:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    702:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));

    704:   PetscCall(DMHasLabel(dm, name, &hasLabel));
    705:   if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
    706:   PetscCall(DMGetLabelSize(dm, name, &numRegions));
    707:   if (numRegions != 2) PetscFunctionReturn(PETSC_SUCCESS);
    708:   /* Get the inner id */
    709:   PetscCall(DMGetLabelIdIS(dm, name, &regionIS));
    710:   PetscCall(ISGetIndices(regionIS, &regions));
    711:   innerRegion = regions[0];
    712:   PetscCall(ISRestoreIndices(regionIS, &regions));
    713:   PetscCall(ISDestroy(&regionIS));
    714:   /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
    715:   PetscCall(DMGetStratumIS(dm, name, innerRegion, &innerIS));
    716:   PetscCall(ISGetLocalSize(innerIS, &numCells));
    717:   PetscCall(ISGetIndices(innerIS, &cells));
    718:   PetscCall(DMCreateLabel(dm, bdname));
    719:   for (c = 0; c < numCells; ++c) {
    720:     const PetscInt  cell = cells[c];
    721:     const PetscInt *faces;
    722:     PetscInt        numFaces, f;

    724:     PetscCheck((cell >= cStart) && (cell < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", cell);
    725:     PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
    726:     PetscCall(DMPlexGetCone(dm, cell, &faces));
    727:     for (f = 0; f < numFaces; ++f) {
    728:       const PetscInt  face = faces[f];
    729:       const PetscInt *neighbors;
    730:       PetscInt        nC, regionA, regionB;

    732:       PetscCheck((face >= fStart) && (face < fEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a face", face);
    733:       PetscCall(DMPlexGetSupportSize(dm, face, &nC));
    734:       if (nC != 2) continue;
    735:       PetscCall(DMPlexGetSupport(dm, face, &neighbors));
    736:       if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue;
    737:       PetscCheck((neighbors[0] >= cStart) && (neighbors[0] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[0]);
    738:       PetscCheck((neighbors[1] >= cStart) && (neighbors[1] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[1]);
    739:       PetscCall(DMGetLabelValue(dm, name, neighbors[0], &regionA));
    740:       PetscCall(DMGetLabelValue(dm, name, neighbors[1], &regionB));
    741:       PetscCheck(regionA >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[0]);
    742:       PetscCheck(regionB >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[1]);
    743:       if (regionA != regionB) PetscCall(DMSetLabelValue(dm, bdname, faces[f], 1));
    744:     }
    745:   }
    746:   PetscCall(ISRestoreIndices(innerIS, &cells));
    747:   PetscCall(ISDestroy(&innerIS));
    748:   {
    749:     DMLabel label;

    751:     PetscCall(DMGetLabel(dm, bdname, &label));
    752:     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
    753:   }
    754:   PetscFunctionReturn(PETSC_SUCCESS);
    755: }

    757: /* Right now, I have just added duplicate faces, which see both cells. We can
    758: - Add duplicate vertices and decouple the face cones
    759: - Disconnect faces from cells across the rotation gap
    760: */
    761: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)
    762: {
    763:   DM              dm = *dmSplit, sdm;
    764:   PetscSF         sfPoint, gsfPoint;
    765:   PetscSection    coordSection, newCoordSection;
    766:   Vec             coordinates;
    767:   IS              idIS;
    768:   const PetscInt *ids;
    769:   PetscInt       *newpoints;
    770:   PetscInt        dim, depth, maxConeSize, maxSupportSize, numLabels, numGhostCells;
    771:   PetscInt        numFS, fs, pStart, pEnd, p, cEnd, cEndInterior, vStart, vEnd, v, fStart, fEnd, newf, d, l;
    772:   PetscBool       hasLabel;

    774:   PetscFunctionBeginUser;
    775:   PetscCall(DMHasLabel(dm, labelName, &hasLabel));
    776:   if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
    777:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm));
    778:   PetscCall(DMSetType(sdm, DMPLEX));
    779:   PetscCall(DMGetDimension(dm, &dim));
    780:   PetscCall(DMSetDimension(sdm, dim));

    782:   PetscCall(DMGetLabelIdIS(dm, labelName, &idIS));
    783:   PetscCall(ISGetLocalSize(idIS, &numFS));
    784:   PetscCall(ISGetIndices(idIS, &ids));

    786:   user->numSplitFaces = 0;
    787:   for (fs = 0; fs < numFS; ++fs) {
    788:     PetscInt numBdFaces;

    790:     PetscCall(DMGetStratumSize(dm, labelName, ids[fs], &numBdFaces));
    791:     user->numSplitFaces += numBdFaces;
    792:   }
    793:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    794:   pEnd += user->numSplitFaces;
    795:   PetscCall(DMPlexSetChart(sdm, pStart, pEnd));
    796:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
    797:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
    798:   numGhostCells = cEnd - cEndInterior;
    799:   /* Set cone and support sizes */
    800:   PetscCall(DMPlexGetDepth(dm, &depth));
    801:   for (d = 0; d <= depth; ++d) {
    802:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    803:     for (p = pStart; p < pEnd; ++p) {
    804:       PetscInt newp = p;
    805:       PetscInt size;

    807:       PetscCall(DMPlexGetConeSize(dm, p, &size));
    808:       PetscCall(DMPlexSetConeSize(sdm, newp, size));
    809:       PetscCall(DMPlexGetSupportSize(dm, p, &size));
    810:       PetscCall(DMPlexSetSupportSize(sdm, newp, size));
    811:     }
    812:   }
    813:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    814:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
    815:     IS              faceIS;
    816:     const PetscInt *faces;
    817:     PetscInt        numFaces, f;

    819:     PetscCall(DMGetStratumIS(dm, labelName, ids[fs], &faceIS));
    820:     PetscCall(ISGetLocalSize(faceIS, &numFaces));
    821:     PetscCall(ISGetIndices(faceIS, &faces));
    822:     for (f = 0; f < numFaces; ++f, ++newf) {
    823:       PetscInt size;

    825:       /* Right now I think that both faces should see both cells */
    826:       PetscCall(DMPlexGetConeSize(dm, faces[f], &size));
    827:       PetscCall(DMPlexSetConeSize(sdm, newf, size));
    828:       PetscCall(DMPlexGetSupportSize(dm, faces[f], &size));
    829:       PetscCall(DMPlexSetSupportSize(sdm, newf, size));
    830:     }
    831:     PetscCall(ISRestoreIndices(faceIS, &faces));
    832:     PetscCall(ISDestroy(&faceIS));
    833:   }
    834:   PetscCall(DMSetUp(sdm));
    835:   /* Set cones and supports */
    836:   PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
    837:   PetscCall(PetscMalloc1(PetscMax(maxConeSize, maxSupportSize), &newpoints));
    838:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    839:   for (p = pStart; p < pEnd; ++p) {
    840:     const PetscInt *points, *orientations;
    841:     PetscInt        size, i, newp = p;

    843:     PetscCall(DMPlexGetConeSize(dm, p, &size));
    844:     PetscCall(DMPlexGetCone(dm, p, &points));
    845:     PetscCall(DMPlexGetConeOrientation(dm, p, &orientations));
    846:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
    847:     PetscCall(DMPlexSetCone(sdm, newp, newpoints));
    848:     PetscCall(DMPlexSetConeOrientation(sdm, newp, orientations));
    849:     PetscCall(DMPlexGetSupportSize(dm, p, &size));
    850:     PetscCall(DMPlexGetSupport(dm, p, &points));
    851:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
    852:     PetscCall(DMPlexSetSupport(sdm, newp, newpoints));
    853:   }
    854:   PetscCall(PetscFree(newpoints));
    855:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
    856:     IS              faceIS;
    857:     const PetscInt *faces;
    858:     PetscInt        numFaces, f;

    860:     PetscCall(DMGetStratumIS(dm, labelName, ids[fs], &faceIS));
    861:     PetscCall(ISGetLocalSize(faceIS, &numFaces));
    862:     PetscCall(ISGetIndices(faceIS, &faces));
    863:     for (f = 0; f < numFaces; ++f, ++newf) {
    864:       const PetscInt *points;

    866:       PetscCall(DMPlexGetCone(dm, faces[f], &points));
    867:       PetscCall(DMPlexSetCone(sdm, newf, points));
    868:       PetscCall(DMPlexGetSupport(dm, faces[f], &points));
    869:       PetscCall(DMPlexSetSupport(sdm, newf, points));
    870:     }
    871:     PetscCall(ISRestoreIndices(faceIS, &faces));
    872:     PetscCall(ISDestroy(&faceIS));
    873:   }
    874:   PetscCall(ISRestoreIndices(idIS, &ids));
    875:   PetscCall(ISDestroy(&idIS));
    876:   PetscCall(DMPlexStratify(sdm));
    877:   /* Convert coordinates */
    878:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    879:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    880:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection));
    881:   PetscCall(PetscSectionSetNumFields(newCoordSection, 1));
    882:   PetscCall(PetscSectionSetFieldComponents(newCoordSection, 0, dim));
    883:   PetscCall(PetscSectionSetChart(newCoordSection, vStart, vEnd));
    884:   for (v = vStart; v < vEnd; ++v) {
    885:     PetscCall(PetscSectionSetDof(newCoordSection, v, dim));
    886:     PetscCall(PetscSectionSetFieldDof(newCoordSection, v, 0, dim));
    887:   }
    888:   PetscCall(PetscSectionSetUp(newCoordSection));
    889:   PetscCall(DMSetCoordinateSection(sdm, PETSC_DETERMINE, newCoordSection));
    890:   PetscCall(PetscSectionDestroy(&newCoordSection)); /* relinquish our reference */
    891:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    892:   PetscCall(DMSetCoordinatesLocal(sdm, coordinates));
    893:   /* Convert labels */
    894:   PetscCall(DMGetNumLabels(dm, &numLabels));
    895:   for (l = 0; l < numLabels; ++l) {
    896:     const char *lname;
    897:     PetscBool   isDepth;

    899:     PetscCall(DMGetLabelName(dm, l, &lname));
    900:     PetscCall(PetscStrcmp(lname, "depth", &isDepth));
    901:     if (isDepth) continue;
    902:     PetscCall(DMCreateLabel(sdm, lname));
    903:     PetscCall(DMGetLabelIdIS(dm, lname, &idIS));
    904:     PetscCall(ISGetLocalSize(idIS, &numFS));
    905:     PetscCall(ISGetIndices(idIS, &ids));
    906:     for (fs = 0; fs < numFS; ++fs) {
    907:       IS              pointIS;
    908:       const PetscInt *points;
    909:       PetscInt        numPoints;

    911:       PetscCall(DMGetStratumIS(dm, lname, ids[fs], &pointIS));
    912:       PetscCall(ISGetLocalSize(pointIS, &numPoints));
    913:       PetscCall(ISGetIndices(pointIS, &points));
    914:       for (p = 0; p < numPoints; ++p) {
    915:         PetscInt newpoint = points[p];

    917:         PetscCall(DMSetLabelValue(sdm, lname, newpoint, ids[fs]));
    918:       }
    919:       PetscCall(ISRestoreIndices(pointIS, &points));
    920:       PetscCall(ISDestroy(&pointIS));
    921:     }
    922:     PetscCall(ISRestoreIndices(idIS, &ids));
    923:     PetscCall(ISDestroy(&idIS));
    924:   }
    925:   /* Convert pointSF */
    926:   const PetscSFNode *remotePoints;
    927:   PetscSFNode       *gremotePoints;
    928:   const PetscInt    *localPoints;
    929:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
    930:   PetscInt           numRoots, numLeaves;
    931:   PetscMPIInt        size;

    933:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    934:   PetscCall(DMGetPointSF(dm, &sfPoint));
    935:   PetscCall(DMGetPointSF(sdm, &gsfPoint));
    936:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    937:   PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
    938:   if (numRoots >= 0) {
    939:     PetscCall(PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation));
    940:     for (l = 0; l < numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? numGhostCells : 0); */
    941:     PetscCall(PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
    942:     PetscCall(PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
    943:     PetscCall(PetscMalloc1(numLeaves, &glocalPoints));
    944:     PetscCall(PetscMalloc1(numLeaves, &gremotePoints));
    945:     for (l = 0; l < numLeaves; ++l) {
    946:       glocalPoints[l]        = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + numGhostCells : localPoints[l]; */
    947:       gremotePoints[l].rank  = remotePoints[l].rank;
    948:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
    949:     }
    950:     PetscCall(PetscFree2(newLocation, newRemoteLocation));
    951:     PetscCall(PetscSFSetGraph(gsfPoint, numRoots + numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER));
    952:   }
    953:   PetscCall(DMDestroy(dmSplit));
    954:   *dmSplit = sdm;
    955:   PetscFunctionReturn(PETSC_SUCCESS);
    956: }

    958: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
    959: {
    960:   PetscSF      sfPoint;
    961:   PetscSection coordSection;
    962:   Vec          coordinates;
    963:   PetscSection sectionCell;
    964:   PetscScalar *part;
    965:   PetscInt     cStart, cEnd, c;
    966:   PetscMPIInt  rank;

    968:   PetscFunctionBeginUser;
    969:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    970:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    971:   PetscCall(DMClone(dm, dmCell));
    972:   PetscCall(DMGetPointSF(dm, &sfPoint));
    973:   PetscCall(DMSetPointSF(*dmCell, sfPoint));
    974:   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
    975:   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
    976:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    977:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
    978:   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
    979:   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
    980:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
    981:   PetscCall(PetscSectionSetUp(sectionCell));
    982:   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
    983:   PetscCall(PetscSectionDestroy(&sectionCell));
    984:   PetscCall(DMCreateLocalVector(*dmCell, partition));
    985:   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
    986:   PetscCall(VecGetArray(*partition, &part));
    987:   for (c = cStart; c < cEnd; ++c) {
    988:     PetscScalar *p;

    990:     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
    991:     p[0] = rank;
    992:   }
    993:   PetscCall(VecRestoreArray(*partition, &part));
    994:   PetscFunctionReturn(PETSC_SUCCESS);
    995: }

    997: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
    998: {
    999:   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
    1000:   PetscSection       coordSection;
    1001:   Vec                coordinates, facegeom, cellgeom;
    1002:   PetscSection       sectionMass;
    1003:   PetscScalar       *m;
    1004:   const PetscScalar *fgeom, *cgeom, *coords;
    1005:   PetscInt           vStart, vEnd, v;

    1007:   PetscFunctionBeginUser;
    1008:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1009:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    1010:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    1011:   PetscCall(DMClone(dm, &dmMass));
    1012:   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
    1013:   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
    1014:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
    1015:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    1016:   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
    1017:   for (v = vStart; v < vEnd; ++v) {
    1018:     PetscInt numFaces;

    1020:     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
    1021:     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
    1022:   }
    1023:   PetscCall(PetscSectionSetUp(sectionMass));
    1024:   PetscCall(DMSetLocalSection(dmMass, sectionMass));
    1025:   PetscCall(PetscSectionDestroy(&sectionMass));
    1026:   PetscCall(DMGetLocalVector(dmMass, massMatrix));
    1027:   PetscCall(VecGetArray(*massMatrix, &m));
    1028:   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
    1029:   PetscCall(VecGetDM(facegeom, &dmFace));
    1030:   PetscCall(VecGetArrayRead(facegeom, &fgeom));
    1031:   PetscCall(VecGetDM(cellgeom, &dmCell));
    1032:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    1033:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
    1034:   PetscCall(VecGetArrayRead(coordinates, &coords));
    1035:   for (v = vStart; v < vEnd; ++v) {
    1036:     const PetscInt        *faces;
    1037:     const PetscFVFaceGeom *fgA, *fgB, *cg;
    1038:     const PetscScalar     *vertex;
    1039:     PetscInt               numFaces, sides[2], f, g;

    1041:     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
    1042:     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
    1043:     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
    1044:     for (f = 0; f < numFaces; ++f) {
    1045:       sides[0] = faces[f];
    1046:       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
    1047:       for (g = 0; g < numFaces; ++g) {
    1048:         const PetscInt *cells = NULL;
    1049:         PetscReal       area  = 0.0;
    1050:         PetscInt        numCells;

    1052:         sides[1] = faces[g];
    1053:         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
    1054:         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
    1055:         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
    1056:         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
    1057:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
    1058:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
    1059:         m[f * numFaces + g] = Dot2(fgA->normal, fgB->normal) * area * 0.5;
    1060:         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
    1061:       }
    1062:     }
    1063:   }
    1064:   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
    1065:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    1066:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
    1067:   PetscCall(VecRestoreArray(*massMatrix, &m));
    1068:   PetscCall(DMDestroy(&dmMass));
    1069:   PetscCall(DMDestroy(&plex));
    1070:   PetscFunctionReturn(PETSC_SUCCESS);
    1071: }

    1073: PetscErrorCode SetUpLocalSpace(DM dm, User user)
    1074: {
    1075:   PetscSection stateSection;
    1076:   Physics      phys;
    1077:   PetscInt     dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, cEndInterior, c, i;

    1079:   PetscFunctionBeginUser;
    1080:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    1081:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
    1082:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection));
    1083:   phys = user->model->physics;
    1084:   PetscCall(PetscSectionSetNumFields(stateSection, phys->nfields));
    1085:   for (i = 0; i < phys->nfields; i++) {
    1086:     PetscCall(PetscSectionSetFieldName(stateSection, i, phys->field_desc[i].name));
    1087:     PetscCall(PetscSectionSetFieldComponents(stateSection, i, phys->field_desc[i].dof));
    1088:   }
    1089:   PetscCall(PetscSectionSetChart(stateSection, cStart, cEnd));
    1090:   for (c = cStart; c < cEnd; ++c) {
    1091:     for (i = 0; i < phys->nfields; i++) PetscCall(PetscSectionSetFieldDof(stateSection, c, i, phys->field_desc[i].dof));
    1092:     PetscCall(PetscSectionSetDof(stateSection, c, dof));
    1093:   }
    1094:   for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintDof(stateSection, c, dof));
    1095:   PetscCall(PetscSectionSetUp(stateSection));
    1096:   PetscCall(PetscMalloc1(dof, &cind));
    1097:   for (d = 0; d < dof; ++d) cind[d] = d;
    1098: #if 0
    1099:   for (c = cStart; c < cEnd; ++c) {
    1100:     PetscInt val;

    1102:     PetscCall(DMGetLabelValue(dm, "vtk", c, &val));
    1103:     if (val < 0) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
    1104:   }
    1105: #endif
    1106:   for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
    1107:   PetscCall(PetscFree(cind));
    1108:   PetscCall(PetscSectionGetStorageSize(stateSection, &stateSize));
    1109:   PetscCall(DMSetLocalSection(dm, stateSection));
    1110:   PetscCall(PetscSectionDestroy(&stateSection));
    1111:   PetscFunctionReturn(PETSC_SUCCESS);
    1112: }

    1114: #if 0
    1115: PetscErrorCode SetUpBoundaries(DM dm, User user)
    1116: {
    1117:   Model          mod = user->model;
    1118:   BoundaryLink   b;

    1120:   PetscFunctionBeginUser;
    1121:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm),NULL,"Boundary condition options","");
    1122:   for (b = mod->boundary; b; b=b->next) {
    1123:     char      optname[512];
    1124:     PetscInt  ids[512],len = 512;
    1125:     PetscBool flg;
    1126:     PetscCall(PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name));
    1127:     PetscCall(PetscMemzero(ids,sizeof(ids)));
    1128:     PetscCall(PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg));
    1129:     if (flg) {
    1130:       /* TODO: check all IDs to make sure they exist in the mesh */
    1131:       PetscCall(PetscFree(b->ids));
    1132:       b->numids = len;
    1133:       PetscCall(PetscMalloc1(len,&b->ids));
    1134:       PetscCall(PetscArraycpy(b->ids,ids,len));
    1135:     }
    1136:   }
    1137:   PetscOptionsEnd();
    1138:   PetscFunctionReturn(PETSC_SUCCESS);
    1139: }
    1140: #endif

    1142: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
    1143: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
    1144: {
    1145:   PetscFunctionBeginUser;
    1146:   mod->solution    = func;
    1147:   mod->solutionctx = ctx;
    1148:   PetscFunctionReturn(PETSC_SUCCESS);
    1149: }

    1151: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
    1152: {
    1153:   FunctionalLink link, *ptr;
    1154:   PetscInt       lastoffset = -1;

    1156:   PetscFunctionBeginUser;
    1157:   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
    1158:   PetscCall(PetscNew(&link));
    1159:   PetscCall(PetscStrallocpy(name, &link->name));
    1160:   link->offset = lastoffset + 1;
    1161:   link->func   = func;
    1162:   link->ctx    = ctx;
    1163:   link->next   = NULL;
    1164:   *ptr         = link;
    1165:   *offset      = link->offset;
    1166:   PetscFunctionReturn(PETSC_SUCCESS);
    1167: }

    1169: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptions *PetscOptionsObject)
    1170: {
    1171:   PetscInt       i, j;
    1172:   FunctionalLink link;
    1173:   char          *names[256];

    1175:   PetscFunctionBeginUser;
    1176:   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
    1177:   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
    1178:   /* Create list of functionals that will be computed somehow */
    1179:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
    1180:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
    1181:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
    1182:   mod->numCall = 0;
    1183:   for (i = 0; i < mod->numMonitored; i++) {
    1184:     for (link = mod->functionalRegistry; link; link = link->next) {
    1185:       PetscBool match;
    1186:       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
    1187:       if (match) break;
    1188:     }
    1189:     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
    1190:     mod->functionalMonitored[i] = link;
    1191:     for (j = 0; j < i; j++) {
    1192:       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
    1193:     }
    1194:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
    1195:   next_name:
    1196:     PetscCall(PetscFree(names[i]));
    1197:   }

    1199:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
    1200:   mod->maxComputed = -1;
    1201:   for (link = mod->functionalRegistry; link; link = link->next) {
    1202:     for (i = 0; i < mod->numCall; i++) {
    1203:       FunctionalLink call = mod->functionalCall[i];
    1204:       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
    1205:     }
    1206:   }
    1207:   PetscFunctionReturn(PETSC_SUCCESS);
    1208: }

    1210: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
    1211: {
    1212:   FunctionalLink l, next;

    1214:   PetscFunctionBeginUser;
    1215:   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
    1216:   l     = *link;
    1217:   *link = NULL;
    1218:   for (; l; l = next) {
    1219:     next = l->next;
    1220:     PetscCall(PetscFree(l->name));
    1221:     PetscCall(PetscFree(l));
    1222:   }
    1223:   PetscFunctionReturn(PETSC_SUCCESS);
    1224: }

    1226: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
    1227: {
    1228:   DM                 plex, dmCell;
    1229:   Model              mod = user->model;
    1230:   Vec                cellgeom;
    1231:   const PetscScalar *cgeom;
    1232:   PetscScalar       *x;
    1233:   PetscInt           cStart, cEnd, c;

    1235:   PetscFunctionBeginUser;
    1236:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1237:   PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
    1238:   PetscCall(VecGetDM(cellgeom, &dmCell));
    1239:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    1240:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    1241:   PetscCall(VecGetArray(X, &x));
    1242:   for (c = cStart; c < cEnd; ++c) {
    1243:     const PetscFVCellGeom *cg;
    1244:     PetscScalar           *xc;

    1246:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
    1247:     PetscCall(DMPlexPointGlobalRef(dm, c, x, &xc));
    1248:     if (xc) PetscCall((*mod->solution)(mod, 0.0, cg->centroid, xc, mod->solutionctx));
    1249:   }
    1250:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    1251:   PetscCall(VecRestoreArray(X, &x));
    1252:   PetscCall(DMDestroy(&plex));
    1253:   PetscFunctionReturn(PETSC_SUCCESS);
    1254: }

    1256: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
    1257: {
    1258:   PetscFunctionBeginUser;
    1259:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
    1260:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
    1261:   PetscCall(PetscViewerFileSetName(*viewer, filename));
    1262:   PetscFunctionReturn(PETSC_SUCCESS);
    1263: }

    1265: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
    1266: {
    1267:   User        user = (User)ctx;
    1268:   DM          dm, plex;
    1269:   PetscViewer viewer;
    1270:   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
    1271:   PetscReal   xnorm;

    1273:   PetscFunctionBeginUser;
    1274:   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
    1275:   PetscCall(VecGetDM(X, &dm));
    1276:   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
    1277:   if (stepnum >= 0) { /* No summary for final time */
    1278:     Model              mod = user->model;
    1279:     Vec                cellgeom;
    1280:     PetscInt           c, cStart, cEnd, fcount, i;
    1281:     size_t             ftableused, ftablealloc;
    1282:     const PetscScalar *cgeom, *x;
    1283:     DM                 dmCell;
    1284:     PetscReal         *fmin, *fmax, *fintegral, *ftmp;

    1286:     PetscCall(DMConvert(dm, DMPLEX, &plex));
    1287:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
    1288:     fcount = mod->maxComputed + 1;
    1289:     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
    1290:     for (i = 0; i < fcount; i++) {
    1291:       fmin[i]      = PETSC_MAX_REAL;
    1292:       fmax[i]      = PETSC_MIN_REAL;
    1293:       fintegral[i] = 0;
    1294:     }
    1295:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    1296:     PetscCall(VecGetDM(cellgeom, &dmCell));
    1297:     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    1298:     PetscCall(VecGetArrayRead(X, &x));
    1299:     for (c = cStart; c < cEnd; ++c) {
    1300:       const PetscFVCellGeom *cg;
    1301:       const PetscScalar     *cx;
    1302:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
    1303:       PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
    1304:       if (!cx) continue; /* not a global cell */
    1305:       for (i = 0; i < mod->numCall; i++) {
    1306:         FunctionalLink flink = mod->functionalCall[i];
    1307:         PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
    1308:       }
    1309:       for (i = 0; i < fcount; i++) {
    1310:         fmin[i] = PetscMin(fmin[i], ftmp[i]);
    1311:         fmax[i] = PetscMax(fmax[i], ftmp[i]);
    1312:         fintegral[i] += cg->volume * ftmp[i];
    1313:       }
    1314:     }
    1315:     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    1316:     PetscCall(VecRestoreArrayRead(X, &x));
    1317:     PetscCall(DMDestroy(&plex));
    1318:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
    1319:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
    1320:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));

    1322:     ftablealloc = fcount * 100;
    1323:     ftableused  = 0;
    1324:     PetscCall(PetscMalloc1(ftablealloc, &ftable));
    1325:     for (i = 0; i < mod->numMonitored; i++) {
    1326:       size_t         countused;
    1327:       char           buffer[256], *p;
    1328:       FunctionalLink flink = mod->functionalMonitored[i];
    1329:       PetscInt       id    = flink->offset;
    1330:       if (i % 3) {
    1331:         PetscCall(PetscArraycpy(buffer, "  ", 2));
    1332:         p = buffer + 2;
    1333:       } else if (i) {
    1334:         char newline[] = "\n";
    1335:         PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
    1336:         p = buffer + sizeof(newline) - 1;
    1337:       } else {
    1338:         p = buffer;
    1339:       }
    1340:       PetscCall(PetscSNPrintfCount(p, sizeof buffer - (p - buffer), "%12s [%10.7g,%10.7g] int %10.7g", &countused, flink->name, (double)fmin[id], (double)fmax[id], (double)fintegral[id]));
    1341:       countused += p - buffer;
    1342:       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
    1343:         char *ftablenew;
    1344:         ftablealloc = 2 * ftablealloc + countused;
    1345:         PetscCall(PetscMalloc(ftablealloc, &ftablenew));
    1346:         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
    1347:         PetscCall(PetscFree(ftable));
    1348:         ftable = ftablenew;
    1349:       }
    1350:       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
    1351:       ftableused += countused;
    1352:       ftable[ftableused] = 0;
    1353:     }
    1354:     PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));

    1356:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
    1357:     PetscCall(PetscFree(ftable));
    1358:   }
    1359:   if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
    1360:   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
    1361:     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
    1362:       PetscCall(TSGetStepNumber(ts, &stepnum));
    1363:     }
    1364:     PetscCall(PetscSNPrintf(filename, sizeof filename, "ex11-%03" PetscInt_FMT ".vtu", stepnum));
    1365:     PetscCall(OutputVTK(dm, filename, &viewer));
    1366:     PetscCall(VecView(X, viewer));
    1367:     PetscCall(PetscViewerDestroy(&viewer));
    1368:   }
    1369:   PetscFunctionReturn(PETSC_SUCCESS);
    1370: }

    1372: static PetscErrorCode OutputBIN(DM dm, const char *filename, PetscViewer *viewer)
    1373: {
    1374:   PetscFunctionBeginUser;
    1375:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
    1376:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERBINARY));
    1377:   PetscCall(PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE));
    1378:   PetscCall(PetscViewerFileSetName(*viewer, filename));
    1379:   PetscFunctionReturn(PETSC_SUCCESS);
    1380: }

    1382: static PetscErrorCode TestMonitor(DM dm, const char *filename, Vec X, PetscReal time)
    1383: {
    1384:   Vec         odesolution;
    1385:   PetscInt    Nr;
    1386:   PetscBool   equal;
    1387:   PetscReal   timeread;
    1388:   PetscViewer viewer;

    1390:   PetscFunctionBeginUser;
    1391:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
    1392:   PetscCall(VecCreate(PETSC_COMM_WORLD, &odesolution));
    1393:   PetscCall(VecLoad(odesolution, viewer));
    1394:   VecEqual(X, odesolution, &equal);
    1395:   PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the vec data from file");
    1396:   else
    1397:   {
    1398:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for Vec\n"));
    1399:   }
    1400:   /*Nr   = 1;
    1401:    PetscCall(PetscRealLoad(Nr,&Nr,&timeread,viewer));*/
    1402:   PetscCall(PetscViewerBinaryRead(viewer, &timeread, 1, NULL, PETSC_REAL));

    1404:   PetscCheck(timeread == time, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the scalar data from file");
    1405:   else
    1406:   {
    1407:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for PetscReal\n"));
    1408:   }

    1410:   PetscCall(PetscViewerDestroy(&viewer));
    1411:   PetscCall(VecDestroy(&odesolution));
    1412:   PetscFunctionReturn(PETSC_SUCCESS);
    1413: }

    1415: static PetscErrorCode MonitorBIN(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
    1416: {
    1417:   User        user = (User)ctx;
    1418:   DM          dm;
    1419:   PetscViewer viewer;
    1420:   char        filename[PETSC_MAX_PATH_LEN];

    1422:   PetscFunctionBeginUser;
    1423:   PetscCall(VecGetDM(X, &dm));
    1424:   PetscCall(PetscSNPrintf(filename, sizeof filename, "ex11-SA-%06d.bin", stepnum));
    1425:   PetscCall(OutputBIN(dm, filename, &viewer));
    1426:   PetscCall(VecView(X, viewer));
    1427:   PetscCall(PetscRealView(1, &time, viewer));
    1428:   /* PetscCall(PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL));*/
    1429:   PetscCall(PetscViewerDestroy(&viewer));
    1430:   PetscCall(TestMonitor(dm, filename, X, time));
    1431:   PetscFunctionReturn(PETSC_SUCCESS);
    1432: }

    1434: int main(int argc, char **argv)
    1435: {
    1436:   MPI_Comm          comm;
    1437:   PetscFV           fvm;
    1438:   User              user;
    1439:   Model             mod;
    1440:   Physics           phys;
    1441:   DM                dm, plex;
    1442:   PetscReal         ftime, cfl, dt, minRadius;
    1443:   PetscInt          dim, nsteps;
    1444:   TS                ts;
    1445:   TSConvergedReason reason;
    1446:   Vec               X;
    1447:   PetscViewer       viewer;
    1448:   PetscMPIInt       rank;
    1449:   PetscBool         vtkCellGeom, splitFaces;
    1450:   PetscInt          overlap, f;
    1451:   char              filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";

    1453:   PetscFunctionBeginUser;
    1454:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    1455:   comm = PETSC_COMM_WORLD;
    1456:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    1458:   PetscCall(PetscNew(&user));
    1459:   PetscCall(PetscNew(&user->model));
    1460:   PetscCall(PetscNew(&user->model->physics));
    1461:   mod       = user->model;
    1462:   phys      = mod->physics;
    1463:   mod->comm = comm;

    1465:   /* Register physical models to be available on the command line */
    1466:   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
    1467:   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
    1468:   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));

    1470:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
    1471:   {
    1472:     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
    1473:     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
    1474:     PetscCall(PetscOptionsString("-f", "Exodus.II filename to read", "", filename, filename, sizeof(filename), NULL));
    1475:     splitFaces = PETSC_FALSE;
    1476:     PetscCall(PetscOptionsBool("-ufv_split_faces", "Split faces between cell sets", "", splitFaces, &splitFaces, NULL));
    1477:     overlap = 1;
    1478:     PetscCall(PetscOptionsInt("-ufv_mesh_overlap", "Number of cells to overlap partitions", "", overlap, &overlap, NULL));
    1479:     user->vtkInterval = 1;
    1480:     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
    1481:     vtkCellGeom = PETSC_FALSE;
    1482:     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
    1483:   }
    1484:   PetscOptionsEnd();
    1485:   PetscCall(DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, &dm));
    1486:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
    1487:   PetscCall(DMGetDimension(dm, &dim));

    1489:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
    1490:   {
    1491:     PetscDS prob;
    1492:     PetscErrorCode (*physcreate)(PetscDS, Model, Physics);
    1493:     char physname[256] = "advect";

    1495:     PetscCall(DMCreateLabel(dm, "Face Sets"));
    1496:     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
    1497:     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
    1498:     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
    1499:     PetscCall(DMGetDS(dm, &prob));
    1500:     PetscCall((*physcreate)(prob, mod, phys));
    1501:     mod->maxspeed = phys->maxspeed;
    1502:     /* Count number of fields and dofs */
    1503:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;

    1505:     PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
    1506:     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
    1507:     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
    1508:   }
    1509:   PetscOptionsEnd();
    1510:   {
    1511:     DM dmDist;

    1513:     PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
    1514:     PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist));
    1515:     if (dmDist) {
    1516:       PetscCall(DMDestroy(&dm));
    1517:       dm = dmDist;
    1518:     }
    1519:   }
    1520:   PetscCall(DMSetFromOptions(dm));
    1521:   {
    1522:     DM gdm;

    1524:     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
    1525:     PetscCall(DMDestroy(&dm));
    1526:     dm = gdm;
    1527:     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
    1528:   }
    1529:   if (splitFaces) PetscCall(ConstructCellBoundary(dm, user));
    1530:   PetscCall(SplitFaces(&dm, "split faces", user));

    1532:   PetscCall(PetscFVCreate(comm, &fvm));
    1533:   PetscCall(PetscFVSetFromOptions(fvm));
    1534:   PetscCall(DMSetNumFields(dm, phys->nfields));
    1535:   for (f = 0; f < phys->nfields; ++f) PetscCall(DMSetField(dm, f, (PetscObject)fvm));
    1536:   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
    1537:   PetscCall(PetscFVSetSpatialDimension(fvm, dim));

    1539:   /* Set up DM with section describing local vector and configure local vector. */
    1540:   PetscCall(SetUpLocalSpace(dm, user));

    1542:   PetscCall(TSCreate(comm, &ts));
    1543:   PetscCall(TSSetType(ts, TSSSP));
    1544:   /* PetscCall(TSSetType(ts, TSRK)); */
    1545:   PetscCall(TSSetDM(ts, dm));
    1546:   /* PetscCall(TSMonitorSet(ts,MonitorVTK,user,NULL)); */
    1547:   PetscCall(TSMonitorSet(ts, MonitorBIN, user, NULL));
    1548:   PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
    1549:   PetscCall(DMCreateGlobalVector(dm, &X));
    1550:   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
    1551:   PetscCall(SetInitialCondition(dm, X, user));
    1552:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1553:   if (vtkCellGeom) {
    1554:     DM  dmCell;
    1555:     Vec cellgeom, partition;

    1557:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
    1558:     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
    1559:     PetscCall(VecView(cellgeom, viewer));
    1560:     PetscCall(PetscViewerDestroy(&viewer));
    1561:     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
    1562:     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
    1563:     PetscCall(VecView(partition, viewer));
    1564:     PetscCall(PetscViewerDestroy(&viewer));
    1565:     PetscCall(VecDestroy(&partition));
    1566:     PetscCall(DMDestroy(&dmCell));
    1567:   }

    1569:   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
    1570:   PetscCall(DMDestroy(&plex));
    1571:   PetscCall(TSSetMaxTime(ts, 2.0));
    1572:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1573:   dt = cfl * minRadius / user->model->maxspeed;
    1574:   PetscCall(TSSetTimeStep(ts, dt));
    1575:   PetscCall(TSSetFromOptions(ts));
    1576:   PetscCall(TSSolve(ts, X));
    1577:   PetscCall(TSGetSolveTime(ts, &ftime));
    1578:   PetscCall(TSGetStepNumber(ts, &nsteps));
    1579:   PetscCall(TSGetConvergedReason(ts, &reason));
    1580:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
    1581:   PetscCall(TSDestroy(&ts));

    1583:   PetscCall(PetscFunctionListDestroy(&PhysicsList));
    1584:   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
    1585:   PetscCall(PetscFree(user->model->functionalMonitored));
    1586:   PetscCall(PetscFree(user->model->functionalCall));
    1587:   PetscCall(PetscFree(user->model->physics->data));
    1588:   PetscCall(PetscFree(user->model->physics));
    1589:   PetscCall(PetscFree(user->model));
    1590:   PetscCall(PetscFree(user));
    1591:   PetscCall(VecDestroy(&X));
    1592:   PetscCall(PetscFVDestroy(&fvm));
    1593:   PetscCall(DMDestroy(&dm));
    1594:   PetscCall(PetscFinalize());
    1595:   return 0;
    1596: }