static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter with semidiscrete time stepping.\n"
                           "  advection   - Constant coefficient scalar advection\n"
                           "                u_t       + (a*u)_x               = 0\n"
                           "  for this toy problem, we choose different meshsizes for different sub-domains, say\n"
                           "                hxs  = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n"
                           "                hxf  = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n"
                           "  with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of coarse\n\n"
                           "  grids and fine grids is hratio.\n"
                           "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
                           "                the states across shocks and rarefactions\n"
                           "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
                           "                also the reference solution should be generated by user and stored in a binary file.\n"
                           "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
                           "Several initial conditions can be chosen with -initial N\n\n"
                           "The problem size should be set with -da_grid_x M\n\n"
                           "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n"
                           "                             u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t))                 \n"
                           "                     limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2))))    \n"
                           "                             r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1)))                                   \n"
                           "                             alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1))                     \n"
                           "                             gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1))                 \n";

#include <petscts.h>
#include <petscdm.h>
#include <petscdmda.h>
#include <petscdraw.h>
#include <petscmath.h>

static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
{
  PetscReal range = xmax - xmin;
  return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
}

/* --------------------------------- Finite Volume data structures ----------------------------------- */

typedef enum {
  FVBC_PERIODIC,
  FVBC_OUTFLOW
} FVBCType;
static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", 0};

typedef struct {
  PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *);
  PetscErrorCode (*flux)(void *, const PetscScalar *, PetscScalar *, PetscReal *);
  PetscErrorCode (*destroy)(void *);
  void    *user;
  PetscInt dof;
  char    *fieldname[16];
} PhysicsCtx;

typedef struct {
  PhysicsCtx physics;
  MPI_Comm   comm;
  char       prefix[256];

  /* Local work arrays */
  PetscScalar *flux;   /* Flux across interface                                                      */
  PetscReal   *speeds; /* Speeds of each wave                                                        */
  PetscScalar *u;      /* value at face                                                              */

  PetscReal cfl_idt; /* Max allowable value of 1/Delta t                                           */
  PetscReal cfl;
  PetscReal xmin, xmax;
  PetscInt  initial;
  PetscBool exact;
  PetscBool simulation;
  FVBCType  bctype;
  PetscInt  hratio; /* hratio = hslow/hfast */
  IS        isf, iss;
  PetscInt  sf, fs; /* slow-fast and fast-slow interfaces */
} FVCtx;

/* --------------------------------- Physics ----------------------------------- */
static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
{
  PetscFunctionBeginUser;
  PetscCall(PetscFree(vctx));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* --------------------------------- Advection ----------------------------------- */
typedef struct {
  PetscReal a; /* advective velocity */
} AdvectCtx;

static PetscErrorCode PhysicsFlux_Advect(void *vctx, const PetscScalar *u, PetscScalar *flux, PetscReal *maxspeed)
{
  AdvectCtx *ctx = (AdvectCtx *)vctx;
  PetscReal  speed;

  PetscFunctionBeginUser;
  speed     = ctx->a;
  flux[0]   = speed * u[0];
  *maxspeed = speed;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
{
  AdvectCtx *ctx = (AdvectCtx *)vctx;
  PetscReal  a   = ctx->a, x0;

  PetscFunctionBeginUser;
  switch (bctype) {
  case FVBC_OUTFLOW:
    x0 = x - a * t;
    break;
  case FVBC_PERIODIC:
    x0 = RangeMod(x - a * t, xmin, xmax);
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
  }
  switch (initial) {
  case 0:
    u[0] = (x0 < 0) ? 1 : -1;
    break;
  case 1:
    u[0] = (x0 < 0) ? -1 : 1;
    break;
  case 2:
    u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
    break;
  case 3:
    u[0] = PetscSinReal(2 * PETSC_PI * x0);
    break;
  case 4:
    u[0] = PetscAbs(x0);
    break;
  case 5:
    u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    break;
  case 6:
    u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    break;
  case 7:
    u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
{
  AdvectCtx *user;

  PetscFunctionBeginUser;
  PetscCall(PetscNew(&user));
  ctx->physics.sample  = PhysicsSample_Advect;
  ctx->physics.flux    = PhysicsFlux_Advect;
  ctx->physics.destroy = PhysicsDestroy_SimpleFree;
  ctx->physics.user    = user;
  ctx->physics.dof     = 1;
  PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
  user->a = 1;
  PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
  {
    PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
  }
  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* --------------------------------- Finite Volume Solver ----------------------------------- */

static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
  PetscReal    hf, hs, cfl_idt = 0;
  PetscScalar *x, *f, *r, *min, *alpha, *gamma;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
  hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
  hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(DMDAVecGetArray(da, F, &f));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
  PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

  if (ctx->bctype == FVBC_OUTFLOW) {
    for (i = xs - 2; i < 0; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    }
    for (i = Mx; i < xs + xm + 2; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *u;
    if (i < sf || i > fs + 1) {
      u        = &ctx->u[0];
      alpha[0] = 1.0 / 6.0;
      gamma[0] = 1.0 / 3.0;
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hs));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
      }
    } else if (i == sf) {
      u        = &ctx->u[0];
      alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
      gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
      }
    } else if (i == sf + 1) {
      u        = &ctx->u[0];
      alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
      gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
      }
    } else if (i > sf + 1 && i < fs) {
      u        = &ctx->u[0];
      alpha[0] = 1.0 / 6.0;
      gamma[0] = 1.0 / 3.0;
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
      }
    } else if (i == fs) {
      u        = &ctx->u[0];
      alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
      gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
      }
    } else if (i == fs + 1) {
      u        = &ctx->u[0];
      alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
      gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(DMDAVecRestoreArray(da, F, &f));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscCall(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
  if (0) {
    /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
    PetscReal dt, tnow;
    PetscCall(TSGetTimeStep(ts, &dt));
    PetscCall(TSGetTime(ts, &tnow));
    if (dt > 0.5 / ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
  }
  PetscCall(PetscFree4(r, min, alpha, gamma));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs;
  PetscReal    hf, hs;
  PetscScalar *x, *f, *r, *min, *alpha, *gamma;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
  hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
  hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(VecGetArray(F, &f));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
  PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

  if (ctx->bctype == FVBC_OUTFLOW) {
    for (i = xs - 2; i < 0; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    }
    for (i = Mx; i < xs + xm + 2; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *u;
    if (i < sf) {
      u        = &ctx->u[0];
      alpha[0] = 1.0 / 6.0;
      gamma[0] = 1.0 / 3.0;
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
        islow++;
      }
    } else if (i == sf) {
      u        = &ctx->u[0];
      alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
      gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
      }
    } else if (i == fs) {
      u        = &ctx->u[0];
      alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
      gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
        islow++;
      }
    } else if (i == fs + 1) {
      u        = &ctx->u[0];
      alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
      gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
        islow++;
      }
    } else if (i > fs + 1) {
      u        = &ctx->u[0];
      alpha[0] = 1.0 / 6.0;
      gamma[0] = 1.0 / 3.0;
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
        islow++;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(VecRestoreArray(F, &f));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscCall(PetscFree4(r, min, alpha, gamma));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
  PetscReal    hf, hs;
  PetscScalar *x, *f, *r, *min, *alpha, *gamma;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
  hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
  hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(VecGetArray(F, &f));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
  PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

  if (ctx->bctype == FVBC_OUTFLOW) {
    for (i = xs - 2; i < 0; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    }
    for (i = Mx; i < xs + xm + 2; i++) {
      for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *u;
    if (i == sf) {
      u        = &ctx->u[0];
      alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
      gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
        ifast++;
      }
    } else if (i == sf + 1) {
      u        = &ctx->u[0];
      alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
      gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
        ifast++;
      }
    } else if (i > sf + 1 && i < fs) {
      u        = &ctx->u[0];
      alpha[0] = 1.0 / 6.0;
      gamma[0] = 1.0 / 3.0;
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
        ifast++;
      }
    } else if (i == fs) {
      u        = &ctx->u[0];
      alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
      gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
      for (j = 0; j < dof; j++) {
        r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
        min[j] = PetscMin(r[j], 2.0);
        u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
      }
      PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(VecRestoreArray(F, &f));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscCall(PetscFree4(r, min, alpha, gamma));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */

PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
{
  PetscScalar   *u, *uj, xj, xi;
  PetscInt       i, j, k, dof, xs, xm, Mx, count_slow, count_fast;
  const PetscInt N = 200;

  PetscFunctionBeginUser;
  PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
  PetscCall(DMDAVecGetArray(da, U, &u));
  PetscCall(PetscMalloc1(dof, &uj));
  const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
  const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
  count_slow         = Mx / (1 + ctx->hratio);
  count_fast         = Mx - count_slow;
  for (i = xs; i < xs + xm; i++) {
    if (i * hs + 0.5 * hs < (ctx->xmax - ctx->xmin) * 0.25) {
      xi = ctx->xmin + 0.5 * hs + i * hs;
      /* Integrate over cell i using trapezoid rule with N points. */
      for (k = 0; k < dof; k++) u[i * dof + k] = 0;
      for (j = 0; j < N + 1; j++) {
        xj = xi + hs * (j - N / 2) / (PetscReal)N;
        PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
        for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
      }
    } else if ((ctx->xmax - ctx->xmin) * 0.25 + (i - count_slow / 2) * hf + 0.5 * hf < (ctx->xmax - ctx->xmin) * 0.75) {
      xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.25 + 0.5 * hf + (i - count_slow / 2) * hf;
      /* Integrate over cell i using trapezoid rule with N points. */
      for (k = 0; k < dof; k++) u[i * dof + k] = 0;
      for (j = 0; j < N + 1; j++) {
        xj = xi + hf * (j - N / 2) / (PetscReal)N;
        PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
        for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
      }
    } else {
      xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.75 + 0.5 * hs + (i - count_slow / 2 - count_fast) * hs;
      /* Integrate over cell i using trapezoid rule with N points. */
      for (k = 0; k < dof; k++) u[i * dof + k] = 0;
      for (j = 0; j < N + 1; j++) {
        xj = xi + hs * (j - N / 2) / (PetscReal)N;
        PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
        for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, U, &u));
  PetscCall(PetscFree(uj));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
{
  PetscReal          xmin, xmax;
  PetscScalar        sum, tvsum, tvgsum;
  const PetscScalar *x;
  PetscInt           imin, imax, Mx, i, j, xs, xm, dof;
  Vec                Xloc;
  PetscBool          iascii;

  PetscFunctionBeginUser;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
    PetscCall(DMGetLocalVector(da, &Xloc));
    PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
    PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    tvsum = 0;
    for (i = xs; i < xs + xm; i++) {
      for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
    }
    PetscCall(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da)));
    PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
    PetscCall(DMRestoreLocalVector(da, &Xloc));

    PetscCall(VecMin(X, &imin, &xmin));
    PetscCall(VecMax(X, &imax, &xmax));
    PetscCall(VecSum(X, &sum));
    PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%g,%g] with minimum at %" PetscInt_FMT ", mean %g, ||x||_TV %g\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvgsum / Mx)));
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
{
  Vec                Y;
  PetscInt           i, Mx, count_slow = 0, count_fast = 0;
  const PetscScalar *ptr_X, *ptr_Y;

  PetscFunctionBeginUser;
  PetscCall(VecGetSize(X, &Mx));
  PetscCall(VecDuplicate(X, &Y));
  PetscCall(FVSample(ctx, da, t, Y));
  const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
  const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
  count_slow         = (PetscReal)Mx / (1.0 + ctx->hratio);
  count_fast         = Mx - count_slow;
  PetscCall(VecGetArrayRead(X, &ptr_X));
  PetscCall(VecGetArrayRead(Y, &ptr_Y));
  for (i = 0; i < Mx; i++) {
    if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
    else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
  }
  PetscCall(VecRestoreArrayRead(X, &ptr_X));
  PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
  PetscCall(VecDestroy(&Y));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char *argv[])
{
  char              physname[256] = "advect", final_fname[256] = "solution.m";
  PetscFunctionList physics = 0;
  MPI_Comm          comm;
  TS                ts;
  DM                da;
  Vec               X, X0, R;
  FVCtx             ctx;
  PetscInt          i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, *index_slow, *index_fast;
  PetscBool         view_final = PETSC_FALSE;
  PetscReal         ptime;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, 0, help));
  comm = PETSC_COMM_WORLD;
  PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

  /* Register physical models to be available on the command line */
  PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));

  ctx.comm   = comm;
  ctx.cfl    = 0.9;
  ctx.bctype = FVBC_PERIODIC;
  ctx.xmin   = -1.0;
  ctx.xmax   = 1.0;
  PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
  PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
  PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
  PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
  PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
  PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
  PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
  PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
  PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
  PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
  PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
  PetscOptionsEnd();

  /* Choose the physics from the list of registered models */
  {
    PetscErrorCode (*r)(FVCtx *);
    PetscCall(PetscFunctionListFind(physics, physname, &r));
    PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
    /* Create the physics, will set the number of fields and their names */
    PetscCall((*r)(&ctx));
  }

  /* Create a DMDA to manage the parallel grid */
  PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
  PetscCall(DMSetFromOptions(da));
  PetscCall(DMSetUp(da));
  /* Inform the DMDA of the field names provided by the physics. */
  /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
  for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

  /* Set coordinates of cell centers */
  PetscCall(DMDASetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax - ctx.xmin) / Mx, 0, 0, 0, 0));

  /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
  PetscCall(PetscMalloc3(dof, &ctx.u, dof, &ctx.flux, dof, &ctx.speeds));

  /* Create a vector to store the solution and to save the initial state */
  PetscCall(DMCreateGlobalVector(da, &X));
  PetscCall(VecDuplicate(X, &X0));
  PetscCall(VecDuplicate(X, &R));

  /* create index for slow parts and fast parts*/
  count_slow = Mx / (1 + ctx.hratio);
  PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio) is even");
  count_fast = Mx - count_slow;
  ctx.sf     = count_slow / 2;
  ctx.fs     = ctx.sf + count_fast;
  PetscCall(PetscMalloc1(xm * dof, &index_slow));
  PetscCall(PetscMalloc1(xm * dof, &index_fast));
  for (i = xs; i < xs + xm; i++) {
    if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1)
      for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    else
      for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
  }
  PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
  PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));

  /* Create a time-stepping object */
  PetscCall(TSCreate(comm, &ts));
  PetscCall(TSSetDM(ts, da));
  PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
  PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
  PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
  PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
  PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));

  PetscCall(TSSetType(ts, TSMPRK));
  PetscCall(TSSetMaxTime(ts, 10));
  PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

  /* Compute initial conditions and starting time step */
  PetscCall(FVSample(&ctx, da, 0, X0));
  PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
  PetscCall(VecCopy(X0, X));                            /* The function value was not used so we set X=X0 again */
  PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
  PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
  PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
  {
    PetscInt           steps;
    PetscScalar        mass_initial, mass_final, mass_difference, mass_differenceg;
    const PetscScalar *ptr_X, *ptr_X0;
    const PetscReal    hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow;
    const PetscReal    hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast;
    PetscCall(TSSolve(ts, X));
    PetscCall(TSGetSolveTime(ts, &ptime));
    PetscCall(TSGetStepNumber(ts, &steps));
    /* calculate the total mass at initial time and final time */
    mass_initial = 0.0;
    mass_final   = 0.0;
    PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
    PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
    for (i = xs; i < xs + xm; i++) {
      if (i < ctx.sf || i > ctx.fs - 1) {
        for (k = 0; k < dof; k++) {
          mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
          mass_final   = mass_final + hs * ptr_X[i * dof + k];
        }
      } else {
        for (k = 0; k < dof; k++) {
          mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
          mass_final   = mass_final + hf * ptr_X[i * dof + k];
        }
      }
    }
    PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
    PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
    mass_difference = mass_final - mass_initial;
    PetscCall(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
    PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
    PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    if (ctx.exact) {
      PetscReal nrm1 = 0;
      PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1));
      PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    }
    if (ctx.simulation) {
      PetscReal          nrm1 = 0;
      PetscViewer        fd;
      char               filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
      Vec                XR;
      PetscBool          flg;
      const PetscScalar *ptr_XR;
      PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
      PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
      PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
      PetscCall(VecDuplicate(X0, &XR));
      PetscCall(VecLoad(XR, fd));
      PetscCall(PetscViewerDestroy(&fd));
      PetscCall(VecGetArrayRead(X, &ptr_X));
      PetscCall(VecGetArrayRead(XR, &ptr_XR));
      for (i = 0; i < Mx; i++) {
        if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] - ptr_XR[i]);
        else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]);
      }
      PetscCall(VecRestoreArrayRead(X, &ptr_X));
      PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
      PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
      PetscCall(VecDestroy(&XR));
    }
  }

  PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
  if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
  if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
  if (draw & 0x4) {
    Vec Y;
    PetscCall(VecDuplicate(X, &Y));
    PetscCall(FVSample(&ctx, da, ptime, Y));
    PetscCall(VecAYPX(Y, -1, X));
    PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    PetscCall(VecDestroy(&Y));
  }

  if (view_final) {
    PetscViewer viewer;
    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    PetscCall(VecView(X, viewer));
    PetscCall(PetscViewerPopFormat(viewer));
    PetscCall(PetscViewerDestroy(&viewer));
  }

  /* Clean up */
  PetscCall((*ctx.physics.destroy)(ctx.physics.user));
  for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
  PetscCall(PetscFree3(ctx.u, ctx.flux, ctx.speeds));
  PetscCall(ISDestroy(&ctx.iss));
  PetscCall(ISDestroy(&ctx.isf));
  PetscCall(VecDestroy(&X));
  PetscCall(VecDestroy(&X0));
  PetscCall(VecDestroy(&R));
  PetscCall(DMDestroy(&da));
  PetscCall(TSDestroy(&ts));
  PetscCall(PetscFree(index_slow));
  PetscCall(PetscFree(index_fast));
  PetscCall(PetscFunctionListDestroy(&physics));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

    build:
      requires: !complex

    test:
      args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0

    test:
      suffix: 2
      args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
      output_file: output/ex7_1.out

    test:
      suffix: 3
      args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0

    test:
      suffix: 4
      args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
      output_file: output/ex7_3.out

    test:
      suffix: 5
      nsize: 2
      args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
      output_file: output/ex7_3.out
TEST*/