/*
  Note:
    -hratio is the ratio between mesh size of coarse grids and fine grids
*/

static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
                           "  advect      - Constant coefficient scalar advection\n"
                           "                u_t       + (a*u)_x               = 0\n"
                           "  shallow     - 1D Shallow water equations (Saint Venant System)\n"
                           "                h_t + (q)_x = 0 \n"
                           "                q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0  \n"
                           "                where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n"
                           "  for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
                           "                hxs  = hratio*hxf \n"
                           "  where hxs and hxf are the grid spacings for coarse and fine grids respectively.\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"
                           "  bc_type     - Boundary condition for the problem, options are: periodic, outflow, inflow "
                           "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n"
                           "The problem size should be set with -da_grid_x M\n\n";

/*
  Example:
     ./ex4 -da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
     ./ex4 -da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 2.5 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
     ./ex4 -da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
     ./ex4 -da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
     ./ex4 -da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 5.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0

  Contributed by: Aidan Hamilton <aidan@udel.edu>
*/

#include <petscts.h>
#include <petscdm.h>
#include <petscdmda.h>
#include <petscdraw.h>
#include "finitevolume1d.h"
#include <petsc/private/kernels/blockinvert.h>

static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
{
  PetscReal range = xmax - xmin;
  return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
}
static inline PetscReal MaxAbs(PetscReal a, PetscReal b)
{
  return (PetscAbs(a) > PetscAbs(b)) ? a : b;
}
/* --------------------------------- Advection ----------------------------------- */
typedef struct {
  PetscReal a; /* advective velocity */
} AdvectCtx;

static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
{
  AdvectCtx *ctx = (AdvectCtx *)vctx;
  PetscReal  speed;

  PetscFunctionBeginUser;
  speed     = ctx->a;
  flux[0]   = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
  *maxspeed = speed;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
{
  AdvectCtx *ctx = (AdvectCtx *)vctx;

  PetscFunctionBeginUser;
  X[0]      = 1.;
  Xi[0]     = 1.;
  speeds[0] = ctx->a;
  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->physics2.sample2         = PhysicsSample_Advect;
  ctx->physics2.riemann2        = PhysicsRiemann_Advect;
  ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
  ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
  ctx->physics2.user            = user;
  ctx->physics2.dof             = 1;

  PetscCall(PetscStrallocpy("u", &ctx->physics2.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);
}
/* --------------------------------- Shallow Water ----------------------------------- */

typedef struct {
  PetscReal gravity;
} ShallowCtx;

static inline void ShallowFlux(ShallowCtx *phys, const PetscScalar *u, PetscScalar *f)
{
  f[0] = u[1];
  f[1] = PetscSqr(u[1]) / u[0] + 0.5 * phys->gravity * PetscSqr(u[0]);
}

static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
{
  ShallowCtx *phys = (ShallowCtx *)vctx;
  PetscScalar g    = phys->gravity, ustar[2], cL, cR, c, cstar;
  struct {
    PetscScalar h, u;
  } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star;
  PetscInt i;

  PetscFunctionBeginUser;
  PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
  cL = PetscSqrtScalar(g * L.h);
  cR = PetscSqrtScalar(g * R.h);
  c  = PetscMax(cL, cR);
  {
    /* Solve for star state */
    const PetscInt maxits = 50;
    PetscScalar    tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */
    h0 = h;
    for (i = 0; i < maxits; i++) {
      PetscScalar fr, fl, dfr, dfl;
      fl  = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */
                      : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h);         /* rarefaction */
      fr  = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */
                      : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h);         /* rarefaction */
      res = R.u - L.u + fr + fl;
      PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
      if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h - h0) < PETSC_SQRT_MACHINE_EPSILON)) {
        star.h = h;
        star.u = L.u - fl;
        goto converged;
      } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
        h = 0.8 * h0 + 0.2 * h;
        continue;
      }
      /* Accept the last step and take another */
      res0 = res;
      h0   = h;
      dfl  = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h);
      dfr  = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h);
      tmp  = h - res / (dfr + dfl);
      if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
      else h = tmp;
      PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h);
    }
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i);
  }
converged:
  cstar = PetscSqrtScalar(g * star.h);
  if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */
    PetscScalar ufan[2];
    ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL);
    ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0];
    ShallowFlux(phys, ufan, flux);
  } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */
    PetscScalar ufan[2];
    ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR);
    ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0];
    ShallowFlux(phys, ufan, flux);
  } else if ((L.h >= star.h && L.u - c >= 0) || (L.h < star.h && (star.h * star.u - L.h * L.u) / (star.h - L.h) > 0)) {
    /* 1-wave is right-travelling shock (supersonic) */
    ShallowFlux(phys, uL, flux);
  } else if ((star.h <= R.h && R.u + c <= 0) || (star.h > R.h && (R.h * R.u - star.h * star.h) / (R.h - star.h) < 0)) {
    /* 2-wave is left-travelling shock (supersonic) */
    ShallowFlux(phys, uR, flux);
  } else {
    ustar[0] = star.h;
    ustar[1] = star.h * star.u;
    ShallowFlux(phys, ustar, flux);
  }
  *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
{
  ShallowCtx *phys = (ShallowCtx *)vctx;
  PetscScalar g    = phys->gravity, fL[2], fR[2], s;
  struct {
    PetscScalar h, u;
  } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};
  PetscReal tol = 1e-6;

  PetscFunctionBeginUser;
  /* Positivity preserving modification*/
  if (L.h < tol) L.u = 0.0;
  if (R.h < tol) R.u = 0.0;

  /*simple pos preserve limiter*/
  if (L.h < 0) L.h = 0;
  if (R.h < 0) R.h = 0;

  ShallowFlux(phys, uL, fL);
  ShallowFlux(phys, uR, fR);

  s         = PetscMax(PetscAbs(L.u) + PetscSqrtScalar(g * L.h), PetscAbs(R.u) + PetscSqrtScalar(g * R.h));
  flux[0]   = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (L.h - R.h);
  flux[1]   = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
  *maxspeed = s;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
{
  PetscInt i, j;

  PetscFunctionBeginUser;
  for (i = 0; i < m; i++) {
    for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j);
    speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
{
  ShallowCtx *phys = (ShallowCtx *)vctx;
  PetscReal   c;
  PetscReal   tol = 1e-6;

  PetscFunctionBeginUser;
  c = PetscSqrtScalar(u[0] * phys->gravity);

  if (u[0] < tol) { /*Use conservative variables*/
    X[0 * 2 + 0] = 1;
    X[0 * 2 + 1] = 0;
    X[1 * 2 + 0] = 0;
    X[1 * 2 + 1] = 1;
  } else {
    speeds[0]    = u[1] / u[0] - c;
    speeds[1]    = u[1] / u[0] + c;
    X[0 * 2 + 0] = 1;
    X[0 * 2 + 1] = speeds[0];
    X[1 * 2 + 0] = 1;
    X[1 * 2 + 1] = speeds[1];
  }

  PetscCall(PetscArraycpy(Xi, X, 4));
  PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsSample_Shallow(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
{
  PetscFunctionBeginUser;
  PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0");
  switch (initial) {
  case 0:
    u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */
    u[1] = (x < 0) ? 0 : 0;
    break;
  case 1:
    u[0] = (x < 10) ? 1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */
    u[1] = (x < 10) ? 2.5 : 0;
    break;
  case 2:
    u[0] = (x < 25) ? 1 : 1;
    u[1] = (x < 25) ? -5 : 5;
    break;
  case 3:
    u[0] = (x < 20) ? 1 : 1e-12;
    u[1] = (x < 20) ? 0 : 0;
    break;
  case 4:
    u[0] = (x < 30) ? 1e-12 : 1;
    u[1] = (x < 30) ? 0 : 0;
    break;
  case 5:
    u[0] = (x < 25) ? 0.1 : 0.1;
    u[1] = (x < 25) ? -0.3 : 0.3;
    break;
  case 6:
    u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x);
    u[1] = 1 * u[0];
    break;
  case 7:
    if (x < -0.1) {
      u[0] = 1e-9;
      u[1] = 0.0;
    } else if (x < 0.1) {
      u[0] = 1.0;
      u[1] = 0.0;
    } else {
      u[0] = 1e-9;
      u[1] = 0.0;
    }
    break;
  case 8:
    if (x < -0.1) {
      u[0] = 2;
      u[1] = 0.0;
    } else if (x < 0.1) {
      u[0] = 3.0;
      u[1] = 0.0;
    } else {
      u[0] = 2;
      u[1] = 0.0;
    }
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* Implements inflow conditions for the given initial conditions. Which conditions are actually enforced depends on
   on the results of PhysicsSetInflowType_Shallow. */
static PetscErrorCode PhysicsInflow_Shallow(void *vctx, PetscReal t, PetscReal x, PetscReal *u)
{
  FVCtx *ctx = (FVCtx *)vctx;

  PetscFunctionBeginUser;
  if (ctx->bctype == FVBC_INFLOW) {
    switch (ctx->initial) {
    case 0:
    case 1:
    case 2:
    case 3:
    case 4:
    case 5:
      u[0] = 0;
      u[1] = 0.0; /* Left boundary conditions */
      u[2] = 0;
      u[3] = 0.0; /* Right boundary conditions */
      break;
    case 6:
      u[0] = 0;
      u[1] = 0.0; /* Left boundary conditions */
      u[2] = 0;
      u[3] = 0.0; /* Right boundary conditions */
      break;
    case 7:
      u[0] = 0;
      u[1] = 0.0; /* Left boundary conditions */
      u[2] = 0;
      u[3] = 0.0; /* Right boundary conditions */
      break;
    case 8:
      u[0] = 0;
      u[1] = 1.0; /* Left boundary conditions */
      u[2] = 0;
      u[3] = -1.0; /* Right boundary conditions */
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */
static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx)
{
  PetscFunctionBeginUser;
  switch (ctx->initial) {
  case 0:
  case 1:
  case 2:
  case 3:
  case 4:
  case 5:
  case 6:
  case 7:
  case 8: /* Fix left and right momentum, height is outflow */
    ctx->physics2.bcinflowindex[0] = PETSC_FALSE;
    ctx->physics2.bcinflowindex[1] = PETSC_TRUE;
    ctx->physics2.bcinflowindex[2] = PETSC_FALSE;
    ctx->physics2.bcinflowindex[3] = PETSC_TRUE;
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
{
  ShallowCtx       *user;
  PetscFunctionList rlist = 0, rclist = 0;
  char              rname[256] = "rusanov", rcname[256] = "characteristic";

  PetscFunctionBeginUser;
  PetscCall(PetscNew(&user));
  ctx->physics2.sample2         = PhysicsSample_Shallow;
  ctx->physics2.inflow          = PhysicsInflow_Shallow;
  ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
  ctx->physics2.riemann2        = PhysicsRiemann_Shallow_Rusanov;
  ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow;
  ctx->physics2.user            = user;
  ctx->physics2.dof             = 2;

  PetscCall(PetscMalloc1(2 * (ctx->physics2.dof), &ctx->physics2.bcinflowindex));
  PetscCall(PhysicsSetInflowType_Shallow(ctx));

  PetscCall(PetscStrallocpy("height", &ctx->physics2.fieldname[0]));
  PetscCall(PetscStrallocpy("momentum", &ctx->physics2.fieldname[1]));

  user->gravity = 9.81;

  PetscCall(RiemannListAdd_2WaySplit(&rlist, "exact", PhysicsRiemann_Shallow_Exact));
  PetscCall(RiemannListAdd_2WaySplit(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov));
  PetscCall(ReconstructListAdd_2WaySplit(&rclist, "characteristic", PhysicsCharacteristic_Shallow));
  PetscCall(ReconstructListAdd_2WaySplit(&rclist, "conservative", PhysicsCharacteristic_Conservative));
  PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", "");
  PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL));
  PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
  PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
  PetscOptionsEnd();
  PetscCall(RiemannListFind_2WaySplit(rlist, rname, &ctx->physics2.riemann2));
  PetscCall(ReconstructListFind_2WaySplit(rclist, rcname, &ctx->physics2.characteristic2));
  PetscCall(PetscFunctionListDestroy(&rlist));
  PetscCall(PetscFunctionListDestroy(&rclist));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCheck(ctx->physics2.sample2, 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));
  hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  for (i = xs; i < xs + xm; i++) {
    if (i < ctx->sf) {
      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->physics2.sample2)(ctx->physics2.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 (i < ctx->fs) {
      xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * 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->physics2.sample2)(ctx->physics2.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->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * 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->physics2.sample2)(ctx->physics2.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 SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
{
  Vec                Y;
  PetscInt           i, Mx;
  const PetscScalar *ptr_X, *ptr_Y;
  PetscReal          hs, hf;

  PetscFunctionBeginUser;
  PetscCall(VecGetSize(X, &Mx));
  PetscCall(VecDuplicate(X, &Y));
  PetscCall(FVSample_2WaySplit(ctx, da, t, Y));
  hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  PetscCall(VecGetArrayRead(X, &ptr_X));
  PetscCall(VecGetArrayRead(Y, &ptr_Y));
  for (i = 0; i < Mx; i++) {
    if (i < ctx->sf || i > ctx->fs - 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);
}

PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
  PetscReal    hxf, hxs, cfl_idt = 0;
  PetscScalar *x, *f, *slope;
  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                              */
  hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  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(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points                                           */
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

  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];
    }
  }

  if (ctx->bctype == FVBC_INFLOW) {
    /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    pages 137-138 for the scheme. */
    if (xs == 0) { /* Left Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[j]) {
          for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
        } else {
          for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
        }
      }
    }
    if (xs + xm == Mx) { /* Right Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[dof + j]) {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
        } else {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
        }
      }
    }
  }

  for (i = xs - 1; i < xs + xm + 1; i++) {
    struct _LimitInfo info;
    PetscScalar      *cjmpL, *cjmpR;
    /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    cjmpL = &ctx->cjmpLR[0];
    cjmpR = &ctx->cjmpLR[dof];
    for (j = 0; j < dof; j++) {
      PetscScalar jmpL, jmpR;
      jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
      jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
      for (k = 0; k < dof; k++) {
        cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
        cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
      }
    }
    /* Apply limiter to the left and right characteristic jumps */
    info.m   = dof;
    info.hxs = hxs;
    info.hxf = hxf;
    (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    for (j = 0; j < dof; j++) {
      PetscScalar tmp = 0;
      for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
      slope[i * dof + j] = tmp;
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *uL, *uR;
    uL = &ctx->uLR[0];
    uR = &ctx->uLR[dof];
    if (i < sf) { /* slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
      }
    } else if (i == sf) { /* interface between the slow region and the fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
      }
    } else if (i > sf && i < fs) { /* fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
      }
    } else if (i == fs) { /* interface between the fast region and the slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
      }
    } else { /* slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(DMDAVecRestoreArray(da, F, &f));
  PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
  if (0) {
    /* We need to 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) {
      if (1) {
        PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
      } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
  PetscReal    hxs, hxf, cfl_idt = 0;
  PetscScalar *x, *f, *slope;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
  hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F));
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(VecGetArray(F, &f));
  PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

  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];
    }
  }
  if (ctx->bctype == FVBC_INFLOW) {
    /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    pages 137-138 for the scheme. */
    if (xs == 0) { /* Left Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
          for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
        } else {
          for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
        }
      }
    }
    if (xs + xm == Mx) { /* Right Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
        } else {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
        }
      }
    }
  }
  for (i = xs - 1; i < xs + xm + 1; i++) {
    struct _LimitInfo info;
    PetscScalar      *cjmpL, *cjmpR;
    if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
      /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
      PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
      /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
      PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
      cjmpL = &ctx->cjmpLR[0];
      cjmpR = &ctx->cjmpLR[dof];
      for (j = 0; j < dof; j++) {
        PetscScalar jmpL, jmpR;
        jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
        jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
        for (k = 0; k < dof; k++) {
          cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
          cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
        }
      }
      /* Apply limiter to the left and right characteristic jumps */
      info.m   = dof;
      info.hxs = hxs;
      info.hxf = hxf;
      (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
      for (j = 0; j < dof; j++) {
        PetscScalar tmp = 0;
        for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
        slope[i * dof + j] = tmp;
      }
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *uL, *uR;
    uL = &ctx->uLR[0];
    uR = &ctx->uLR[dof];
    if (i < sf - lsbwidth) { /* slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
    }
    if (i == fs + rsbwidth) { /* slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i > fs + rsbwidth) { /* slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(VecRestoreArray(F, &f));
  PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
  PetscReal    hxs, hxf;
  PetscScalar *x, *f, *slope;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
  hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F));
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(VecGetArray(F, &f));
  PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

  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];
    }
  }
  if (ctx->bctype == FVBC_INFLOW) {
    /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    pages 137-138 for the scheme. */
    if (xs == 0) { /* Left Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
          for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
        } else {
          for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
        }
      }
    }
    if (xs + xm == Mx) { /* Right Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
        } else {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
        }
      }
    }
  }
  for (i = xs - 1; i < xs + xm + 1; i++) {
    struct _LimitInfo info;
    PetscScalar      *cjmpL, *cjmpR;
    if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
      /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
      PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
      /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
      PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
      cjmpL = &ctx->cjmpLR[0];
      cjmpR = &ctx->cjmpLR[dof];
      for (j = 0; j < dof; j++) {
        PetscScalar jmpL, jmpR;
        jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
        jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
        for (k = 0; k < dof; k++) {
          cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
          cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
        }
      }
      /* Apply limiter to the left and right characteristic jumps */
      info.m   = dof;
      info.hxs = hxs;
      info.hxf = hxf;
      (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
      for (j = 0; j < dof; j++) {
        PetscScalar tmp = 0;
        for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
        slope[i * dof + j] = tmp;
      }
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *uL, *uR;
    uL = &ctx->uLR[0];
    uR = &ctx->uLR[dof];
    if (i == sf - lsbwidth) {
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i > sf - lsbwidth && i < sf) {
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i == sf) { /* interface between the slow region and the fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
    }
    if (i == fs) { /* interface between the fast region and the slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i > fs && i < fs + rsbwidth) {
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
        islow++;
      }
    }
    if (i == fs + rsbwidth) {
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(VecRestoreArray(F, &f));
  PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
{
  FVCtx       *ctx = (FVCtx *)vctx;
  PetscInt     i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
  PetscReal    hxs, hxf;
  PetscScalar *x, *f, *slope;
  Vec          Xloc;
  DM           da;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &Xloc));
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
  hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
  hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
  PetscCall(VecZeroEntries(F));
  PetscCall(DMDAVecGetArray(da, Xloc, &x));
  PetscCall(VecGetArray(F, &f));
  PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
  PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

  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];
    }
  }
  if (ctx->bctype == FVBC_INFLOW) {
    /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    pages 137-138 for the scheme. */
    if (xs == 0) { /* Left Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
          for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
        } else {
          for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
        }
      }
    }
    if (xs + xm == Mx) { /* Right Boundary */
      PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
      for (j = 0; j < dof; j++) {
        if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
        } else {
          for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
        }
      }
    }
  }
  for (i = xs - 1; i < xs + xm + 1; i++) {
    struct _LimitInfo info;
    PetscScalar      *cjmpL, *cjmpR;
    if (i > sf - 2 && i < fs + 1) {
      PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
      PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
      cjmpL = &ctx->cjmpLR[0];
      cjmpR = &ctx->cjmpLR[dof];
      for (j = 0; j < dof; j++) {
        PetscScalar jmpL, jmpR;
        jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
        jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
        for (k = 0; k < dof; k++) {
          cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
          cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
        }
      }
      /* Apply limiter to the left and right characteristic jumps */
      info.m   = dof;
      info.hxs = hxs;
      info.hxf = hxf;
      (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
      for (j = 0; j < dof; j++) {
        PetscScalar tmp = 0;
        for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
        slope[i * dof + j] = tmp;
      }
    }
  }

  for (i = xs; i < xs + xm + 1; i++) {
    PetscReal    maxspeed;
    PetscScalar *uL, *uR;
    uL = &ctx->uLR[0];
    uR = &ctx->uLR[dof];
    if (i == sf) { /* interface between the slow region and the fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
        ifast++;
      }
    }
    if (i > sf && i < fs) { /* fast region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
      }
      if (i < xs + xm) {
        for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
        ifast++;
      }
    }
    if (i == fs) { /* interface between the fast region and the slow region */
      for (j = 0; j < dof; j++) {
        uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
        uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
      }
      PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
      if (i > xs) {
        for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
      }
    }
  }
  PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
  PetscCall(VecRestoreArray(F, &f));
  PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
  PetscCall(DMRestoreLocalVector(da, &Xloc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char *argv[])
{
  char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
  PetscFunctionList limiters = 0, 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, islowbuffer = 0, *index_slow, *index_fast, *index_slowbuffer;
  PetscBool         view_final = PETSC_FALSE;
  PetscReal         ptime, maxtime;

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

  /* Register limiters to be available on the command line */
  PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
  PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
  PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
  PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
  PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
  PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
  PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
  PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));

  /* Register physical models to be available on the command line */
  PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow));
  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;
  ctx.initial    = 1;
  ctx.hratio     = 2;
  maxtime        = 10.0;
  ctx.simulation = PETSC_FALSE;
  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(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
  PetscCall(PetscOptionsFList("-physics", "Name of physics model to use", "", physics, physname, physname, sizeof(physname), 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 limiter from the list of registered limiters */
  PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2));
  PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);

  /* 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.physics2.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.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.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(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
  PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
  PetscCall(PetscMalloc1(2 * dof, &ctx.ub));

  /* 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 + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
  count_slow = Mx * 3 / (3 + ctx.hratio); // compute Mx / (1.0 + ctx.hratio / 3.0);
  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+hratio/3) 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));
  PetscCall(PetscMalloc1(8 * dof, &index_slowbuffer));
  ctx.lsbwidth = 4;
  ctx.rsbwidth = 4;

  for (i = xs; i < xs + xm; i++) {
    if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
      for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
      for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = 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));
  PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));

  /* Create a time-stepping object */
  PetscCall(TSCreate(comm, &ts));
  PetscCall(TSSetDM(ts, da));
  PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
  PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
  PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
  PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
  PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
  PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
  PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));

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

  /* Compute initial conditions and starting time step */
  PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
  PetscCall(FVRHSFunction_2WaySplit(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) * 3.0 / 4.0 / count_slow;
    const PetscReal    hf = (ctx.xmax - ctx.xmin) / 4.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;
    PetscCallMPI(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));
    PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1.0 / ctx.cfl_idt)));
    if (ctx.exact) {
      PetscReal nrm1 = 0;
      PetscCall(SolutionErrorNorms_2WaySplit(&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 = xs; i < xs + xm; i++) {
        if (i < ctx.sf || i > ctx.fs - 1)
          for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
        else
          for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
      }
      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_2WaySplit(&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.physics2.destroy)(ctx.physics2.user));
  for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
  PetscCall(PetscFree(ctx.physics2.bcinflowindex));
  PetscCall(PetscFree(ctx.ub));
  PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
  PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
  PetscCall(VecDestroy(&X));
  PetscCall(VecDestroy(&X0));
  PetscCall(VecDestroy(&R));
  PetscCall(DMDestroy(&da));
  PetscCall(TSDestroy(&ts));
  PetscCall(ISDestroy(&ctx.iss));
  PetscCall(ISDestroy(&ctx.isf));
  PetscCall(ISDestroy(&ctx.issb));
  PetscCall(PetscFree(index_slow));
  PetscCall(PetscFree(index_fast));
  PetscCall(PetscFree(index_slowbuffer));
  PetscCall(PetscFunctionListDestroy(&limiters));
  PetscCall(PetscFunctionListDestroy(&physics));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

    build:
      requires: !complex !single
      depends: finitevolume1d.c

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

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

    test:
      suffix: 3
      args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
      output_file: output/ex4_3.out

    test:
      suffix: 4
      nsize: 2
      args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
      output_file: output/ex4_3.out

    test:
      suffix: 5
      nsize: 4
      args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
      output_file: output/ex4_3.out
TEST*/