Actual source code: ex9.c

      1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
      2:                            "Solves scalar and vector problems, choose the physical model with -physics\n"
      3:                            "  advection   - Constant coefficient scalar advection\n"
      4:                            "                u_t       + (a*u)_x               = 0\n"
      5:                            "  burgers     - Burgers equation\n"
      6:                            "                u_t       + (u^2/2)_x             = 0\n"
      7:                            "  traffic     - Traffic equation\n"
      8:                            "                u_t       + (u*(1-u))_x           = 0\n"
      9:                            "  acoustics   - Acoustic wave propagation\n"
     10:                            "                u_t       + (c*z*v)_x             = 0\n"
     11:                            "                v_t       + (c/z*u)_x             = 0\n"
     12:                            "  isogas      - Isothermal gas dynamics\n"
     13:                            "                rho_t     + (rho*u)_x             = 0\n"
     14:                            "                (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
     15:                            "  shallow     - Shallow water equations\n"
     16:                            "                h_t       + (h*u)_x               = 0\n"
     17:                            "                (h*u)_t   + (h*u^2 + g*h^2/2)_x   = 0\n"
     18:                            "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
     19:                            "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
     20:                            "                the states across shocks and rarefactions\n"
     21:                            "  roe         - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
     22:                            "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
     23:                            "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
     24:                            "  conservative   - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
     25:                            "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
     26:                            "  upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
     27:                            "  and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
     28:                            "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
     29:                            "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
     30:                            "Several initial conditions can be chosen with -initial N\n\n"
     31:                            "The problem size should be set with -da_grid_x M\n\n";

     33: #include <petscts.h>
     34: #include <petscdm.h>
     35: #include <petscdmda.h>
     36: #include <petscdraw.h>

     38: #include <petsc/private/kernels/blockinvert.h>

     40: static inline PetscReal Sgn(PetscReal a)
     41: {
     42:   return (a < 0) ? -1 : 1;
     43: }
     44: static inline PetscReal Abs(PetscReal a)
     45: {
     46:   return (a < 0) ? 0 : a;
     47: }
     48: static inline PetscReal Sqr(PetscReal a)
     49: {
     50:   return a * a;
     51: }
     52: static inline PetscReal MaxAbs(PetscReal a, PetscReal b)
     53: {
     54:   return (PetscAbs(a) > PetscAbs(b)) ? a : b;
     55: }
     56: PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a, PetscReal b)
     57: {
     58:   return (PetscAbs(a) < PetscAbs(b)) ? a : b;
     59: }
     60: static inline PetscReal MinMod2(PetscReal a, PetscReal b)
     61: {
     62:   return (a * b < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscAbs(b));
     63: }
     64: static inline PetscReal MaxMod2(PetscReal a, PetscReal b)
     65: {
     66:   return (a * b < 0) ? 0 : Sgn(a) * PetscMax(PetscAbs(a), PetscAbs(b));
     67: }
     68: static inline PetscReal MinMod3(PetscReal a, PetscReal b, PetscReal c)
     69: {
     70:   return (a * b < 0 || a * c < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscMin(PetscAbs(b), PetscAbs(c)));
     71: }

     73: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
     74: {
     75:   PetscReal range = xmax - xmin;
     76:   return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
     77: }

     79: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
     80: typedef struct _LimitInfo {
     81:   PetscReal hx;
     82:   PetscInt  m;
     83: }          *LimitInfo;
     84: static void Limit_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     85: {
     86:   PetscInt i;
     87:   for (i = 0; i < info->m; i++) lmt[i] = 0;
     88: }
     89: static void Limit_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     90: {
     91:   PetscInt i;
     92:   for (i = 0; i < info->m; i++) lmt[i] = jR[i];
     93: }
     94: static void Limit_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     95: {
     96:   PetscInt i;
     97:   for (i = 0; i < info->m; i++) lmt[i] = jL[i];
     98: }
     99: static void Limit_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    100: {
    101:   PetscInt i;
    102:   for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]);
    103: }
    104: static void Limit_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    105: {
    106:   PetscInt i;
    107:   for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]);
    108: }
    109: static void Limit_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    110: {
    111:   PetscInt i;
    112:   for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i]));
    113: }
    114: static void Limit_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    115: {
    116:   PetscInt i;
    117:   for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]);
    118: }
    119: static void Limit_VanLeer(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    120: { /* phi = (t + abs(t)) / (1 + abs(t)) */
    121:   PetscInt i;
    122:   for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Abs(jR[i]) + Abs(jL[i]) * jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
    123: }
    124: static void Limit_VanAlbada(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
    125: {                                                                                                           /* phi = (t + t^2) / (1 + t^2) */
    126:   PetscInt i;
    127:   for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
    128: }
    129: static void Limit_VanAlbadaTVD(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    130: { /* phi = (t + t^2) / (1 + t^2) */
    131:   PetscInt i;
    132:   for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * jR[i] < 0) ? 0 : (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
    133: }
    134: static void Limit_Koren(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
    135: {                                                                                                       /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
    136:   PetscInt i;
    137:   for (i = 0; i < info->m; i++) lmt[i] = ((jL[i] * Sqr(jR[i]) + 2 * Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
    138: }
    139: static void Limit_KorenSym(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
    140: {                                                                                                          /* Symmetric version of above */
    141:   PetscInt i;
    142:   for (i = 0; i < info->m; i++) lmt[i] = (1.5 * (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
    143: }
    144: static void Limit_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    145: { /* Eq 11 of Cada-Torrilhon 2009 */
    146:   PetscInt i;
    147:   for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]);
    148: }
    149: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L, PetscReal R)
    150: {
    151:   return PetscMax(0, PetscMin((L + 2 * R) / 3, PetscMax(-0.5 * L, PetscMin(2 * L, PetscMin((L + 2 * R) / 3, 1.6 * R)))));
    152: }
    153: static void Limit_CadaTorrilhon2(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    154: { /* Cada-Torrilhon 2009, Eq 13 */
    155:   PetscInt i;
    156:   for (i = 0; i < info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]);
    157: }
    158: static void Limit_CadaTorrilhon3R(PetscReal r, LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    159: { /* Cada-Torrilhon 2009, Eq 22 */
    160:   /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
    161:   const PetscReal eps = 1e-7, hx = info->hx;
    162:   PetscInt        i;
    163:   for (i = 0; i < info->m; i++) {
    164:     const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r * hx);
    165:     lmt[i] = ((eta < 1 - eps) ? (jL[i] + 2 * jR[i]) / 3 : ((eta > 1 + eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]) : 0.5 * ((1 - (eta - 1) / eps) * (jL[i] + 2 * jR[i]) / 3 + (1 + (eta + 1) / eps) * CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]))));
    166:   }
    167: }
    168: static void Limit_CadaTorrilhon3R0p1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    169: {
    170:   Limit_CadaTorrilhon3R(0.1, info, jL, jR, lmt);
    171: }
    172: static void Limit_CadaTorrilhon3R1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    173: {
    174:   Limit_CadaTorrilhon3R(1, info, jL, jR, lmt);
    175: }
    176: static void Limit_CadaTorrilhon3R10(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    177: {
    178:   Limit_CadaTorrilhon3R(10, info, jL, jR, lmt);
    179: }
    180: static void Limit_CadaTorrilhon3R100(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    181: {
    182:   Limit_CadaTorrilhon3R(100, info, jL, jR, lmt);
    183: }

    185: /* --------------------------------- Finite Volume data structures ----------------------------------- */

    187: typedef enum {
    188:   FVBC_PERIODIC,
    189:   FVBC_OUTFLOW
    190: } FVBCType;
    191: static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", 0};
    192: typedef PetscErrorCode (*RiemannFunction)(void *, PetscInt, const PetscScalar *, const PetscScalar *, PetscScalar *, PetscReal *);
    193: typedef PetscErrorCode (*ReconstructFunction)(void *, PetscInt, const PetscScalar *, PetscScalar *, PetscScalar *, PetscReal *);

    195: typedef struct {
    196:   PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *);
    197:   RiemannFunction     riemann;
    198:   ReconstructFunction characteristic;
    199:   PetscErrorCode (*destroy)(void *);
    200:   void    *user;
    201:   PetscInt dof;
    202:   char    *fieldname[16];
    203: } PhysicsCtx;

    205: typedef struct {
    206:   void (*limit)(LimitInfo, const PetscScalar *, const PetscScalar *, PetscScalar *);
    207:   PhysicsCtx physics;
    208:   MPI_Comm   comm;
    209:   char       prefix[256];

    211:   /* Local work arrays */
    212:   PetscScalar *R, *Rinv; /* Characteristic basis, and it's inverse.  COLUMN-MAJOR */
    213:   PetscScalar *cjmpLR;   /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
    214:   PetscScalar *cslope;   /* Limited slope, written in characteristic basis */
    215:   PetscScalar *uLR;      /* Solution at left and right of interface, conservative variables, len=2*dof */
    216:   PetscScalar *flux;     /* Flux across interface */
    217:   PetscReal   *speeds;   /* Speeds of each wave */

    219:   PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
    220:   PetscReal cfl;
    221:   PetscReal xmin, xmax;
    222:   PetscInt  initial;
    223:   PetscBool exact;
    224:   FVBCType  bctype;
    225: } FVCtx;

    227: PetscErrorCode RiemannListAdd(PetscFunctionList *flist, const char *name, RiemannFunction rsolve)
    228: {
    229:   PetscFunctionBeginUser;
    230:   PetscCall(PetscFunctionListAdd(flist, name, rsolve));
    231:   PetscFunctionReturn(PETSC_SUCCESS);
    232: }

    234: PetscErrorCode RiemannListFind(PetscFunctionList flist, const char *name, RiemannFunction *rsolve)
    235: {
    236:   PetscFunctionBeginUser;
    237:   PetscCall(PetscFunctionListFind(flist, name, rsolve));
    238:   PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
    239:   PetscFunctionReturn(PETSC_SUCCESS);
    240: }

    242: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist, const char *name, ReconstructFunction r)
    243: {
    244:   PetscFunctionBeginUser;
    245:   PetscCall(PetscFunctionListAdd(flist, name, r));
    246:   PetscFunctionReturn(PETSC_SUCCESS);
    247: }

    249: PetscErrorCode ReconstructListFind(PetscFunctionList flist, const char *name, ReconstructFunction *r)
    250: {
    251:   PetscFunctionBeginUser;
    252:   PetscCall(PetscFunctionListFind(flist, name, r));
    253:   PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

    257: /* --------------------------------- Physics ----------------------------------- */
    258: /*
    259:   Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction.  These
    260:   are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
    261:   number of fields and their names, and a function to deallocate private storage.
    262: */

    264: /* First a few functions useful to several different physics */
    265: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
    266: {
    267:   PetscInt i, j;

    269:   PetscFunctionBeginUser;
    270:   for (i = 0; i < m; i++) {
    271:     for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j);
    272:     speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
    273:   }
    274:   PetscFunctionReturn(PETSC_SUCCESS);
    275: }

    277: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
    278: {
    279:   PetscFunctionBeginUser;
    280:   PetscCall(PetscFree(vctx));
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }

    284: /* --------------------------------- Advection ----------------------------------- */

    286: typedef struct {
    287:   PetscReal a; /* advective velocity */
    288: } AdvectCtx;

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

    295:   PetscFunctionBeginUser;
    296:   speed     = ctx->a;
    297:   flux[0]   = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
    298:   *maxspeed = speed;
    299:   PetscFunctionReturn(PETSC_SUCCESS);
    300: }

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

    306:   PetscFunctionBeginUser;
    307:   X[0]      = 1.;
    308:   Xi[0]     = 1.;
    309:   speeds[0] = ctx->a;
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

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

    318:   PetscFunctionBeginUser;
    319:   switch (bctype) {
    320:   case FVBC_OUTFLOW:
    321:     x0 = x - a * t;
    322:     break;
    323:   case FVBC_PERIODIC:
    324:     x0 = RangeMod(x - a * t, xmin, xmax);
    325:     break;
    326:   default:
    327:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
    328:   }
    329:   switch (initial) {
    330:   case 0:
    331:     u[0] = (x0 < 0) ? 1 : -1;
    332:     break;
    333:   case 1:
    334:     u[0] = (x0 < 0) ? -1 : 1;
    335:     break;
    336:   case 2:
    337:     u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
    338:     break;
    339:   case 3:
    340:     u[0] = PetscSinReal(2 * PETSC_PI * x0);
    341:     break;
    342:   case 4:
    343:     u[0] = PetscAbs(x0);
    344:     break;
    345:   case 5:
    346:     u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    347:     break;
    348:   case 6:
    349:     u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    350:     break;
    351:   case 7:
    352:     u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    353:     break;
    354:   default:
    355:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    356:   }
    357:   PetscFunctionReturn(PETSC_SUCCESS);
    358: }

    360: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
    361: {
    362:   AdvectCtx *user;

    364:   PetscFunctionBeginUser;
    365:   PetscCall(PetscNew(&user));
    366:   ctx->physics.sample         = PhysicsSample_Advect;
    367:   ctx->physics.riemann        = PhysicsRiemann_Advect;
    368:   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
    369:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
    370:   ctx->physics.user           = user;
    371:   ctx->physics.dof            = 1;
    372:   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
    373:   user->a = 1;
    374:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    375:   {
    376:     PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
    377:   }
    378:   PetscOptionsEnd();
    379:   PetscFunctionReturn(PETSC_SUCCESS);
    380: }

    382: /* --------------------------------- Burgers ----------------------------------- */

    384: typedef struct {
    385:   PetscReal lxf_speed;
    386: } BurgersCtx;

    388: static PetscErrorCode PhysicsSample_Burgers(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    389: {
    390:   PetscFunctionBeginUser;
    391:   PetscCheck(bctype != FVBC_PERIODIC || t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solution not implemented for periodic");
    392:   switch (initial) {
    393:   case 0:
    394:     u[0] = (x < 0) ? 1 : -1;
    395:     break;
    396:   case 1:
    397:     if (x < -t) u[0] = -1;
    398:     else if (x < t) u[0] = x / t;
    399:     else u[0] = 1;
    400:     break;
    401:   case 2:
    402:     if (x <= 0) u[0] = 0;
    403:     else if (x < t) u[0] = x / t;
    404:     else if (x < 1 + 0.5 * t) u[0] = 1;
    405:     else u[0] = 0;
    406:     break;
    407:   case 3:
    408:     if (x < 0.2 * t) u[0] = 0.2;
    409:     else if (x < t) u[0] = x / t;
    410:     else u[0] = 1;
    411:     break;
    412:   case 4:
    413:     PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only initial condition available");
    414:     u[0] = 0.7 + 0.3 * PetscSinReal(2 * PETSC_PI * ((x - xmin) / (xmax - xmin)));
    415:     break;
    416:   case 5: /* Pure shock solution */
    417:     if (x < 0.5 * t) u[0] = 1;
    418:     else u[0] = 0;
    419:     break;
    420:   default:
    421:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    422:   }
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    427: {
    428:   PetscFunctionBeginUser;
    429:   if (uL[0] < uR[0]) {                /* rarefaction */
    430:     flux[0] = (uL[0] * uR[0] < 0) ? 0 /* sonic rarefaction */
    431:                                   : 0.5 * PetscMin(PetscSqr(uL[0]), PetscSqr(uR[0]));
    432:   } else { /* shock */
    433:     flux[0] = 0.5 * PetscMax(PetscSqr(uL[0]), PetscSqr(uR[0]));
    434:   }
    435:   *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
    436:   PetscFunctionReturn(PETSC_SUCCESS);
    437: }

    439: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    440: {
    441:   PetscReal speed;

    443:   PetscFunctionBeginUser;
    444:   speed   = 0.5 * (uL[0] + uR[0]);
    445:   flux[0] = 0.25 * (PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5 * PetscAbs(speed) * (uR[0] - uL[0]);
    446:   if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
    447:   *maxspeed = speed;
    448:   PetscFunctionReturn(PETSC_SUCCESS);
    449: }

    451: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    452: {
    453:   PetscReal   c;
    454:   PetscScalar fL, fR;

    456:   PetscFunctionBeginUser;
    457:   c         = ((BurgersCtx *)vctx)->lxf_speed;
    458:   fL        = 0.5 * PetscSqr(uL[0]);
    459:   fR        = 0.5 * PetscSqr(uR[0]);
    460:   flux[0]   = 0.5 * (fL + fR) - 0.5 * c * (uR[0] - uL[0]);
    461:   *maxspeed = c;
    462:   PetscFunctionReturn(PETSC_SUCCESS);
    463: }

    465: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    466: {
    467:   PetscReal   c;
    468:   PetscScalar fL, fR;

    470:   PetscFunctionBeginUser;
    471:   c         = PetscMax(PetscAbs(uL[0]), PetscAbs(uR[0]));
    472:   fL        = 0.5 * PetscSqr(uL[0]);
    473:   fR        = 0.5 * PetscSqr(uR[0]);
    474:   flux[0]   = 0.5 * (fL + fR) - 0.5 * c * (uR[0] - uL[0]);
    475:   *maxspeed = c;
    476:   PetscFunctionReturn(PETSC_SUCCESS);
    477: }

    479: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
    480: {
    481:   BurgersCtx       *user;
    482:   RiemannFunction   r;
    483:   PetscFunctionList rlist      = 0;
    484:   char              rname[256] = "exact";

    486:   PetscFunctionBeginUser;
    487:   PetscCall(PetscNew(&user));

    489:   ctx->physics.sample         = PhysicsSample_Burgers;
    490:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
    491:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
    492:   ctx->physics.user           = user;
    493:   ctx->physics.dof            = 1;

    495:   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
    496:   PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Burgers_Exact));
    497:   PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_Burgers_Roe));
    498:   PetscCall(RiemannListAdd(&rlist, "lxf", PhysicsRiemann_Burgers_LxF));
    499:   PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Burgers_Rusanov));
    500:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    501:   {
    502:     PetscCall(PetscOptionsFList("-physics_burgers_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    503:   }
    504:   PetscOptionsEnd();
    505:   PetscCall(RiemannListFind(rlist, rname, &r));
    506:   PetscCall(PetscFunctionListDestroy(&rlist));
    507:   ctx->physics.riemann = r;

    509:   /* *
    510:   * Hack to deal with LxF in semi-discrete form
    511:   * max speed is 1 for the basic initial conditions (where |u| <= 1)
    512:   * */
    513:   if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
    514:   PetscFunctionReturn(PETSC_SUCCESS);
    515: }

    517: /* --------------------------------- Traffic ----------------------------------- */

    519: typedef struct {
    520:   PetscReal lxf_speed;
    521:   PetscReal a;
    522: } TrafficCtx;

    524: static inline PetscScalar TrafficFlux(PetscScalar a, PetscScalar u)
    525: {
    526:   return a * u * (1 - u);
    527: }

    529: static PetscErrorCode PhysicsSample_Traffic(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    530: {
    531:   PetscReal a = ((TrafficCtx *)vctx)->a;

    533:   PetscFunctionBeginUser;
    534:   PetscCheck(bctype != FVBC_PERIODIC || t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solution not implemented for periodic");
    535:   switch (initial) {
    536:   case 0:
    537:     u[0] = (-a * t < x) ? 2 : 0;
    538:     break;
    539:   case 1:
    540:     if (x < PetscMin(2 * a * t, 0.5 + a * t)) u[0] = -1;
    541:     else if (x < 1) u[0] = 0;
    542:     else u[0] = 1;
    543:     break;
    544:   case 2:
    545:     PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only initial condition available");
    546:     u[0] = 0.7 + 0.3 * PetscSinReal(2 * PETSC_PI * ((x - xmin) / (xmax - xmin)));
    547:     break;
    548:   default:
    549:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    550:   }
    551:   PetscFunctionReturn(PETSC_SUCCESS);
    552: }

    554: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    555: {
    556:   PetscReal a = ((TrafficCtx *)vctx)->a;

    558:   PetscFunctionBeginUser;
    559:   if (uL[0] < uR[0]) {
    560:     flux[0] = PetscMin(TrafficFlux(a, uL[0]), TrafficFlux(a, uR[0]));
    561:   } else {
    562:     flux[0] = (uR[0] < 0.5 && 0.5 < uL[0]) ? TrafficFlux(a, 0.5) : PetscMax(TrafficFlux(a, uL[0]), TrafficFlux(a, uR[0]));
    563:   }
    564:   *maxspeed = a * MaxAbs(1 - 2 * uL[0], 1 - 2 * uR[0]);
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: }

    568: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    569: {
    570:   PetscReal a = ((TrafficCtx *)vctx)->a;
    571:   PetscReal speed;

    573:   PetscFunctionBeginUser;
    574:   speed     = a * (1 - (uL[0] + uR[0]));
    575:   flux[0]   = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * PetscAbs(speed) * (uR[0] - uL[0]);
    576:   *maxspeed = speed;
    577:   PetscFunctionReturn(PETSC_SUCCESS);
    578: }

    580: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    581: {
    582:   TrafficCtx *phys = (TrafficCtx *)vctx;
    583:   PetscReal   a    = phys->a;
    584:   PetscReal   speed;

    586:   PetscFunctionBeginUser;
    587:   speed     = a * (1 - (uL[0] + uR[0]));
    588:   flux[0]   = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * phys->lxf_speed * (uR[0] - uL[0]);
    589:   *maxspeed = speed;
    590:   PetscFunctionReturn(PETSC_SUCCESS);
    591: }

    593: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    594: {
    595:   PetscReal a = ((TrafficCtx *)vctx)->a;
    596:   PetscReal speed;

    598:   PetscFunctionBeginUser;
    599:   speed     = a * PetscMax(PetscAbs(1 - 2 * uL[0]), PetscAbs(1 - 2 * uR[0]));
    600:   flux[0]   = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * speed * (uR[0] - uL[0]);
    601:   *maxspeed = speed;
    602:   PetscFunctionReturn(PETSC_SUCCESS);
    603: }

    605: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
    606: {
    607:   TrafficCtx       *user;
    608:   RiemannFunction   r;
    609:   PetscFunctionList rlist      = 0;
    610:   char              rname[256] = "exact";

    612:   PetscFunctionBeginUser;
    613:   PetscCall(PetscNew(&user));
    614:   ctx->physics.sample         = PhysicsSample_Traffic;
    615:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
    616:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
    617:   ctx->physics.user           = user;
    618:   ctx->physics.dof            = 1;

    620:   PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
    621:   user->a = 0.5;
    622:   PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Traffic_Exact));
    623:   PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_Traffic_Roe));
    624:   PetscCall(RiemannListAdd(&rlist, "lxf", PhysicsRiemann_Traffic_LxF));
    625:   PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Traffic_Rusanov));
    626:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Traffic", "");
    627:   PetscCall(PetscOptionsReal("-physics_traffic_a", "Flux = a*u*(1-u)", "", user->a, &user->a, NULL));
    628:   PetscCall(PetscOptionsFList("-physics_traffic_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    629:   PetscOptionsEnd();

    631:   PetscCall(RiemannListFind(rlist, rname, &r));
    632:   PetscCall(PetscFunctionListDestroy(&rlist));

    634:   ctx->physics.riemann = r;

    636:   /* *
    637:   * Hack to deal with LxF in semi-discrete form
    638:   * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
    639:   * */
    640:   if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3 * user->a;
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }

    644: /* --------------------------------- Linear Acoustics ----------------------------------- */

    646: /* Flux: u_t + (A u)_x
    647:  * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
    648:  * Spectral decomposition: A = R * D * Rinv
    649:  * [    cz] = [-z   z] [-c    ] [-1/2z  1/2]
    650:  * [c/z   ] = [ 1   1] [     c] [ 1/2z  1/2]
    651:  *
    652:  * We decompose this into the left-traveling waves Al = R * D^- Rinv
    653:  * and the right-traveling waves Ar = R * D^+ * Rinv
    654:  * Multiplying out these expressions produces the following two matrices
    655:  */

    657: typedef struct {
    658:   PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
    659:   PetscReal z; /* impedance: z = sqrt(rho*bulk) */
    660: } AcousticsCtx;

    662: PETSC_UNUSED static inline void AcousticsFlux(AcousticsCtx *ctx, const PetscScalar *u, PetscScalar *f)
    663: {
    664:   f[0] = ctx->c * ctx->z * u[1];
    665:   f[1] = ctx->c / ctx->z * u[0];
    666: }

    668: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
    669: {
    670:   AcousticsCtx *phys = (AcousticsCtx *)vctx;
    671:   PetscReal     z = phys->z, c = phys->c;

    673:   PetscFunctionBeginUser;
    674:   X[0 * 2 + 0]  = -z;
    675:   X[0 * 2 + 1]  = z;
    676:   X[1 * 2 + 0]  = 1;
    677:   X[1 * 2 + 1]  = 1;
    678:   Xi[0 * 2 + 0] = -1. / (2 * z);
    679:   Xi[0 * 2 + 1] = 1. / 2;
    680:   Xi[1 * 2 + 0] = 1. / (2 * z);
    681:   Xi[1 * 2 + 1] = 1. / 2;
    682:   speeds[0]     = -c;
    683:   speeds[1]     = c;
    684:   PetscFunctionReturn(PETSC_SUCCESS);
    685: }

    687: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys, PetscInt initial, PetscReal xmin, PetscReal xmax, PetscReal x, PetscReal *u)
    688: {
    689:   PetscFunctionBeginUser;
    690:   switch (initial) {
    691:   case 0:
    692:     u[0] = (PetscAbs((x - xmin) / (xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
    693:     u[1] = (PetscAbs((x - xmin) / (xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
    694:     break;
    695:   case 1:
    696:     u[0] = PetscCosReal(3 * 2 * PETSC_PI * x / (xmax - xmin));
    697:     u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin) / 2) / (2 * PetscSqr(0.2 * (xmax - xmin)))) - 0.5;
    698:     break;
    699:   default:
    700:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    701:   }
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

    705: static PetscErrorCode PhysicsSample_Acoustics(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    706: {
    707:   AcousticsCtx *phys = (AcousticsCtx *)vctx;
    708:   PetscReal     c    = phys->c;
    709:   PetscReal     x0a, x0b, u0a[2], u0b[2], tmp[2];
    710:   PetscReal     X[2][2], Xi[2][2], dummy[2];

    712:   PetscFunctionBeginUser;
    713:   switch (bctype) {
    714:   case FVBC_OUTFLOW:
    715:     x0a = x + c * t;
    716:     x0b = x - c * t;
    717:     break;
    718:   case FVBC_PERIODIC:
    719:     x0a = RangeMod(x + c * t, xmin, xmax);
    720:     x0b = RangeMod(x - c * t, xmin, xmax);
    721:     break;
    722:   default:
    723:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
    724:   }
    725:   PetscCall(PhysicsSample_Acoustics_Initial(phys, initial, xmin, xmax, x0a, u0a));
    726:   PetscCall(PhysicsSample_Acoustics_Initial(phys, initial, xmin, xmax, x0b, u0b));
    727:   PetscCall(PhysicsCharacteristic_Acoustics(vctx, 2, u, &X[0][0], &Xi[0][0], dummy));
    728:   tmp[0] = Xi[0][0] * u0a[0] + Xi[0][1] * u0a[1];
    729:   tmp[1] = Xi[1][0] * u0b[0] + Xi[1][1] * u0b[1];
    730:   u[0]   = X[0][0] * tmp[0] + X[0][1] * tmp[1];
    731:   u[1]   = X[1][0] * tmp[0] + X[1][1] * tmp[1];
    732:   PetscFunctionReturn(PETSC_SUCCESS);
    733: }

    735: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    736: {
    737:   AcousticsCtx *phys = (AcousticsCtx *)vctx;
    738:   PetscReal     c = phys->c, z = phys->z;
    739:   PetscReal     Al[2][2] =
    740:     {
    741:       {-c / 2,      c * z / 2},
    742:       {c / (2 * z), -c / 2   }
    743:   },         /* Left traveling waves */
    744:     Ar[2][2] = {{c / 2, c * z / 2}, {c / (2 * z), c / 2}}; /* Right traveling waves */

    746:   PetscFunctionBeginUser;
    747:   flux[0]   = Al[0][0] * uR[0] + Al[0][1] * uR[1] + Ar[0][0] * uL[0] + Ar[0][1] * uL[1];
    748:   flux[1]   = Al[1][0] * uR[0] + Al[1][1] * uR[1] + Ar[1][0] * uL[0] + Ar[1][1] * uL[1];
    749:   *maxspeed = c;
    750:   PetscFunctionReturn(PETSC_SUCCESS);
    751: }

    753: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
    754: {
    755:   AcousticsCtx     *user;
    756:   PetscFunctionList rlist = 0, rclist = 0;
    757:   char              rname[256] = "exact", rcname[256] = "characteristic";

    759:   PetscFunctionBeginUser;
    760:   PetscCall(PetscNew(&user));
    761:   ctx->physics.sample  = PhysicsSample_Acoustics;
    762:   ctx->physics.destroy = PhysicsDestroy_SimpleFree;
    763:   ctx->physics.user    = user;
    764:   ctx->physics.dof     = 2;

    766:   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
    767:   PetscCall(PetscStrallocpy("v", &ctx->physics.fieldname[1]));

    769:   user->c = 1;
    770:   user->z = 1;

    772:   PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Acoustics_Exact));
    773:   PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_Acoustics));
    774:   PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
    775:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for linear Acoustics", "");
    776:   {
    777:     PetscCall(PetscOptionsReal("-physics_acoustics_c", "c = sqrt(bulk/rho)", "", user->c, &user->c, NULL));
    778:     PetscCall(PetscOptionsReal("-physics_acoustics_z", "z = sqrt(bulk*rho)", "", user->z, &user->z, NULL));
    779:     PetscCall(PetscOptionsFList("-physics_acoustics_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    780:     PetscCall(PetscOptionsFList("-physics_acoustics_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
    781:   }
    782:   PetscOptionsEnd();
    783:   PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
    784:   PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
    785:   PetscCall(PetscFunctionListDestroy(&rlist));
    786:   PetscCall(PetscFunctionListDestroy(&rclist));
    787:   PetscFunctionReturn(PETSC_SUCCESS);
    788: }

    790: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */

    792: typedef struct {
    793:   PetscReal acoustic_speed;
    794: } IsoGasCtx;

    796: static inline void IsoGasFlux(PetscReal c, const PetscScalar *u, PetscScalar *f)
    797: {
    798:   f[0] = u[1];
    799:   f[1] = PetscSqr(u[1]) / u[0] + c * c * u[0];
    800: }

    802: static PetscErrorCode PhysicsSample_IsoGas(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    803: {
    804:   PetscFunctionBeginUser;
    805:   PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0");
    806:   switch (initial) {
    807:   case 0:
    808:     u[0] = (x < 0) ? 1 : 0.5;
    809:     u[1] = (x < 0) ? 1 : 0.7;
    810:     break;
    811:   case 1:
    812:     u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x);
    813:     u[1] = 1 * u[0];
    814:     break;
    815:   default:
    816:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    817:   }
    818:   PetscFunctionReturn(PETSC_SUCCESS);
    819: }

    821: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    822: {
    823:   IsoGasCtx  *phys = (IsoGasCtx *)vctx;
    824:   PetscReal   c    = phys->acoustic_speed;
    825:   PetscScalar ubar, du[2], a[2], fL[2], fR[2], lam[2], ustar[2], R[2][2];
    826:   PetscInt    i;

    828:   PetscFunctionBeginUser;
    829:   ubar = (uL[1] / PetscSqrtScalar(uL[0]) + uR[1] / PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
    830:   /* write fluxuations in characteristic basis */
    831:   du[0] = uR[0] - uL[0];
    832:   du[1] = uR[1] - uL[1];
    833:   a[0]  = (1 / (2 * c)) * ((ubar + c) * du[0] - du[1]);
    834:   a[1]  = (1 / (2 * c)) * ((-ubar + c) * du[0] + du[1]);
    835:   /* wave speeds */
    836:   lam[0] = ubar - c;
    837:   lam[1] = ubar + c;
    838:   /* Right eigenvectors */
    839:   R[0][0] = 1;
    840:   R[0][1] = ubar - c;
    841:   R[1][0] = 1;
    842:   R[1][1] = ubar + c;
    843:   /* Compute state in star region (between the 1-wave and 2-wave) */
    844:   for (i = 0; i < 2; i++) ustar[i] = uL[i] + a[0] * R[0][i];
    845:   if (uL[1] / uL[0] < c && c < ustar[1] / ustar[0]) { /* 1-wave is sonic rarefaction */
    846:     PetscScalar ufan[2];
    847:     ufan[0] = uL[0] * PetscExpScalar(uL[1] / (uL[0] * c) - 1);
    848:     ufan[1] = c * ufan[0];
    849:     IsoGasFlux(c, ufan, flux);
    850:   } else if (ustar[1] / ustar[0] < -c && -c < uR[1] / uR[0]) { /* 2-wave is sonic rarefaction */
    851:     PetscScalar ufan[2];
    852:     ufan[0] = uR[0] * PetscExpScalar(-uR[1] / (uR[0] * c) - 1);
    853:     ufan[1] = -c * ufan[0];
    854:     IsoGasFlux(c, ufan, flux);
    855:   } else { /* Centered form */
    856:     IsoGasFlux(c, uL, fL);
    857:     IsoGasFlux(c, uR, fR);
    858:     for (i = 0; i < 2; i++) {
    859:       PetscScalar absdu = PetscAbsScalar(lam[0]) * a[0] * R[0][i] + PetscAbsScalar(lam[1]) * a[1] * R[1][i];
    860:       flux[i]           = 0.5 * (fL[i] + fR[i]) - 0.5 * absdu;
    861:     }
    862:   }
    863:   *maxspeed = MaxAbs(lam[0], lam[1]);
    864:   PetscFunctionReturn(PETSC_SUCCESS);
    865: }

    867: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    868: {
    869:   IsoGasCtx  *phys = (IsoGasCtx *)vctx;
    870:   PetscReal   c    = phys->acoustic_speed;
    871:   PetscScalar ustar[2];
    872:   struct {
    873:     PetscScalar rho, u;
    874:   } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star;
    875:   PetscInt i;

    877:   PetscFunctionBeginUser;
    878:   PetscCheck(L.rho > 0 && R.rho > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
    879:   {
    880:     /* Solve for star state */
    881:     PetscScalar res, tmp, rho = 0.5 * (L.rho + R.rho); /* initial guess */
    882:     for (i = 0; i < 20; i++) {
    883:       PetscScalar fr, fl, dfr, dfl;
    884:       fl  = (L.rho < rho) ? (rho - L.rho) / PetscSqrtScalar(L.rho * rho) /* shock */
    885:                           : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
    886:       fr  = (R.rho < rho) ? (rho - R.rho) / PetscSqrtScalar(R.rho * rho) /* shock */
    887:                           : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
    888:       res = R.u - L.u + c * (fr + fl);
    889:       PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
    890:       if (PetscAbsScalar(res) < 1e-10) {
    891:         star.rho = rho;
    892:         star.u   = L.u - c * fl;
    893:         goto converged;
    894:       }
    895:       dfl = (L.rho < rho) ? 1 / PetscSqrtScalar(L.rho * rho) * (1 - 0.5 * (rho - L.rho) / rho) : 1 / rho;
    896:       dfr = (R.rho < rho) ? 1 / PetscSqrtScalar(R.rho * rho) * (1 - 0.5 * (rho - R.rho) / rho) : 1 / rho;
    897:       tmp = rho - res / (c * (dfr + dfl));
    898:       if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
    899:       else rho = tmp;
    900:       PetscCheck((rho > 0) && PetscIsNormalScalar(rho), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate rho=%g", (double)PetscRealPart(rho));
    901:     }
    902:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Newton iteration for star.rho diverged after %" PetscInt_FMT " iterations", i);
    903:   }
    904: converged:
    905:   if (L.u - c < 0 && 0 < star.u - c) { /* 1-wave is sonic rarefaction */
    906:     PetscScalar ufan[2];
    907:     ufan[0] = L.rho * PetscExpScalar(L.u / c - 1);
    908:     ufan[1] = c * ufan[0];
    909:     IsoGasFlux(c, ufan, flux);
    910:   } else if (star.u + c < 0 && 0 < R.u + c) { /* 2-wave is sonic rarefaction */
    911:     PetscScalar ufan[2];
    912:     ufan[0] = R.rho * PetscExpScalar(-R.u / c - 1);
    913:     ufan[1] = -c * ufan[0];
    914:     IsoGasFlux(c, ufan, flux);
    915:   } else if ((L.rho >= star.rho && L.u - c >= 0) || (L.rho < star.rho && (star.rho * star.u - L.rho * L.u) / (star.rho - L.rho) > 0)) {
    916:     /* 1-wave is supersonic rarefaction, or supersonic shock */
    917:     IsoGasFlux(c, uL, flux);
    918:   } else if ((star.rho <= R.rho && R.u + c <= 0) || (star.rho > R.rho && (R.rho * R.u - star.rho * star.u) / (R.rho - star.rho) < 0)) {
    919:     /* 2-wave is supersonic rarefaction or supersonic shock */
    920:     IsoGasFlux(c, uR, flux);
    921:   } else {
    922:     ustar[0] = star.rho;
    923:     ustar[1] = star.rho * star.u;
    924:     IsoGasFlux(c, ustar, flux);
    925:   }
    926:   *maxspeed = MaxAbs(MaxAbs(star.u - c, star.u + c), MaxAbs(L.u - c, R.u + c));
    927:   PetscFunctionReturn(PETSC_SUCCESS);
    928: }

    930: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    931: {
    932:   IsoGasCtx  *phys = (IsoGasCtx *)vctx;
    933:   PetscScalar c    = phys->acoustic_speed, fL[2], fR[2], s;
    934:   struct {
    935:     PetscScalar rho, u;
    936:   } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};

    938:   PetscFunctionBeginUser;
    939:   PetscCheck(L.rho > 0 && R.rho > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
    940:   IsoGasFlux(c, uL, fL);
    941:   IsoGasFlux(c, uR, fR);
    942:   s         = PetscMax(PetscAbs(L.u), PetscAbs(R.u)) + c;
    943:   flux[0]   = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (uL[0] - uR[0]);
    944:   flux[1]   = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
    945:   *maxspeed = s;
    946:   PetscFunctionReturn(PETSC_SUCCESS);
    947: }

    949: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
    950: {
    951:   IsoGasCtx *phys = (IsoGasCtx *)vctx;
    952:   PetscReal  c    = phys->acoustic_speed;

    954:   PetscFunctionBeginUser;
    955:   speeds[0]    = u[1] / u[0] - c;
    956:   speeds[1]    = u[1] / u[0] + c;
    957:   X[0 * 2 + 0] = 1;
    958:   X[0 * 2 + 1] = speeds[0];
    959:   X[1 * 2 + 0] = 1;
    960:   X[1 * 2 + 1] = speeds[1];
    961:   PetscCall(PetscArraycpy(Xi, X, 4));
    962:   PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
    963:   PetscFunctionReturn(PETSC_SUCCESS);
    964: }

    966: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
    967: {
    968:   IsoGasCtx        *user;
    969:   PetscFunctionList rlist = 0, rclist = 0;
    970:   char              rname[256] = "exact", rcname[256] = "characteristic";

    972:   PetscFunctionBeginUser;
    973:   PetscCall(PetscNew(&user));
    974:   ctx->physics.sample  = PhysicsSample_IsoGas;
    975:   ctx->physics.destroy = PhysicsDestroy_SimpleFree;
    976:   ctx->physics.user    = user;
    977:   ctx->physics.dof     = 2;

    979:   PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
    980:   PetscCall(PetscStrallocpy("momentum", &ctx->physics.fieldname[1]));

    982:   user->acoustic_speed = 1;

    984:   PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_IsoGas_Exact));
    985:   PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_IsoGas_Roe));
    986:   PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_IsoGas_Rusanov));
    987:   PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_IsoGas));
    988:   PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
    989:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for IsoGas", "");
    990:   PetscCall(PetscOptionsReal("-physics_isogas_acoustic_speed", "Acoustic speed", "", user->acoustic_speed, &user->acoustic_speed, NULL));
    991:   PetscCall(PetscOptionsFList("-physics_isogas_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    992:   PetscCall(PetscOptionsFList("-physics_isogas_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
    993:   PetscOptionsEnd();
    994:   PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
    995:   PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
    996:   PetscCall(PetscFunctionListDestroy(&rlist));
    997:   PetscCall(PetscFunctionListDestroy(&rclist));
    998:   PetscFunctionReturn(PETSC_SUCCESS);
    999: }

    1001: /* --------------------------------- Shallow Water ----------------------------------- */
    1002: typedef struct {
    1003:   PetscReal gravity;
    1004: } ShallowCtx;

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

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

    1021:   PetscFunctionBeginUser;
    1022:   PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
    1023:   cL = PetscSqrtScalar(g * L.h);
    1024:   cR = PetscSqrtScalar(g * R.h);
    1025:   c  = PetscMax(cL, cR);
    1026:   {
    1027:     /* Solve for star state */
    1028:     const PetscInt maxits = 50;
    1029:     PetscScalar    tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */
    1030:     h0 = h;
    1031:     for (i = 0; i < maxits; i++) {
    1032:       PetscScalar fr, fl, dfr, dfl;
    1033:       fl  = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */
    1034:                       : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h);         /* rarefaction */
    1035:       fr  = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */
    1036:                       : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h);         /* rarefaction */
    1037:       res = R.u - L.u + fr + fl;
    1038:       PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
    1039:       if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h - h0) < 1e-8)) {
    1040:         star.h = h;
    1041:         star.u = L.u - fl;
    1042:         goto converged;
    1043:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
    1044:         h = 0.8 * h0 + 0.2 * h;
    1045:         continue;
    1046:       }
    1047:       /* Accept the last step and take another */
    1048:       res0 = res;
    1049:       h0   = h;
    1050:       dfl  = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h);
    1051:       dfr  = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h);
    1052:       tmp  = h - res / (dfr + dfl);
    1053:       if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
    1054:       else h = tmp;
    1055:       PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h);
    1056:     }
    1057:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i);
    1058:   }
    1059: converged:
    1060:   cstar = PetscSqrtScalar(g * star.h);
    1061:   if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */
    1062:     PetscScalar ufan[2];
    1063:     ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL);
    1064:     ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0];
    1065:     ShallowFlux(phys, ufan, flux);
    1066:   } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */
    1067:     PetscScalar ufan[2];
    1068:     ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR);
    1069:     ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0];
    1070:     ShallowFlux(phys, ufan, flux);
    1071:   } 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)) {
    1072:     /* 1-wave is right-travelling shock (supersonic) */
    1073:     ShallowFlux(phys, uL, flux);
    1074:   } 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)) {
    1075:     /* 2-wave is left-travelling shock (supersonic) */
    1076:     ShallowFlux(phys, uR, flux);
    1077:   } else {
    1078:     ustar[0] = star.h;
    1079:     ustar[1] = star.h * star.u;
    1080:     ShallowFlux(phys, ustar, flux);
    1081:   }
    1082:   *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR));
    1083:   PetscFunctionReturn(PETSC_SUCCESS);
    1084: }

    1086: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
    1087: {
    1088:   ShallowCtx *phys = (ShallowCtx *)vctx;
    1089:   PetscScalar g    = phys->gravity, fL[2], fR[2], s;
    1090:   struct {
    1091:     PetscScalar h, u;
    1092:   } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};

    1094:   PetscFunctionBeginUser;
    1095:   PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
    1096:   ShallowFlux(phys, uL, fL);
    1097:   ShallowFlux(phys, uR, fR);
    1098:   s         = PetscMax(PetscAbs(L.u) + PetscSqrtScalar(g * L.h), PetscAbs(R.u) + PetscSqrtScalar(g * R.h));
    1099:   flux[0]   = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (uL[0] - uR[0]);
    1100:   flux[1]   = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
    1101:   *maxspeed = s;
    1102:   PetscFunctionReturn(PETSC_SUCCESS);
    1103: }

    1105: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
    1106: {
    1107:   ShallowCtx *phys = (ShallowCtx *)vctx;
    1108:   PetscReal   c;

    1110:   PetscFunctionBeginUser;
    1111:   c            = PetscSqrtScalar(u[0] * phys->gravity);
    1112:   speeds[0]    = u[1] / u[0] - c;
    1113:   speeds[1]    = u[1] / u[0] + c;
    1114:   X[0 * 2 + 0] = 1;
    1115:   X[0 * 2 + 1] = speeds[0];
    1116:   X[1 * 2 + 0] = 1;
    1117:   X[1 * 2 + 1] = speeds[1];
    1118:   PetscCall(PetscArraycpy(Xi, X, 4));
    1119:   PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
    1120:   PetscFunctionReturn(PETSC_SUCCESS);
    1121: }

    1123: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
    1124: {
    1125:   ShallowCtx       *user;
    1126:   PetscFunctionList rlist = 0, rclist = 0;
    1127:   char              rname[256] = "exact", rcname[256] = "characteristic";

    1129:   PetscFunctionBeginUser;
    1130:   PetscCall(PetscNew(&user));
    1131:   /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
    1132:   ctx->physics.sample  = PhysicsSample_IsoGas;
    1133:   ctx->physics.destroy = PhysicsDestroy_SimpleFree;
    1134:   ctx->physics.user    = user;
    1135:   ctx->physics.dof     = 2;

    1137:   PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
    1138:   PetscCall(PetscStrallocpy("momentum", &ctx->physics.fieldname[1]));

    1140:   user->gravity = 1;

    1142:   PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Shallow_Exact));
    1143:   PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov));
    1144:   PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_Shallow));
    1145:   PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
    1146:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", "");
    1147:   PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL));
    1148:   PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    1149:   PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
    1150:   PetscOptionsEnd();
    1151:   PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
    1152:   PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
    1153:   PetscCall(PetscFunctionListDestroy(&rlist));
    1154:   PetscCall(PetscFunctionListDestroy(&rclist));
    1155:   PetscFunctionReturn(PETSC_SUCCESS);
    1156: }

    1158: /* --------------------------------- Finite Volume Solver ----------------------------------- */

    1160: static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    1161: {
    1162:   FVCtx       *ctx = (FVCtx *)vctx;
    1163:   PetscInt     i, j, k, Mx, dof, xs, xm;
    1164:   PetscReal    hx, cfl_idt = 0;
    1165:   PetscScalar *x, *f, *slope;
    1166:   Vec          Xloc;
    1167:   DM           da;

    1169:   PetscFunctionBeginUser;
    1170:   PetscCall(TSGetDM(ts, &da));
    1171:   PetscCall(DMGetLocalVector(da, &Xloc));
    1172:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    1173:   hx = (ctx->xmax - ctx->xmin) / Mx;
    1174:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    1175:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

    1177:   PetscCall(VecZeroEntries(F));

    1179:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    1180:   PetscCall(DMDAVecGetArray(da, F, &f));
    1181:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));

    1183:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    1185:   if (ctx->bctype == FVBC_OUTFLOW) {
    1186:     for (i = xs - 2; i < 0; i++) {
    1187:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    1188:     }
    1189:     for (i = Mx; i < xs + xm + 2; i++) {
    1190:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    1191:     }
    1192:   }
    1193:   for (i = xs - 1; i < xs + xm + 1; i++) {
    1194:     struct _LimitInfo info;
    1195:     PetscScalar      *cjmpL, *cjmpR;
    1196:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    1197:     PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    1198:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    1199:     PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    1200:     cjmpL = &ctx->cjmpLR[0];
    1201:     cjmpR = &ctx->cjmpLR[dof];
    1202:     for (j = 0; j < dof; j++) {
    1203:       PetscScalar jmpL, jmpR;
    1204:       jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    1205:       jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    1206:       for (k = 0; k < dof; k++) {
    1207:         cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    1208:         cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    1209:       }
    1210:     }
    1211:     /* Apply limiter to the left and right characteristic jumps */
    1212:     info.m  = dof;
    1213:     info.hx = hx;
    1214:     (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    1215:     for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    1216:     for (j = 0; j < dof; j++) {
    1217:       PetscScalar tmp = 0;
    1218:       for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    1219:       slope[i * dof + j] = tmp;
    1220:     }
    1221:   }

    1223:   for (i = xs; i < xs + xm + 1; i++) {
    1224:     PetscReal    maxspeed;
    1225:     PetscScalar *uL, *uR;
    1226:     uL = &ctx->uLR[0];
    1227:     uR = &ctx->uLR[dof];
    1228:     for (j = 0; j < dof; j++) {
    1229:       uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    1230:       uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    1231:     }
    1232:     PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed));
    1233:     cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */

    1235:     if (i > xs) {
    1236:       for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
    1237:     }
    1238:     if (i < xs + xm) {
    1239:       for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
    1240:     }
    1241:   }

    1243:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    1244:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    1245:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    1246:   PetscCall(DMRestoreLocalVector(da, &Xloc));

    1248:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    1249:   if (0) {
    1250:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
    1251:     PetscReal dt, tnow;
    1252:     PetscCall(TSGetTimeStep(ts, &dt));
    1253:     PetscCall(TSGetTime(ts, &tnow));
    1254:     if (dt > 0.5 / ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
    1255:   }
    1256:   PetscFunctionReturn(PETSC_SUCCESS);
    1257: }

    1259: static PetscErrorCode SmallMatMultADB(PetscScalar *C, PetscInt bs, const PetscScalar *A, const PetscReal *D, const PetscScalar *B)
    1260: {
    1261:   PetscInt i, j, k;

    1263:   PetscFunctionBeginUser;
    1264:   for (i = 0; i < bs; i++) {
    1265:     for (j = 0; j < bs; j++) {
    1266:       PetscScalar tmp = 0;
    1267:       for (k = 0; k < bs; k++) tmp += A[i * bs + k] * D[k] * B[k * bs + j];
    1268:       C[i * bs + j] = tmp;
    1269:     }
    1270:   }
    1271:   PetscFunctionReturn(PETSC_SUCCESS);
    1272: }

    1274: static PetscErrorCode FVIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat A, Mat B, void *vctx)
    1275: {
    1276:   FVCtx             *ctx = (FVCtx *)vctx;
    1277:   PetscInt           i, j, dof = ctx->physics.dof;
    1278:   PetscScalar       *J;
    1279:   const PetscScalar *x;
    1280:   PetscReal          hx;
    1281:   DM                 da;
    1282:   DMDALocalInfo      dainfo;

    1284:   PetscFunctionBeginUser;
    1285:   PetscCall(TSGetDM(ts, &da));
    1286:   PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
    1287:   PetscCall(DMDAGetLocalInfo(da, &dainfo));
    1288:   hx = (ctx->xmax - ctx->xmin) / dainfo.mx;
    1289:   PetscCall(PetscMalloc1(dof * dof, &J));
    1290:   for (i = dainfo.xs; i < dainfo.xs + dainfo.xm; i++) {
    1291:     PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    1292:     for (j = 0; j < dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
    1293:     PetscCall(SmallMatMultADB(J, dof, ctx->R, ctx->speeds, ctx->Rinv));
    1294:     for (j = 0; j < dof * dof; j++) J[j] = J[j] / hx + shift * (j / dof == j % dof);
    1295:     PetscCall(MatSetValuesBlocked(B, 1, &i, 1, &i, J, INSERT_VALUES));
    1296:   }
    1297:   PetscCall(PetscFree(J));
    1298:   PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));

    1300:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    1301:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    1302:   if (A != B) {
    1303:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1304:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    1305:   }
    1306:   PetscFunctionReturn(PETSC_SUCCESS);
    1307: }

    1309: static PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
    1310: {
    1311:   PetscScalar *u, *uj;
    1312:   PetscInt     i, j, k, dof, xs, xm, Mx;

    1314:   PetscFunctionBeginUser;
    1315:   PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
    1316:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    1317:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    1318:   PetscCall(DMDAVecGetArray(da, U, &u));
    1319:   PetscCall(PetscMalloc1(dof, &uj));
    1320:   for (i = xs; i < xs + xm; i++) {
    1321:     const PetscReal h = (ctx->xmax - ctx->xmin) / Mx, xi = ctx->xmin + h / 2 + i * h;
    1322:     const PetscInt  N = 200;
    1323:     /* Integrate over cell i using trapezoid rule with N points. */
    1324:     for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    1325:     for (j = 0; j < N + 1; j++) {
    1326:       PetscScalar xj = xi + h * (j - N / 2) / (PetscReal)N;
    1327:       PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    1328:       for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    1329:     }
    1330:   }
    1331:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    1332:   PetscCall(PetscFree(uj));
    1333:   PetscFunctionReturn(PETSC_SUCCESS);
    1334: }

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

    1345:   PetscFunctionBeginUser;
    1346:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    1347:   if (iascii) {
    1348:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
    1349:     PetscCall(DMGetLocalVector(da, &Xloc));
    1350:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    1351:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    1352:     PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
    1353:     PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    1354:     PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    1355:     tvsum = 0;
    1356:     for (i = xs; i < xs + xm; i++) {
    1357:       for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
    1358:     }
    1359:     PetscCallMPI(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)da)));
    1360:     PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
    1361:     PetscCall(DMRestoreLocalVector(da, &Xloc));

    1363:     PetscCall(VecMin(X, &imin, &xmin));
    1364:     PetscCall(VecMax(X, &imax, &xmax));
    1365:     PetscCall(VecSum(X, &sum));
    1366:     PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%8.5f,%8.5f] with extrema at %" PetscInt_FMT " and %" PetscInt_FMT ", mean %8.5f, ||x||_TV %8.5f\n", (double)xmin, (double)xmax, imin, imax, (double)(sum / Mx), (double)(tvgsum / Mx)));
    1367:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
    1368:   PetscFunctionReturn(PETSC_SUCCESS);
    1369: }

    1371: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1, PetscReal *nrmsup)
    1372: {
    1373:   Vec      Y;
    1374:   PetscInt Mx;

    1376:   PetscFunctionBeginUser;
    1377:   PetscCall(VecGetSize(X, &Mx));
    1378:   PetscCall(VecDuplicate(X, &Y));
    1379:   PetscCall(FVSample(ctx, da, t, Y));
    1380:   PetscCall(VecAYPX(Y, -1, X));
    1381:   PetscCall(VecNorm(Y, NORM_1, nrm1));
    1382:   PetscCall(VecNorm(Y, NORM_INFINITY, nrmsup));
    1383:   *nrm1 /= Mx;
    1384:   PetscCall(VecDestroy(&Y));
    1385:   PetscFunctionReturn(PETSC_SUCCESS);
    1386: }

    1388: int main(int argc, char *argv[])
    1389: {
    1390:   char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
    1391:   PetscFunctionList limiters = 0, physics = 0;
    1392:   MPI_Comm          comm;
    1393:   TS                ts;
    1394:   DM                da;
    1395:   Vec               X, X0, R;
    1396:   Mat               B;
    1397:   FVCtx             ctx;
    1398:   PetscInt          i, dof, xs, xm, Mx, draw = 0;
    1399:   PetscBool         view_final = PETSC_FALSE;
    1400:   PetscReal         ptime;

    1402:   PetscFunctionBeginUser;
    1403:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
    1404:   comm = PETSC_COMM_WORLD;
    1405:   PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

    1407:   /* Register limiters to be available on the command line */
    1408:   PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind));
    1409:   PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff));
    1410:   PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming));
    1411:   PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm));
    1412:   PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod));
    1413:   PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee));
    1414:   PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit_MC));
    1415:   PetscCall(PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer));
    1416:   PetscCall(PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada));
    1417:   PetscCall(PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD));
    1418:   PetscCall(PetscFunctionListAdd(&limiters, "koren", Limit_Koren));
    1419:   PetscCall(PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym));
    1420:   PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3));
    1421:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2));
    1422:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1));
    1423:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1));
    1424:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10));
    1425:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100));

    1427:   /* Register physical models to be available on the command line */
    1428:   PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
    1429:   PetscCall(PetscFunctionListAdd(&physics, "burgers", PhysicsCreate_Burgers));
    1430:   PetscCall(PetscFunctionListAdd(&physics, "traffic", PhysicsCreate_Traffic));
    1431:   PetscCall(PetscFunctionListAdd(&physics, "acoustics", PhysicsCreate_Acoustics));
    1432:   PetscCall(PetscFunctionListAdd(&physics, "isogas", PhysicsCreate_IsoGas));
    1433:   PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow));

    1435:   ctx.comm   = comm;
    1436:   ctx.cfl    = 0.9;
    1437:   ctx.bctype = FVBC_PERIODIC;
    1438:   ctx.xmin   = -1;
    1439:   ctx.xmax   = 1;
    1440:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
    1441:   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
    1442:   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
    1443:   PetscCall(PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
    1444:   PetscCall(PetscOptionsFList("-physics", "Name of physics (Riemann solver and characteristics) to use", "", physics, physname, physname, sizeof(physname), NULL));
    1445:   PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
    1446:   PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
    1447:   PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
    1448:   PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
    1449:   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
    1450:   PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
    1451:   PetscOptionsEnd();

    1453:   /* Choose the limiter from the list of registered limiters */
    1454:   PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit));
    1455:   PetscCheck(ctx.limit, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);

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

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

    1476:   /* Set coordinates of cell centers */
    1477:   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));

    1479:   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
    1480:   PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
    1481:   PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));

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

    1488:   PetscCall(DMCreateMatrix(da, &B));

    1490:   /* Create a time-stepping object */
    1491:   PetscCall(TSCreate(comm, &ts));
    1492:   PetscCall(TSSetDM(ts, da));
    1493:   PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
    1494:   PetscCall(TSSetIJacobian(ts, B, B, FVIJacobian, &ctx));
    1495:   PetscCall(TSSetType(ts, TSSSP));
    1496:   PetscCall(TSSetMaxTime(ts, 10));
    1497:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    1499:   /* Compute initial conditions and starting time step */
    1500:   PetscCall(FVSample(&ctx, da, 0, X0));
    1501:   PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
    1502:   PetscCall(VecCopy(X0, X));                            /* The function value was not used so we set X=X0 again */
    1503:   PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
    1504:   PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
    1505:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    1506:   {
    1507:     PetscReal nrm1, nrmsup;
    1508:     PetscInt  steps;

    1510:     PetscCall(TSSolve(ts, X));
    1511:     PetscCall(TSGetSolveTime(ts, &ptime));
    1512:     PetscCall(TSGetStepNumber(ts, &steps));

    1514:     PetscCall(PetscPrintf(comm, "Final time %8.5f, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    1515:     if (ctx.exact) {
    1516:       PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1, &nrmsup));
    1517:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %8.4e  ||x-x_e||_sup %8.4e\n", (double)nrm1, (double)nrmsup));
    1518:     }
    1519:   }

    1521:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    1522:   if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
    1523:   if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
    1524:   if (draw & 0x4) {
    1525:     Vec Y;
    1526:     PetscCall(VecDuplicate(X, &Y));
    1527:     PetscCall(FVSample(&ctx, da, ptime, Y));
    1528:     PetscCall(VecAYPX(Y, -1, X));
    1529:     PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    1530:     PetscCall(VecDestroy(&Y));
    1531:   }

    1533:   if (view_final) {
    1534:     PetscViewer viewer;
    1535:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    1536:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    1537:     PetscCall(VecView(X, viewer));
    1538:     PetscCall(PetscViewerPopFormat(viewer));
    1539:     PetscCall(PetscViewerDestroy(&viewer));
    1540:   }

    1542:   /* Clean up */
    1543:   PetscCall((*ctx.physics.destroy)(ctx.physics.user));
    1544:   for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
    1545:   PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
    1546:   PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
    1547:   PetscCall(VecDestroy(&X));
    1548:   PetscCall(VecDestroy(&X0));
    1549:   PetscCall(VecDestroy(&R));
    1550:   PetscCall(MatDestroy(&B));
    1551:   PetscCall(DMDestroy(&da));
    1552:   PetscCall(TSDestroy(&ts));
    1553:   PetscCall(PetscFunctionListDestroy(&limiters));
    1554:   PetscCall(PetscFunctionListDestroy(&physics));
    1555:   PetscCall(PetscFinalize());
    1556:   return 0;
    1557: }

    1559: /*TEST

    1561:     build:
    1562:       requires: !complex

    1564:     test:
    1565:       args: -da_grid_x 100 -initial 1 -xmin -2 -xmax 5 -exact -limit mc
    1566:       requires: !complex !single

    1568:     test:
    1569:       suffix: 2
    1570:       args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
    1571:       filter:  sed "s/at 48/at 0/g"
    1572:       requires: !complex !single

    1574:     test:
    1575:       suffix: 3
    1576:       args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
    1577:       nsize: 3
    1578:       filter:  sed "s/at 48/at 0/g"
    1579:       requires: !complex !single

    1581: TEST*/