Actual source code: finitevolume1d.c

      1: #include "finitevolume1d.h"
      2: #include <petscdmda.h>
      3: #include <petscdraw.h>
      4: #include <petsc/private/tsimpl.h>

      6: #include <petsc/private/kernels/blockinvert.h>
      7: const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "INFLOW", "FVBCType", "FVBC_", 0};

      9: static inline PetscReal Sgn(PetscReal a)
     10: {
     11:   return (a < 0) ? -1 : 1;
     12: }
     13: static inline PetscReal Abs(PetscReal a)
     14: {
     15:   return (a < 0) ? 0 : a;
     16: }
     17: static inline PetscReal Sqr(PetscReal a)
     18: {
     19:   return a * a;
     20: }

     22: PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a, PetscReal b)
     23: {
     24:   return (PetscAbs(a) < PetscAbs(b)) ? a : b;
     25: }
     26: static inline PetscReal MinMod2(PetscReal a, PetscReal b)
     27: {
     28:   return (a * b < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscAbs(b));
     29: }
     30: static inline PetscReal MaxMod2(PetscReal a, PetscReal b)
     31: {
     32:   return (a * b < 0) ? 0 : Sgn(a) * PetscMax(PetscAbs(a), PetscAbs(b));
     33: }
     34: static inline PetscReal MinMod3(PetscReal a, PetscReal b, PetscReal c)
     35: {
     36:   return (a * b < 0 || a * c < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscMin(PetscAbs(b), PetscAbs(c)));
     37: }

     39: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
     40: void Limit_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     41: {
     42:   PetscInt i;
     43:   for (i = 0; i < info->m; i++) lmt[i] = 0;
     44: }
     45: void Limit_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     46: {
     47:   PetscInt i;
     48:   for (i = 0; i < info->m; i++) lmt[i] = jR[i];
     49: }
     50: void Limit_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     51: {
     52:   PetscInt i;
     53:   for (i = 0; i < info->m; i++) lmt[i] = jL[i];
     54: }
     55: void Limit_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     56: {
     57:   PetscInt i;
     58:   for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]);
     59: }
     60: void Limit_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     61: {
     62:   PetscInt i;
     63:   for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]);
     64: }
     65: void Limit_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     66: {
     67:   PetscInt i;
     68:   for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i]));
     69: }
     70: void Limit_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     71: {
     72:   PetscInt i;
     73:   for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]);
     74: }
     75: void Limit_VanLeer(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     76: { /* phi = (t + abs(t)) / (1 + abs(t)) */
     77:   PetscInt i;
     78:   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);
     79: }
     80: void Limit_VanAlbada(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
     81: {                                                                                                    /* phi = (t + t^2) / (1 + t^2) */
     82:   PetscInt i;
     83:   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);
     84: }
     85: void Limit_VanAlbadaTVD(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
     86: { /* phi = (t + t^2) / (1 + t^2) */
     87:   PetscInt i;
     88:   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);
     89: }
     90: void Limit_Koren(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
     91: {                                                                                                /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
     92:   PetscInt i;
     93:   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));
     94: }
     95: void Limit_KorenSym(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
     96: {                                                                                                   /* Symmetric version of above */
     97:   PetscInt i;
     98:   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));
     99: }
    100: void Limit_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    101: { /* Eq 11 of Cada-Torrilhon 2009 */
    102:   PetscInt i;
    103:   for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]);
    104: }
    105: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L, PetscReal R)
    106: {
    107:   return PetscMax(0, PetscMin((L + 2 * R) / 3, PetscMax(-0.5 * L, PetscMin(2 * L, PetscMin((L + 2 * R) / 3, 1.6 * R)))));
    108: }
    109: void Limit_CadaTorrilhon2(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    110: { /* Cada-Torrilhon 2009, Eq 13 */
    111:   PetscInt i;
    112:   for (i = 0; i < info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]);
    113: }
    114: void Limit_CadaTorrilhon3R(PetscReal r, LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    115: { /* Cada-Torrilhon 2009, Eq 22 */
    116:   /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
    117:   const PetscReal eps = 1e-7, hx = info->hx;
    118:   PetscInt        i;
    119:   for (i = 0; i < info->m; i++) {
    120:     const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r * hx);
    121:     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]))));
    122:   }
    123: }
    124: void Limit_CadaTorrilhon3R0p1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    125: {
    126:   Limit_CadaTorrilhon3R(0.1, info, jL, jR, lmt);
    127: }
    128: void Limit_CadaTorrilhon3R1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    129: {
    130:   Limit_CadaTorrilhon3R(1, info, jL, jR, lmt);
    131: }
    132: void Limit_CadaTorrilhon3R10(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    133: {
    134:   Limit_CadaTorrilhon3R(10, info, jL, jR, lmt);
    135: }
    136: void Limit_CadaTorrilhon3R100(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
    137: {
    138:   Limit_CadaTorrilhon3R(100, info, jL, jR, lmt);
    139: }

    141: /* ----------------------- Limiters for split systems ------------------------- */

    143: void Limit2_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    144: {
    145:   PetscInt i;
    146:   for (i = 0; i < info->m; i++) lmt[i] = 0;
    147: }
    148: void Limit2_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    149: {
    150:   PetscInt i;
    151:   if (n < sf - 1) { /* slow components */
    152:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
    153:   } else if (n == sf - 1) { /* slow component which is next to fast components */
    154:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxs / 2.0 + info->hxf / 2.0);
    155:   } else if (n == sf) { /* fast component which is next to slow components */
    156:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
    157:   } else if (n > sf && n < fs - 1) { /* fast components */
    158:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
    159:   } else if (n == fs - 1) { /* fast component next to slow components */
    160:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxf / 2.0 + info->hxs / 2.0);
    161:   } else if (n == fs) { /* slow component next to fast components */
    162:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
    163:   } else { /* slow components */
    164:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
    165:   }
    166: }
    167: void Limit2_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    168: {
    169:   PetscInt i;
    170:   if (n < sf - 1) {
    171:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
    172:   } else if (n == sf - 1) {
    173:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
    174:   } else if (n == sf) {
    175:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxs / 2.0 + info->hxf / 2.0);
    176:   } else if (n > sf && n < fs - 1) {
    177:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
    178:   } else if (n == fs - 1) {
    179:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
    180:   } else if (n == fs) {
    181:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxf / 2.0 + info->hxs / 2.0);
    182:   } else {
    183:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
    184:   }
    185: }
    186: void Limit2_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    187: {
    188:   PetscInt i;
    189:   if (n < sf - 1) {
    190:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
    191:   } else if (n == sf - 1) {
    192:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
    193:   } else if (n == sf) {
    194:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf);
    195:   } else if (n > sf && n < fs - 1) {
    196:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
    197:   } else if (n == fs - 1) {
    198:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    199:   } else if (n == fs) {
    200:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs);
    201:   } else {
    202:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
    203:   }
    204: }
    205: void Limit2_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    206: {
    207:   PetscInt i;
    208:   if (n < sf - 1) {
    209:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
    210:   } else if (n == sf - 1) {
    211:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
    212:   } else if (n == sf) {
    213:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
    214:   } else if (n > sf && n < fs - 1) {
    215:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxf;
    216:   } else if (n == fs - 1) {
    217:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    218:   } else if (n == fs) {
    219:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), jR[i] / info->hxs);
    220:   } else {
    221:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
    222:   }
    223: }
    224: void Limit2_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    225: {
    226:   PetscInt i;
    227:   if (n < sf - 1) {
    228:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
    229:   } else if (n == sf - 1) {
    230:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxs, 2 * jR[i] / (info->hxs / 2.0 + info->hxf / 2.0)), MinMod2(2 * jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0)));
    231:   } else if (n == sf) {
    232:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), 2 * jR[i] / info->hxf), MinMod2(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf));
    233:   } else if (n > sf && n < fs - 1) {
    234:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxf;
    235:   } else if (n == fs - 1) {
    236:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxf, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), MinMod2(2 * jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)));
    237:   } else if (n == fs) {
    238:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), 2 * jR[i] / info->hxs), MinMod2(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), jR[i] / info->hxs));
    239:   } else {
    240:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
    241:   }
    242: }
    243: void Limit2_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    244: {
    245:   PetscInt i;
    246:   if (n < sf - 1) {
    247:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
    248:   } else if (n == sf - 1) {
    249:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    250:   } else if (n == sf) {
    251:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf), 2 * jR[i] / info->hxf);
    252:   } else if (n > sf && n < fs - 1) {
    253:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxf;
    254:   } else if (n == fs - 1) {
    255:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    256:   } else if (n == fs) {
    257:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs), 2 * jR[i] / info->hxs);
    258:   } else {
    259:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
    260:   }
    261: }
    262: void Limit2_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
    263: { /* Eq 11 of Cada-Torrilhon 2009 */
    264:   PetscInt i;
    265:   if (n < sf - 1) {
    266:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
    267:   } else if (n == sf - 1) {
    268:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, (jL[i] / info->hxs + 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    269:   } else if (n == sf) {
    270:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + 2 * jR[i] / info->hxf) / 3, 2 * jR[i] / info->hxf);
    271:   } else if (n > sf && n < fs - 1) {
    272:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxf;
    273:   } else if (n == fs - 1) {
    274:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, (jL[i] / info->hxf + 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
    275:   } else if (n == fs) {
    276:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + 2 * jR[i] / info->hxs) / 3, 2 * jR[i] / info->hxs);
    277:   } else {
    278:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
    279:   }
    280: }

    282: /* ---- Three-way splitting ---- */
    283: void Limit3_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    284: {
    285:   PetscInt i;
    286:   for (i = 0; i < info->m; i++) lmt[i] = 0;
    287: }
    288: void Limit3_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    289: {
    290:   PetscInt i;
    291:   if (n < sm - 1 || n > ms) { /* slow components */
    292:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
    293:   } else if (n == sm - 1 || n == ms - 1) { /* slow component which is next to medium components */
    294:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxs / 2.0 + info->hxm / 2.0);
    295:   } else if (n < mf - 1 || n > fm) { /* medium components */
    296:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxm;
    297:   } else if (n == mf - 1 || n == fm) { /* medium component next to fast components */
    298:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxm / 2.0 + info->hxf / 2.0);
    299:   } else { /* fast components */
    300:     for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
    301:   }
    302: }
    303: void Limit3_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    304: {
    305:   PetscInt i;
    306:   if (n < sm || n > ms) {
    307:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
    308:   } else if (n == sm || n == ms) {
    309:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxs / 2.0 + info->hxf / 2.0);
    310:   } else if (n < mf || n > fm) {
    311:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxm;
    312:   } else if (n == mf || n == fm) {
    313:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxm / 2.0 + info->hxf / 2.0);
    314:   } else {
    315:     for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
    316:   }
    317: }
    318: void Limit3_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    319: {
    320:   PetscInt i;
    321:   if (n < sm - 1 || n > ms) {
    322:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
    323:   } else if (n == sm - 1) {
    324:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
    325:   } else if (n == sm) {
    326:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm);
    327:   } else if (n == ms - 1) {
    328:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
    329:   } else if (n == ms) {
    330:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs);
    331:   } else if (n < mf - 1 || n > fm) {
    332:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxm;
    333:   } else if (n == mf - 1) {
    334:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
    335:   } else if (n == mf) {
    336:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf);
    337:   } else if (n == fm - 1) {
    338:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
    339:   } else if (n == fm) {
    340:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm);
    341:   } else {
    342:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
    343:   }
    344: }
    345: void Limit3_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    346: {
    347:   PetscInt i;
    348:   if (n < sm - 1 || n > ms) {
    349:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
    350:   } else if (n == sm - 1) {
    351:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
    352:   } else if (n == sm) {
    353:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
    354:   } else if (n == ms - 1) {
    355:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
    356:   } else if (n == ms) {
    357:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), jR[i] / info->hxs);
    358:   } else if (n < mf - 1 || n > fm) {
    359:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxm;
    360:   } else if (n == mf - 1) {
    361:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
    362:   } else if (n == mf) {
    363:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
    364:   } else if (n == fm - 1) {
    365:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
    366:   } else if (n == fm) {
    367:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), jR[i] / info->hxm);
    368:   } else {
    369:     for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
    370:   }
    371: }
    372: void Limit3_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    373: {
    374:   PetscInt i;
    375:   if (n < sm - 1 || n > ms) {
    376:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
    377:   } else if (n == sm - 1) {
    378:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxs, 2 * jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)), MinMod2(2 * jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)));
    379:   } else if (n == sm) {
    380:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), 2 * jR[i] / info->hxm), MinMod2(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), jR[i] / info->hxm));
    381:   } else if (n == ms - 1) {
    382:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxm, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)), MinMod2(2 * jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)));
    383:   } else if (n == ms) {
    384:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), 2 * jR[i] / info->hxs), MinMod2(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), jR[i] / info->hxs));
    385:   } else if (n < mf - 1 || n > fm) {
    386:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxm;
    387:   } else if (n == mf - 1) {
    388:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxm, 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)), MinMod2(2 * jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)));
    389:   } else if (n == mf) {
    390:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), 2 * jR[i] / info->hxf), MinMod2(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), jR[i] / info->hxf));
    391:   } else if (n == fm - 1) {
    392:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxf, 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)), MinMod2(2 * jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)));
    393:   } else if (n == fm) {
    394:     for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), 2 * jR[i] / info->hxm), MinMod2(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), jR[i] / info->hxm));
    395:   } else {
    396:     for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxf;
    397:   }
    398: }
    399: void Limit3_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    400: {
    401:   PetscInt i;
    402:   if (n < sm - 1 || n > ms) {
    403:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
    404:   } else if (n == sm - 1) {
    405:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)), 2 * jR[i] / (info->hxs / 2.0 + info->hxm / 2.0));
    406:   } else if (n == sm) {
    407:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm), 2 * jR[i] / info->hxm);
    408:   } else if (n == ms - 1) {
    409:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
    410:   } else if (n == ms) {
    411:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs), 2 * jR[i] / info->hxs);
    412:   } else if (n < mf - 1 || n > fm) {
    413:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxm;
    414:   } else if (n == mf - 1) {
    415:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)), 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
    416:   } else if (n == mf) {
    417:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf), 2 * jR[i] / info->hxf);
    418:   } else if (n == fm - 1) {
    419:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
    420:   } else if (n == fm) {
    421:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm), 2 * jR[i] / info->hxm);
    422:   } else {
    423:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxf;
    424:   }
    425: }
    426: void Limit3_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
    427: { /* Eq 11 of Cada-Torrilhon 2009 */
    428:   PetscInt i;
    429:   if (n < sm - 1 || n > ms) {
    430:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
    431:   } else if (n == sm - 1) {
    432:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, (jL[i] / info->hxs + 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
    433:   } else if (n == sm) {
    434:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + 2 * jR[i] / info->hxm) / 3, 2 * jR[i] / info->hxm);
    435:   } else if (n == ms - 1) {
    436:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, (jL[i] / info->hxm + 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
    437:   } else if (n == ms) {
    438:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + 2 * jR[i] / info->hxs) / 3, 2 * jR[i] / info->hxs);
    439:   } else if (n < mf - 1 || n > fm) {
    440:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxm;
    441:   } else if (n == mf - 1) {
    442:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, (jL[i] / info->hxm + 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
    443:   } else if (n == mf) {
    444:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + 2 * jR[i] / info->hxf) / 3, 2 * jR[i] / info->hxf);
    445:   } else if (n == fm - 1) {
    446:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, (jL[i] / info->hxf + 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
    447:   } else if (n == fm) {
    448:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + 2 * jR[i] / info->hxm) / 3, 2 * jR[i] / info->hxm);
    449:   } else {
    450:     for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
    451:   }
    452: }

    454: PetscErrorCode RiemannListAdd(PetscFunctionList *flist, const char *name, RiemannFunction rsolve)
    455: {
    456:   PetscFunctionBeginUser;
    457:   PetscCall(PetscFunctionListAdd(flist, name, rsolve));
    458:   PetscFunctionReturn(PETSC_SUCCESS);
    459: }

    461: PetscErrorCode RiemannListFind(PetscFunctionList flist, const char *name, RiemannFunction *rsolve)
    462: {
    463:   PetscFunctionBeginUser;
    464:   PetscCall(PetscFunctionListFind(flist, name, rsolve));
    465:   PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist, const char *name, ReconstructFunction r)
    470: {
    471:   PetscFunctionBeginUser;
    472:   PetscCall(PetscFunctionListAdd(flist, name, r));
    473:   PetscFunctionReturn(PETSC_SUCCESS);
    474: }

    476: PetscErrorCode ReconstructListFind(PetscFunctionList flist, const char *name, ReconstructFunction *r)
    477: {
    478:   PetscFunctionBeginUser;
    479:   PetscCall(PetscFunctionListFind(flist, name, r));
    480:   PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
    481:   PetscFunctionReturn(PETSC_SUCCESS);
    482: }

    484: PetscErrorCode RiemannListAdd_2WaySplit(PetscFunctionList *flist, const char *name, RiemannFunction_2WaySplit rsolve)
    485: {
    486:   PetscFunctionBeginUser;
    487:   PetscCall(PetscFunctionListAdd(flist, name, rsolve));
    488:   PetscFunctionReturn(PETSC_SUCCESS);
    489: }

    491: PetscErrorCode RiemannListFind_2WaySplit(PetscFunctionList flist, const char *name, RiemannFunction_2WaySplit *rsolve)
    492: {
    493:   PetscFunctionBeginUser;
    494:   PetscCall(PetscFunctionListFind(flist, name, rsolve));
    495:   PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }

    499: PetscErrorCode ReconstructListAdd_2WaySplit(PetscFunctionList *flist, const char *name, ReconstructFunction_2WaySplit r)
    500: {
    501:   PetscFunctionBeginUser;
    502:   PetscCall(PetscFunctionListAdd(flist, name, r));
    503:   PetscFunctionReturn(PETSC_SUCCESS);
    504: }

    506: PetscErrorCode ReconstructListFind_2WaySplit(PetscFunctionList flist, const char *name, ReconstructFunction_2WaySplit *r)
    507: {
    508:   PetscFunctionBeginUser;
    509:   PetscCall(PetscFunctionListFind(flist, name, r));
    510:   PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
    511:   PetscFunctionReturn(PETSC_SUCCESS);
    512: }

    514: /* --------------------------------- Physics ------- */
    515: PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
    516: {
    517:   PetscFunctionBeginUser;
    518:   PetscCall(PetscFree(vctx));
    519:   PetscFunctionReturn(PETSC_SUCCESS);
    520: }

    522: /* --------------------------------- Finite Volume Solver --------------- */
    523: PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    524: {
    525:   FVCtx       *ctx = (FVCtx *)vctx;
    526:   PetscInt     i, j, k, Mx, dof, xs, xm;
    527:   PetscReal    hx, cfl_idt = 0;
    528:   PetscScalar *x, *f, *slope;
    529:   Vec          Xloc;
    530:   DM           da;

    532:   PetscFunctionBeginUser;
    533:   PetscCall(TSGetDM(ts, &da));
    534:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points */
    535:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */
    536:   hx = (ctx->xmax - ctx->xmin) / Mx;
    537:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
    538:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    539:   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
    540:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    541:   PetscCall(DMDAVecGetArray(da, F, &f));
    542:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
    543:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    545:   if (ctx->bctype == FVBC_OUTFLOW) {
    546:     for (i = xs - 2; i < 0; i++) {
    547:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    548:     }
    549:     for (i = Mx; i < xs + xm + 2; i++) {
    550:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    551:     }
    552:   }

    554:   for (i = xs - 1; i < xs + xm + 1; i++) {
    555:     struct _LimitInfo info;
    556:     PetscScalar      *cjmpL, *cjmpR;
    557:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    558:     PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
    559:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    560:     PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    561:     cjmpL = &ctx->cjmpLR[0];
    562:     cjmpR = &ctx->cjmpLR[dof];
    563:     for (j = 0; j < dof; j++) {
    564:       PetscScalar jmpL, jmpR;
    565:       jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    566:       jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    567:       for (k = 0; k < dof; k++) {
    568:         cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    569:         cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    570:       }
    571:     }
    572:     /* Apply limiter to the left and right characteristic jumps */
    573:     info.m  = dof;
    574:     info.hx = hx;
    575:     (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    576:     for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    577:     for (j = 0; j < dof; j++) {
    578:       PetscScalar tmp = 0;
    579:       for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    580:       slope[i * dof + j] = tmp;
    581:     }
    582:   }

    584:   for (i = xs; i < xs + xm + 1; i++) {
    585:     PetscReal    maxspeed;
    586:     PetscScalar *uL, *uR;
    587:     uL = &ctx->uLR[0];
    588:     uR = &ctx->uLR[dof];
    589:     for (j = 0; j < dof; j++) {
    590:       uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    591:       uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    592:     }
    593:     PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    594:     cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    595:     if (i > xs) {
    596:       for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
    597:     }
    598:     if (i < xs + xm) {
    599:       for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
    600:     }
    601:   }
    602:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    603:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    604:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    605:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    606:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    607:   if (0) {
    608:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
    609:     PetscReal dt, tnow;
    610:     PetscCall(TSGetTimeStep(ts, &dt));
    611:     PetscCall(TSGetTime(ts, &tnow));
    612:     if (dt > 0.5 / ctx->cfl_idt) {
    613:       if (1) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(1 / (2 * ctx->cfl_idt))));
    614:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
    615:     }
    616:   }
    617:   PetscFunctionReturn(PETSC_SUCCESS);
    618: }

    620: PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
    621: {
    622:   PetscScalar *u, *uj;
    623:   PetscInt     i, j, k, dof, xs, xm, Mx;

    625:   PetscFunctionBeginUser;
    626:   PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
    627:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    628:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    629:   PetscCall(DMDAVecGetArray(da, U, &u));
    630:   PetscCall(PetscMalloc1(dof, &uj));
    631:   for (i = xs; i < xs + xm; i++) {
    632:     const PetscReal h = (ctx->xmax - ctx->xmin) / Mx, xi = ctx->xmin + h / 2 + i * h;
    633:     const PetscInt  N = 200;
    634:     /* Integrate over cell i using trapezoid rule with N points. */
    635:     for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    636:     for (j = 0; j < N + 1; j++) {
    637:       PetscScalar xj = xi + h * (j - N / 2) / (PetscReal)N;
    638:       PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    639:       for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    640:     }
    641:   }
    642:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    643:   PetscCall(PetscFree(uj));
    644:   PetscFunctionReturn(PETSC_SUCCESS);
    645: }

    647: PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
    648: {
    649:   PetscReal          xmin, xmax;
    650:   PetscScalar        sum, tvsum, tvgsum;
    651:   const PetscScalar *x;
    652:   PetscInt           imin, imax, Mx, i, j, xs, xm, dof;
    653:   Vec                Xloc;
    654:   PetscBool          isascii;

    656:   PetscFunctionBeginUser;
    657:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    658:   PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
    659:   /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
    660:   PetscCall(DMGetLocalVector(da, &Xloc));
    661:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    662:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    663:   PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
    664:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    665:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    666:   tvsum = 0;
    667:   for (i = xs; i < xs + xm; i++) {
    668:     for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
    669:   }
    670:   PetscCallMPI(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)da)));
    671:   PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
    672:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    673:   PetscCall(VecMin(X, &imin, &xmin));
    674:   PetscCall(VecMax(X, &imax, &xmax));
    675:   PetscCall(VecSum(X, &sum));
    676:   PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%8.5f,%8.5f] with minimum at %" PetscInt_FMT ", mean %8.5f, ||x||_TV %8.5f\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvgsum / Mx)));
    677:   PetscFunctionReturn(PETSC_SUCCESS);
    678: }