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

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

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

A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
\begin{eqnarray}
  f^{L,R}_h    &=& uh^{L,R} \cdot \hat n \\
  f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
  c^{L,R}      &=& \sqrt{g h^{L,R}} \\
  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) \\
  f_i          &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
\end{eqnarray}
where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.

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

The mesh is read in from an ExodusII file, usually generated by Cubit.
F*/
#include <petscts.h>
#include <petscfv.h>
#include <petscdmplex.h>
#include <petscsf.h>
#include <petscblaslapack.h>

#define DIM 2 /* Geometric dimension */

static PetscFunctionList PhysicsList;

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

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

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

typedef PetscErrorCode (*RiemannFunction)(const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscScalar *, void *);
typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);

struct FieldDescription {
  const char *name;
  PetscInt    dof;
};

typedef struct _n_FunctionalLink *FunctionalLink;
struct _n_FunctionalLink {
  char              *name;
  FunctionalFunction func;
  void              *ctx;
  PetscInt           offset;
  FunctionalLink     next;
};

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

struct _n_Model {
  MPI_Comm         comm; /* Does not do collective communication, but some error conditions can be collective */
  Physics          physics;
  FunctionalLink   functionalRegistry;
  PetscInt         maxComputed;
  PetscInt         numMonitored;
  FunctionalLink  *functionalMonitored;
  PetscInt         numCall;
  FunctionalLink  *functionalCall;
  SolutionFunction solution;
  void            *solutionctx;
  PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
};

struct _n_User {
  PetscInt numSplitFaces;
  PetscInt vtkInterval; /* For monitor */
  Model    model;
};

static inline PetscScalar DotDIM(const PetscScalar *x, const PetscScalar *y)
{
  PetscInt    i;
  PetscScalar prod = 0.0;

  for (i = 0; i < DIM; i++) prod += x[i] * y[i];
  return prod;
}
static inline PetscReal NormDIM(const PetscScalar *x)
{
  return PetscSqrtReal(PetscAbsScalar(DotDIM(x, x)));
}
static inline void axDIM(const PetscScalar a, PetscScalar *x)
{
  PetscInt i;
  for (i = 0; i < DIM; i++) x[i] *= a;
}
static inline void waxDIM(const PetscScalar a, const PetscScalar *x, PetscScalar *w)
{
  PetscInt i;
  for (i = 0; i < DIM; i++) w[i] = x[i] * a;
}
static inline void NormalSplitDIM(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
{ /* Split x into normal and tangential components */
  PetscInt    i;
  PetscScalar c;
  c = DotDIM(x, n) / DotDIM(n, n);
  for (i = 0; i < DIM; i++) {
    xn[i] = c * n[i];
    xt[i] = x[i] - xn[i];
  }
}

static inline PetscScalar Dot2(const PetscScalar *x, const PetscScalar *y)
{
  return x[0] * y[0] + x[1] * y[1];
}
static inline PetscReal Norm2(const PetscScalar *x)
{
  return PetscSqrtReal(PetscAbsScalar(Dot2(x, x)));
}
static inline void Normalize2(PetscScalar *x)
{
  PetscReal a = 1. / Norm2(x);
  x[0] *= a;
  x[1] *= a;
}
static inline void Waxpy2(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
{
  w[0] = a * x[0] + y[0];
  w[1] = a * x[1] + y[1];
}
static inline void Scale2(PetscScalar a, const PetscScalar *x, PetscScalar *y)
{
  y[0] = a * x[0];
  y[1] = a * x[1];
}

static inline void WaxpyD(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
{
  PetscInt d;
  for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
}
static inline PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
{
  PetscScalar sum = 0.0;
  PetscInt    d;
  for (d = 0; d < dim; ++d) sum += x[d] * y[d];
  return sum;
}
static inline PetscReal NormD(PetscInt dim, const PetscScalar *x)
{
  return PetscSqrtReal(PetscAbsScalar(DotD(dim, x, x)));
}

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

/******************* Advect ********************/
typedef enum {
  ADVECT_SOL_TILTED,
  ADVECT_SOL_BUMP
} AdvectSolType;
static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "AdvectSolType", "ADVECT_SOL_", 0};
typedef enum {
  ADVECT_SOL_BUMP_CONE,
  ADVECT_SOL_BUMP_COS
} AdvectSolBumpType;
static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};

typedef struct {
  PetscReal wind[DIM];
} Physics_Advect_Tilted;
typedef struct {
  PetscReal         center[DIM];
  PetscReal         radius;
  AdvectSolBumpType type;
} Physics_Advect_Bump;

typedef struct {
  PetscReal     inflowState;
  AdvectSolType soltype;
  union
  {
    Physics_Advect_Tilted tilted;
    Physics_Advect_Bump   bump;
  } sol;
  struct {
    PetscInt Error;
  } functional;
} Physics_Advect;

static const struct FieldDescription PhysicsFields_Advect[] = {
  {"U",  1},
  {NULL, 0}
};

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

  PetscFunctionBeginUser;
  xG[0] = advect->inflowState;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
{
  PetscFunctionBeginUser;
  xG[0] = xI[0];
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  switch (advect->soltype) {
  case ADVECT_SOL_TILTED: {
    Physics_Advect_Tilted *tilted = &advect->sol.tilted;
    wind[0]                       = tilted->wind[0];
    wind[1]                       = tilted->wind[1];
  } break;
  case ADVECT_SOL_BUMP:
    wind[0] = -qp[1];
    wind[1] = qp[0];
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for solution type %s", AdvectSolBumpTypes[advect->soltype]);
  }
  wn      = Dot2(wind, n);
  flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  switch (advect->soltype) {
  case ADVECT_SOL_TILTED: {
    PetscReal              x0[DIM];
    Physics_Advect_Tilted *tilted = &advect->sol.tilted;
    Waxpy2(-time, tilted->wind, x, x0);
    if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
    else u[0] = advect->inflowState;
  } break;
  case ADVECT_SOL_BUMP: {
    Physics_Advect_Bump *bump = &advect->sol.bump;
    PetscReal            x0[DIM], v[DIM], r, cost, sint;
    cost  = PetscCosReal(time);
    sint  = PetscSinReal(time);
    x0[0] = cost * x[0] + sint * x[1];
    x0[1] = -sint * x[0] + cost * x[1];
    Waxpy2(-1, bump->center, x0, v);
    r = Norm2(v);
    switch (bump->type) {
    case ADVECT_SOL_BUMP_CONE:
      u[0] = PetscMax(1 - r / bump->radius, 0);
      break;
    case ADVECT_SOL_BUMP_COS:
      u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
      break;
    }
  } break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

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

  PetscFunctionBeginUser;
  phys->field_desc = PhysicsFields_Advect;
  phys->riemann    = (RiemannFunction)PhysicsRiemann_Advect;
  PetscCall(PetscNew(&advect));
  phys->data = advect;
  PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
  {
    PetscInt two = 2, dof = 1;
    advect->soltype = ADVECT_SOL_TILTED;
    PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
    switch (advect->soltype) {
    case ADVECT_SOL_TILTED: {
      Physics_Advect_Tilted *tilted = &advect->sol.tilted;
      two                           = 2;
      tilted->wind[0]               = 0.0;
      tilted->wind[1]               = 1.0;
      PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
      advect->inflowState = -2.0;
      PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
      phys->maxspeed = Norm2(tilted->wind);
    } break;
    case ADVECT_SOL_BUMP: {
      Physics_Advect_Bump *bump = &advect->sol.bump;
      two                       = 2;
      bump->center[0]           = 2.;
      bump->center[1]           = 0.;
      PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
      bump->radius = 0.9;
      PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
      bump->type = ADVECT_SOL_BUMP_CONE;
      PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
      phys->maxspeed = 3.; /* radius of mesh, kludge */
    } break;
    }
  }
  PetscOptionsHeadEnd();
  {
    const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
    DMLabel        label;

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

/******************* Shallow Water ********************/
typedef struct {
  PetscReal gravity;
  PetscReal boundaryHeight;
  struct {
    PetscInt Height;
    PetscInt Speed;
    PetscInt Energy;
  } functional;
} Physics_SW;
typedef struct {
  PetscScalar vals[0];
  PetscScalar h;
  PetscScalar uh[DIM];
} SWNode;

static const struct FieldDescription PhysicsFields_SW[] = {
  {"Height",   1  },
  {"Momentum", DIM},
  {NULL,       0  }
};

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

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

static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
{
  PetscFunctionBeginUser;
  xG[0] = xI[0];
  xG[1] = -xI[1];
  xG[2] = -xI[2];
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCheck(uL->h >= 0 && uR->h >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
  nn[0] = n[0];
  nn[1] = n[1];
  Normalize2(nn);
  SWFlux(phys, nn, uL, &fL);
  SWFlux(phys, nn, uR, &fR);
  cL    = PetscSqrtReal(sw->gravity * PetscRealPart(uL->h));
  cR    = PetscSqrtReal(sw->gravity * PetscRealPart(uR->h)); /* gravity wave speed */
  speed = PetscMax(PetscAbsScalar(Dot2(uL->uh, nn) / uL->h) + cL, PetscAbsScalar(Dot2(uR->uh, nn) / uR->h) + cR);
  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);
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
  dx[0] = x[0] - 1.5;
  dx[1] = x[1] - 1.0;
  r     = Norm2(dx);
  sigma = 0.5;
  u[0]  = 1 + 2 * PetscExpScalar(-PetscSqr(r) / (2 * PetscSqr(sigma)));
  u[1]  = 0.0;
  u[2]  = 0.0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

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

  PetscFunctionBeginUser;
  phys->field_desc = PhysicsFields_SW;
  phys->riemann    = (RiemannFunction)PhysicsRiemann_SW;
  PetscCall(PetscNew(&sw));
  phys->data = sw;
  PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
  {
    sw->gravity = 1.0;
    PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
  }
  PetscOptionsHeadEnd();
  phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */

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

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

/******************* Euler ********************/
typedef struct {
  PetscScalar vals[0];
  PetscScalar r;
  PetscScalar ru[DIM];
  PetscScalar e;
} EulerNode;
typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode *, PetscScalar *);
typedef struct {
  PetscInt        npars;
  PetscReal       pars[DIM];
  EquationOfState pressure;
  EquationOfState sound;
  struct {
    PetscInt Density;
    PetscInt Momentum;
    PetscInt Energy;
    PetscInt Pressure;
    PetscInt Speed;
  } monitor;
} Physics_Euler;

static const struct FieldDescription PhysicsFields_Euler[] = {
  {"Density",  1  },
  {"Momentum", DIM},
  {"Energy",   1  },
  {NULL,       0  }
};

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

  PetscFunctionBeginUser;
  ru2 = DotDIM(x->ru, x->ru);
  ru2 /= x->r;
  /* kinematic dof = params[0] */
  (*p) = 2.0 * (x->e - 0.5 * ru2) / pars[0];
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  /* TODO remove direct usage of Pressure_PG */
  Pressure_PG(pars, x, &p);
  /* TODO check the sign of p */
  /* pars[1] = heat capacity ratio */
  (*c) = PetscSqrtScalar(pars[1] * p / x->r);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
 * x = (rho,rho*(u_1),...,rho*e)^T
 * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
 *
 * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
 *
 * */
static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
{
  Physics_Euler *eu = (Physics_Euler *)phys->data;
  PetscScalar    u, nu, p;
  PetscInt       i;

  PetscFunctionBeginUser;
  u = DotDIM(x->ru, x->ru);
  u /= (x->r * x->r);
  nu = DotDIM(x->ru, n);
  /* TODO check the sign of p */
  eu->pressure(eu->pars, x, &p);
  f->r = nu * x->r;
  for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p;
  f->e = nu * (x->e + p);
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  xG[0] = xI[0];
  NormalSplitDIM(n, xI + 1, xn, xt);
  for (i = 0; i < DIM; i++) xG[i + 1] = -xn[i] + xt[i];
  xG[DIM + 1] = xI[DIM + 1];
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

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

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

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

  PetscFunctionBeginUser;
  f[eu->monitor.Density]  = x->r;
  f[eu->monitor.Momentum] = NormDIM(x->ru);
  f[eu->monitor.Energy]   = x->e;
  f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
  eu->pressure(eu->pars, x, &p);
  f[eu->monitor.Pressure] = p;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  phys->field_desc = PhysicsFields_Euler;
  phys->riemann    = (RiemannFunction)PhysicsRiemann_Euler_Rusanov;
  PetscCall(PetscNew(&eu));
  phys->data = eu;
  PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
  {
    eu->pars[0] = 3.0;
    eu->pars[1] = 1.67;
    PetscCall(PetscOptionsReal("-eu_f", "Degrees of freedom", "", eu->pars[0], &eu->pars[0], NULL));
    PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[1], &eu->pars[1], NULL));
  }
  PetscOptionsHeadEnd();
  eu->pressure   = Pressure_PG;
  eu->sound      = SpeedOfSound_PG;
  phys->maxspeed = 1.0;
  {
    const PetscInt wallids[] = {100, 101, 200, 300};
    DMLabel        label;

    PetscCall(DMGetLabel(dm, "Face Sets", &label));
    PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
    PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
    PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
    PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
    PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
    PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
    PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode ConstructCellBoundary(DM dm, User user)
{
  const char     *name   = "Cell Sets";
  const char     *bdname = "split faces";
  IS              regionIS, innerIS;
  const PetscInt *regions, *cells;
  PetscInt        numRegions, innerRegion, numCells, c;
  PetscInt        cStart, cEnd, cEndInterior, fStart, fEnd;
  PetscBool       hasLabel;

  PetscFunctionBeginUser;
  PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
  PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
  PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));

  PetscCall(DMHasLabel(dm, name, &hasLabel));
  if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(DMGetLabelSize(dm, name, &numRegions));
  if (numRegions != 2) PetscFunctionReturn(PETSC_SUCCESS);
  /* Get the inner id */
  PetscCall(DMGetLabelIdIS(dm, name, &regionIS));
  PetscCall(ISGetIndices(regionIS, &regions));
  innerRegion = regions[0];
  PetscCall(ISRestoreIndices(regionIS, &regions));
  PetscCall(ISDestroy(&regionIS));
  /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
  PetscCall(DMGetStratumIS(dm, name, innerRegion, &innerIS));
  PetscCall(ISGetLocalSize(innerIS, &numCells));
  PetscCall(ISGetIndices(innerIS, &cells));
  PetscCall(DMCreateLabel(dm, bdname));
  for (c = 0; c < numCells; ++c) {
    const PetscInt  cell = cells[c];
    const PetscInt *faces;
    PetscInt        numFaces, f;

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

      PetscCheck((face >= fStart) && (face < fEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a face", face);
      PetscCall(DMPlexGetSupportSize(dm, face, &nC));
      if (nC != 2) continue;
      PetscCall(DMPlexGetSupport(dm, face, &neighbors));
      if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue;
      PetscCheck((neighbors[0] >= cStart) && (neighbors[0] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[0]);
      PetscCheck((neighbors[1] >= cStart) && (neighbors[1] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[1]);
      PetscCall(DMGetLabelValue(dm, name, neighbors[0], &regionA));
      PetscCall(DMGetLabelValue(dm, name, neighbors[1], &regionB));
      PetscCheck(regionA >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[0]);
      PetscCheck(regionB >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[1]);
      if (regionA != regionB) PetscCall(DMSetLabelValue(dm, bdname, faces[f], 1));
    }
  }
  PetscCall(ISRestoreIndices(innerIS, &cells));
  PetscCall(ISDestroy(&innerIS));
  {
    DMLabel label;

    PetscCall(DMGetLabel(dm, bdname, &label));
    PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCall(DMHasLabel(dm, labelName, &hasLabel));
  if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm));
  PetscCall(DMSetType(sdm, DMPLEX));
  PetscCall(DMGetDimension(dm, &dim));
  PetscCall(DMSetDimension(sdm, dim));

  PetscCall(DMGetLabelIdIS(dm, labelName, &idIS));
  PetscCall(ISGetLocalSize(idIS, &numFS));
  PetscCall(ISGetIndices(idIS, &ids));

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

    PetscCall(DMGetStratumSize(dm, labelName, ids[fs], &numBdFaces));
    user->numSplitFaces += numBdFaces;
  }
  PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
  pEnd += user->numSplitFaces;
  PetscCall(DMPlexSetChart(sdm, pStart, pEnd));
  PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
  PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
  numGhostCells = cEnd - cEndInterior;
  /* Set cone and support sizes */
  PetscCall(DMPlexGetDepth(dm, &depth));
  for (d = 0; d <= depth; ++d) {
    PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    for (p = pStart; p < pEnd; ++p) {
      PetscInt newp = p;
      PetscInt size;

      PetscCall(DMPlexGetConeSize(dm, p, &size));
      PetscCall(DMPlexSetConeSize(sdm, newp, size));
      PetscCall(DMPlexGetSupportSize(dm, p, &size));
      PetscCall(DMPlexSetSupportSize(sdm, newp, size));
    }
  }
  PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
  for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
    IS              faceIS;
    const PetscInt *faces;
    PetscInt        numFaces, f;

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

      /* Right now I think that both faces should see both cells */
      PetscCall(DMPlexGetConeSize(dm, faces[f], &size));
      PetscCall(DMPlexSetConeSize(sdm, newf, size));
      PetscCall(DMPlexGetSupportSize(dm, faces[f], &size));
      PetscCall(DMPlexSetSupportSize(sdm, newf, size));
    }
    PetscCall(ISRestoreIndices(faceIS, &faces));
    PetscCall(ISDestroy(&faceIS));
  }
  PetscCall(DMSetUp(sdm));
  /* Set cones and supports */
  PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
  PetscCall(PetscMalloc1(PetscMax(maxConeSize, maxSupportSize), &newpoints));
  PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
  for (p = pStart; p < pEnd; ++p) {
    const PetscInt *points, *orientations;
    PetscInt        size, i, newp = p;

    PetscCall(DMPlexGetConeSize(dm, p, &size));
    PetscCall(DMPlexGetCone(dm, p, &points));
    PetscCall(DMPlexGetConeOrientation(dm, p, &orientations));
    for (i = 0; i < size; ++i) newpoints[i] = points[i];
    PetscCall(DMPlexSetCone(sdm, newp, newpoints));
    PetscCall(DMPlexSetConeOrientation(sdm, newp, orientations));
    PetscCall(DMPlexGetSupportSize(dm, p, &size));
    PetscCall(DMPlexGetSupport(dm, p, &points));
    for (i = 0; i < size; ++i) newpoints[i] = points[i];
    PetscCall(DMPlexSetSupport(sdm, newp, newpoints));
  }
  PetscCall(PetscFree(newpoints));
  for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
    IS              faceIS;
    const PetscInt *faces;
    PetscInt        numFaces, f;

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

      PetscCall(DMPlexGetCone(dm, faces[f], &points));
      PetscCall(DMPlexSetCone(sdm, newf, points));
      PetscCall(DMPlexGetSupport(dm, faces[f], &points));
      PetscCall(DMPlexSetSupport(sdm, newf, points));
    }
    PetscCall(ISRestoreIndices(faceIS, &faces));
    PetscCall(ISDestroy(&faceIS));
  }
  PetscCall(ISRestoreIndices(idIS, &ids));
  PetscCall(ISDestroy(&idIS));
  PetscCall(DMPlexStratify(sdm));
  /* Convert coordinates */
  PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
  PetscCall(DMGetCoordinateSection(dm, &coordSection));
  PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection));
  PetscCall(PetscSectionSetNumFields(newCoordSection, 1));
  PetscCall(PetscSectionSetFieldComponents(newCoordSection, 0, dim));
  PetscCall(PetscSectionSetChart(newCoordSection, vStart, vEnd));
  for (v = vStart; v < vEnd; ++v) {
    PetscCall(PetscSectionSetDof(newCoordSection, v, dim));
    PetscCall(PetscSectionSetFieldDof(newCoordSection, v, 0, dim));
  }
  PetscCall(PetscSectionSetUp(newCoordSection));
  PetscCall(DMSetCoordinateSection(sdm, PETSC_DETERMINE, newCoordSection));
  PetscCall(PetscSectionDestroy(&newCoordSection)); /* relinquish our reference */
  PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
  PetscCall(DMSetCoordinatesLocal(sdm, coordinates));
  /* Convert labels */
  PetscCall(DMGetNumLabels(dm, &numLabels));
  for (l = 0; l < numLabels; ++l) {
    const char *lname;
    PetscBool   isDepth;

    PetscCall(DMGetLabelName(dm, l, &lname));
    PetscCall(PetscStrcmp(lname, "depth", &isDepth));
    if (isDepth) continue;
    PetscCall(DMCreateLabel(sdm, lname));
    PetscCall(DMGetLabelIdIS(dm, lname, &idIS));
    PetscCall(ISGetLocalSize(idIS, &numFS));
    PetscCall(ISGetIndices(idIS, &ids));
    for (fs = 0; fs < numFS; ++fs) {
      IS              pointIS;
      const PetscInt *points;
      PetscInt        numPoints;

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

        PetscCall(DMSetLabelValue(sdm, lname, newpoint, ids[fs]));
      }
      PetscCall(ISRestoreIndices(pointIS, &points));
      PetscCall(ISDestroy(&pointIS));
    }
    PetscCall(ISRestoreIndices(idIS, &ids));
    PetscCall(ISDestroy(&idIS));
  }
  /* Convert pointSF */
  const PetscSFNode *remotePoints;
  PetscSFNode       *gremotePoints;
  const PetscInt    *localPoints;
  PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
  PetscInt           numRoots, numLeaves;
  PetscMPIInt        size;

  PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
  PetscCall(DMGetPointSF(dm, &sfPoint));
  PetscCall(DMGetPointSF(sdm, &gsfPoint));
  PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
  PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
  if (numRoots >= 0) {
    PetscCall(PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation));
    for (l = 0; l < numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? numGhostCells : 0); */
    PetscCall(PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
    PetscCall(PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
    PetscCall(PetscMalloc1(numLeaves, &glocalPoints));
    PetscCall(PetscMalloc1(numLeaves, &gremotePoints));
    for (l = 0; l < numLeaves; ++l) {
      glocalPoints[l]        = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + numGhostCells : localPoints[l]; */
      gremotePoints[l].rank  = remotePoints[l].rank;
      gremotePoints[l].index = newRemoteLocation[localPoints[l]];
    }
    PetscCall(PetscFree2(newLocation, newRemoteLocation));
    PetscCall(PetscSFSetGraph(gsfPoint, numRoots + numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER));
  }
  PetscCall(DMDestroy(dmSplit));
  *dmSplit = sdm;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
{
  PetscSF      sfPoint;
  PetscSection coordSection;
  Vec          coordinates;
  PetscSection sectionCell;
  PetscScalar *part;
  PetscInt     cStart, cEnd, c;
  PetscMPIInt  rank;

  PetscFunctionBeginUser;
  PetscCall(DMGetCoordinateSection(dm, &coordSection));
  PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
  PetscCall(DMClone(dm, dmCell));
  PetscCall(DMGetPointSF(dm, &sfPoint));
  PetscCall(DMSetPointSF(*dmCell, sfPoint));
  PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
  PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
  PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
  PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
  PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
  PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
  for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
  PetscCall(PetscSectionSetUp(sectionCell));
  PetscCall(DMSetLocalSection(*dmCell, sectionCell));
  PetscCall(PetscSectionDestroy(&sectionCell));
  PetscCall(DMCreateLocalVector(*dmCell, partition));
  PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
  PetscCall(VecGetArray(*partition, &part));
  for (c = cStart; c < cEnd; ++c) {
    PetscScalar *p;

    PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
    p[0] = rank;
  }
  PetscCall(VecRestoreArray(*partition, &part));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
{
  DM                 plex, dmMass, dmFace, dmCell, dmCoord;
  PetscSection       coordSection;
  Vec                coordinates, facegeom, cellgeom;
  PetscSection       sectionMass;
  PetscScalar       *m;
  const PetscScalar *fgeom, *cgeom, *coords;
  PetscInt           vStart, vEnd, v;

  PetscFunctionBeginUser;
  PetscCall(DMConvert(dm, DMPLEX, &plex));
  PetscCall(DMGetCoordinateSection(dm, &coordSection));
  PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
  PetscCall(DMClone(dm, &dmMass));
  PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
  PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
  PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
  PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
  PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
  for (v = vStart; v < vEnd; ++v) {
    PetscInt numFaces;

    PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
    PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
  }
  PetscCall(PetscSectionSetUp(sectionMass));
  PetscCall(DMSetLocalSection(dmMass, sectionMass));
  PetscCall(PetscSectionDestroy(&sectionMass));
  PetscCall(DMGetLocalVector(dmMass, massMatrix));
  PetscCall(VecGetArray(*massMatrix, &m));
  PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
  PetscCall(VecGetDM(facegeom, &dmFace));
  PetscCall(VecGetArrayRead(facegeom, &fgeom));
  PetscCall(VecGetDM(cellgeom, &dmCell));
  PetscCall(VecGetArrayRead(cellgeom, &cgeom));
  PetscCall(DMGetCoordinateDM(dm, &dmCoord));
  PetscCall(VecGetArrayRead(coordinates, &coords));
  for (v = vStart; v < vEnd; ++v) {
    const PetscInt        *faces;
    const PetscFVFaceGeom *fgA, *fgB, *cg;
    const PetscScalar     *vertex;
    PetscInt               numFaces, sides[2], f, g;

    PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
    PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
    PetscCall(DMPlexGetSupport(dmMass, v, &faces));
    for (f = 0; f < numFaces; ++f) {
      sides[0] = faces[f];
      PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
      for (g = 0; g < numFaces; ++g) {
        const PetscInt *cells = NULL;
        PetscReal       area  = 0.0;
        PetscInt        numCells;

        sides[1] = faces[g];
        PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
        PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
        PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
        PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
        area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
        area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
        m[f * numFaces + g] = Dot2(fgA->normal, fgB->normal) * area * 0.5;
        PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
      }
    }
  }
  PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
  PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
  PetscCall(VecRestoreArrayRead(coordinates, &coords));
  PetscCall(VecRestoreArray(*massMatrix, &m));
  PetscCall(DMDestroy(&dmMass));
  PetscCall(DMDestroy(&plex));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
  PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
  PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection));
  phys = user->model->physics;
  PetscCall(PetscSectionSetNumFields(stateSection, phys->nfields));
  for (i = 0; i < phys->nfields; i++) {
    PetscCall(PetscSectionSetFieldName(stateSection, i, phys->field_desc[i].name));
    PetscCall(PetscSectionSetFieldComponents(stateSection, i, phys->field_desc[i].dof));
  }
  PetscCall(PetscSectionSetChart(stateSection, cStart, cEnd));
  for (c = cStart; c < cEnd; ++c) {
    for (i = 0; i < phys->nfields; i++) PetscCall(PetscSectionSetFieldDof(stateSection, c, i, phys->field_desc[i].dof));
    PetscCall(PetscSectionSetDof(stateSection, c, dof));
  }
  for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintDof(stateSection, c, dof));
  PetscCall(PetscSectionSetUp(stateSection));
  PetscCall(PetscMalloc1(dof, &cind));
  for (d = 0; d < dof; ++d) cind[d] = d;
#if 0
  for (c = cStart; c < cEnd; ++c) {
    PetscInt val;

    PetscCall(DMGetLabelValue(dm, "vtk", c, &val));
    if (val < 0) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
  }
#endif
  for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
  PetscCall(PetscFree(cind));
  PetscCall(PetscSectionGetStorageSize(stateSection, &stateSize));
  PetscCall(DMSetLocalSection(dm, stateSection));
  PetscCall(PetscSectionDestroy(&stateSection));
  PetscFunctionReturn(PETSC_SUCCESS);
}

#if 0
PetscErrorCode SetUpBoundaries(DM dm, User user)
{
  Model          mod = user->model;
  BoundaryLink   b;

  PetscFunctionBeginUser;
  PetscOptionsBegin(PetscObjectComm((PetscObject)dm),NULL,"Boundary condition options","");
  for (b = mod->boundary; b; b=b->next) {
    char      optname[512];
    PetscInt  ids[512],len = 512;
    PetscBool flg;
    PetscCall(PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name));
    PetscCall(PetscMemzero(ids,sizeof(ids)));
    PetscCall(PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg));
    if (flg) {
      /* TODO: check all IDs to make sure they exist in the mesh */
      PetscCall(PetscFree(b->ids));
      b->numids = len;
      PetscCall(PetscMalloc1(len,&b->ids));
      PetscCall(PetscArraycpy(b->ids,ids,len));
    }
  }
  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}
#endif

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

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

  PetscFunctionBeginUser;
  for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
  PetscCall(PetscNew(&link));
  PetscCall(PetscStrallocpy(name, &link->name));
  link->offset = lastoffset + 1;
  link->func   = func;
  link->ctx    = ctx;
  link->next   = NULL;
  *ptr         = link;
  *offset      = link->offset;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptions *PetscOptionsObject)
{
  PetscInt       i, j;
  FunctionalLink link;
  char          *names[256];

  PetscFunctionBeginUser;
  mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
  PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
  /* Create list of functionals that will be computed somehow */
  PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
  /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
  PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
  mod->numCall = 0;
  for (i = 0; i < mod->numMonitored; i++) {
    for (link = mod->functionalRegistry; link; link = link->next) {
      PetscBool match;
      PetscCall(PetscStrcasecmp(names[i], link->name, &match));
      if (match) break;
    }
    PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
    mod->functionalMonitored[i] = link;
    for (j = 0; j < i; j++) {
      if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
    }
    mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
  next_name:
    PetscCall(PetscFree(names[i]));
  }

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

static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
{
  FunctionalLink l, next;

  PetscFunctionBeginUser;
  if (!link) PetscFunctionReturn(PETSC_SUCCESS);
  l     = *link;
  *link = NULL;
  for (; l; l = next) {
    next = l->next;
    PetscCall(PetscFree(l->name));
    PetscCall(PetscFree(l));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
{
  DM                 plex, dmCell;
  Model              mod = user->model;
  Vec                cellgeom;
  const PetscScalar *cgeom;
  PetscScalar       *x;
  PetscInt           cStart, cEnd, c;

  PetscFunctionBeginUser;
  PetscCall(DMConvert(dm, DMPLEX, &plex));
  PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
  PetscCall(VecGetDM(cellgeom, &dmCell));
  PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
  PetscCall(VecGetArrayRead(cellgeom, &cgeom));
  PetscCall(VecGetArray(X, &x));
  for (c = cStart; c < cEnd; ++c) {
    const PetscFVCellGeom *cg;
    PetscScalar           *xc;

    PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
    PetscCall(DMPlexPointGlobalRef(dm, c, x, &xc));
    if (xc) PetscCall((*mod->solution)(mod, 0.0, cg->centroid, xc, mod->solutionctx));
  }
  PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
  PetscCall(VecRestoreArray(X, &x));
  PetscCall(DMDestroy(&plex));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
{
  PetscFunctionBeginUser;
  PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
  PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
  PetscCall(PetscViewerFileSetName(*viewer, filename));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

    PetscCall(DMConvert(dm, DMPLEX, &plex));
    PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
    fcount = mod->maxComputed + 1;
    PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
    for (i = 0; i < fcount; i++) {
      fmin[i]      = PETSC_MAX_REAL;
      fmax[i]      = PETSC_MIN_REAL;
      fintegral[i] = 0;
    }
    PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    PetscCall(VecGetDM(cellgeom, &dmCell));
    PetscCall(VecGetArrayRead(cellgeom, &cgeom));
    PetscCall(VecGetArrayRead(X, &x));
    for (c = cStart; c < cEnd; ++c) {
      const PetscFVCellGeom *cg;
      const PetscScalar     *cx;
      PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
      PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
      if (!cx) continue; /* not a global cell */
      for (i = 0; i < mod->numCall; i++) {
        FunctionalLink flink = mod->functionalCall[i];
        PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
      }
      for (i = 0; i < fcount; i++) {
        fmin[i] = PetscMin(fmin[i], ftmp[i]);
        fmax[i] = PetscMax(fmax[i], ftmp[i]);
        fintegral[i] += cg->volume * ftmp[i];
      }
    }
    PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
    PetscCall(VecRestoreArrayRead(X, &x));
    PetscCall(DMDestroy(&plex));
    PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
    PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
    PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));

    ftablealloc = fcount * 100;
    ftableused  = 0;
    PetscCall(PetscMalloc1(ftablealloc, &ftable));
    for (i = 0; i < mod->numMonitored; i++) {
      size_t         countused;
      char           buffer[256], *p;
      FunctionalLink flink = mod->functionalMonitored[i];
      PetscInt       id    = flink->offset;
      if (i % 3) {
        PetscCall(PetscArraycpy(buffer, "  ", 2));
        p = buffer + 2;
      } else if (i) {
        char newline[] = "\n";
        PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
        p = buffer + sizeof(newline) - 1;
      } else {
        p = buffer;
      }
      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]));
      countused += p - buffer;
      if (countused > ftablealloc - ftableused - 1) { /* reallocate */
        char *ftablenew;
        ftablealloc = 2 * ftablealloc + countused;
        PetscCall(PetscMalloc(ftablealloc, &ftablenew));
        PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
        PetscCall(PetscFree(ftable));
        ftable = ftablenew;
      }
      PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
      ftableused += countused;
      ftable[ftableused] = 0;
    }
    PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));

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

static PetscErrorCode OutputBIN(DM dm, const char *filename, PetscViewer *viewer)
{
  PetscFunctionBeginUser;
  PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
  PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERBINARY));
  PetscCall(PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE));
  PetscCall(PetscViewerFileSetName(*viewer, filename));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TestMonitor(DM dm, const char *filename, Vec X, PetscReal time)
{
  Vec         odesolution;
  PetscInt    Nr;
  PetscBool   equal;
  PetscReal   timeread;
  PetscViewer viewer;

  PetscFunctionBeginUser;
  PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
  PetscCall(VecCreate(PETSC_COMM_WORLD, &odesolution));
  PetscCall(VecLoad(odesolution, viewer));
  VecEqual(X, odesolution, &equal);
  PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the vec data from file");
  else
  {
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for Vec\n"));
  }
  /*Nr   = 1;
   PetscCall(PetscRealLoad(Nr,&Nr,&timeread,viewer));*/
  PetscCall(PetscViewerBinaryRead(viewer, &timeread, 1, NULL, PETSC_REAL));

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

  PetscCall(PetscViewerDestroy(&viewer));
  PetscCall(VecDestroy(&odesolution));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCall(VecGetDM(X, &dm));
  PetscCall(PetscSNPrintf(filename, sizeof filename, "ex11-SA-%06d.bin", stepnum));
  PetscCall(OutputBIN(dm, filename, &viewer));
  PetscCall(VecView(X, viewer));
  PetscCall(PetscRealView(1, &time, viewer));
  /* PetscCall(PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL));*/
  PetscCall(PetscViewerDestroy(&viewer));
  PetscCall(TestMonitor(dm, filename, X, time));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  MPI_Comm          comm;
  PetscFV           fvm;
  User              user;
  Model             mod;
  Physics           phys;
  DM                dm, plex;
  PetscReal         ftime, cfl, dt, minRadius;
  PetscInt          dim, nsteps;
  TS                ts;
  TSConvergedReason reason;
  Vec               X;
  PetscViewer       viewer;
  PetscMPIInt       rank;
  PetscBool         vtkCellGeom, splitFaces;
  PetscInt          overlap, f;
  char              filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  comm = PETSC_COMM_WORLD;
  PetscCallMPI(MPI_Comm_rank(comm, &rank));

  PetscCall(PetscNew(&user));
  PetscCall(PetscNew(&user->model));
  PetscCall(PetscNew(&user->model->physics));
  mod       = user->model;
  phys      = mod->physics;
  mod->comm = comm;

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

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

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

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

    PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
    PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
    PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
  }
  PetscOptionsEnd();
  {
    DM dmDist;

    PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
    PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist));
    if (dmDist) {
      PetscCall(DMDestroy(&dm));
      dm = dmDist;
    }
  }
  PetscCall(DMSetFromOptions(dm));
  {
    DM gdm;

    PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
    PetscCall(DMDestroy(&dm));
    dm = gdm;
    PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
  }
  if (splitFaces) PetscCall(ConstructCellBoundary(dm, user));
  PetscCall(SplitFaces(&dm, "split faces", user));

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

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

  PetscCall(TSCreate(comm, &ts));
  PetscCall(TSSetType(ts, TSSSP));
  /* PetscCall(TSSetType(ts, TSRK)); */
  PetscCall(TSSetDM(ts, dm));
  /* PetscCall(TSMonitorSet(ts,MonitorVTK,user,NULL)); */
  PetscCall(TSMonitorSet(ts, MonitorBIN, user, NULL));
  PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
  PetscCall(DMCreateGlobalVector(dm, &X));
  PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
  PetscCall(SetInitialCondition(dm, X, user));
  PetscCall(DMConvert(dm, DMPLEX, &plex));
  if (vtkCellGeom) {
    DM  dmCell;
    Vec cellgeom, partition;

    PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
    PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
    PetscCall(VecView(cellgeom, viewer));
    PetscCall(PetscViewerDestroy(&viewer));
    PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
    PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
    PetscCall(VecView(partition, viewer));
    PetscCall(PetscViewerDestroy(&viewer));
    PetscCall(VecDestroy(&partition));
    PetscCall(DMDestroy(&dmCell));
  }

  PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
  PetscCall(DMDestroy(&plex));
  PetscCall(TSSetMaxTime(ts, 2.0));
  PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
  dt = cfl * minRadius / user->model->maxspeed;
  PetscCall(TSSetTimeStep(ts, dt));
  PetscCall(TSSetFromOptions(ts));
  PetscCall(TSSolve(ts, X));
  PetscCall(TSGetSolveTime(ts, &ftime));
  PetscCall(TSGetStepNumber(ts, &nsteps));
  PetscCall(TSGetConvergedReason(ts, &reason));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
  PetscCall(TSDestroy(&ts));

  PetscCall(PetscFunctionListDestroy(&PhysicsList));
  PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
  PetscCall(PetscFree(user->model->functionalMonitored));
  PetscCall(PetscFree(user->model->functionalCall));
  PetscCall(PetscFree(user->model->physics->data));
  PetscCall(PetscFree(user->model->physics));
  PetscCall(PetscFree(user->model));
  PetscCall(PetscFree(user));
  PetscCall(VecDestroy(&X));
  PetscCall(PetscFVDestroy(&fvm));
  PetscCall(DMDestroy(&dm));
  PetscCall(PetscFinalize());
  return 0;
}