/*
  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(0);
}

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(0);
}

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(0);
}

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(0);
}
/* --------------------------------- 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,1,"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(0);
}

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(0);
}

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(0);
}

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(0);
}

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(0);
}

/* 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(0);
}

/* 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(0);
}

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;

  PetscMalloc1(2*(ctx->physics2.dof),&ctx->physics2.bcinflowindex);
  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(0);
}

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(0);
}

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(0);
}

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 */
      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 */
      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(MPI_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,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
    }
  }
  PetscFunctionReturn(0);
}

/* --------------------------------- 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 */
      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 */
      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(MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da)));
  PetscFunctionReturn(0);
}

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 */
      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 */
      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(0);
}

/* --------------------------------- 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 */
      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 */
      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(0);
}

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;

  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/(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(MPI_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: 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*/