Actual source code: ex4.c

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

      6: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
      7:                            "  advect      - Constant coefficient scalar advection\n"
      8:                            "                u_t       + (a*u)_x               = 0\n"
      9:                            "  shallow     - 1D Shallow water equations (Saint Venant System)\n"
     10:                            "                h_t + (q)_x = 0 \n"
     11:                            "                q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0  \n"
     12:                            "                where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n"
     13:                            "  for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
     14:                            "                hxs  = hratio*hxf \n"
     15:                            "  where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
     16:                            "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
     17:                            "                the states across shocks and rarefactions\n"
     18:                            "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
     19:                            "                also the reference solution should be generated by user and stored in a binary file.\n"
     20:                            "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
     21:                            "  bc_type     - Boundary condition for the problem, options are: periodic, outflow, inflow "
     22:                            "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n"
     23:                            "The problem size should be set with -da_grid_x M\n\n";

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

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

     36: #include <petscts.h>
     37: #include <petscdm.h>
     38: #include <petscdmda.h>
     39: #include <petscdraw.h>
     40: #include "finitevolume1d.h"
     41: #include <petsc/private/kernels/blockinvert.h>

     43: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
     44: {
     45:   PetscReal range = xmax - xmin;
     46:   return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
     47: }
     48: static inline PetscReal MaxAbs(PetscReal a, PetscReal b)
     49: {
     50:   return (PetscAbs(a) > PetscAbs(b)) ? a : b;
     51: }
     52: /* --------------------------------- Advection ----------------------------------- */
     53: typedef struct {
     54:   PetscReal a; /* advective velocity */
     55: } AdvectCtx;

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

     62:   PetscFunctionBeginUser;
     63:   speed     = ctx->a;
     64:   flux[0]   = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
     65:   *maxspeed = speed;
     66:   PetscFunctionReturn(PETSC_SUCCESS);
     67: }

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

     73:   PetscFunctionBeginUser;
     74:   X[0]      = 1.;
     75:   Xi[0]     = 1.;
     76:   speeds[0] = ctx->a;
     77:   PetscFunctionReturn(PETSC_SUCCESS);
     78: }

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

     85:   PetscFunctionBeginUser;
     86:   switch (bctype) {
     87:   case FVBC_OUTFLOW:
     88:     x0 = x - a * t;
     89:     break;
     90:   case FVBC_PERIODIC:
     91:     x0 = RangeMod(x - a * t, xmin, xmax);
     92:     break;
     93:   default:
     94:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
     95:   }
     96:   switch (initial) {
     97:   case 0:
     98:     u[0] = (x0 < 0) ? 1 : -1;
     99:     break;
    100:   case 1:
    101:     u[0] = (x0 < 0) ? -1 : 1;
    102:     break;
    103:   case 2:
    104:     u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
    105:     break;
    106:   case 3:
    107:     u[0] = PetscSinReal(2 * PETSC_PI * x0);
    108:     break;
    109:   case 4:
    110:     u[0] = PetscAbs(x0);
    111:     break;
    112:   case 5:
    113:     u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    114:     break;
    115:   case 6:
    116:     u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    117:     break;
    118:   case 7:
    119:     u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    120:     break;
    121:   default:
    122:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    123:   }
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
    128: {
    129:   AdvectCtx *user;

    131:   PetscFunctionBeginUser;
    132:   PetscCall(PetscNew(&user));
    133:   ctx->physics2.sample2         = PhysicsSample_Advect;
    134:   ctx->physics2.riemann2        = PhysicsRiemann_Advect;
    135:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
    136:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
    137:   ctx->physics2.user            = user;
    138:   ctx->physics2.dof             = 1;

    140:   PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0]));
    141:   user->a = 1;
    142:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    143:   PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
    144:   PetscOptionsEnd();
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }
    147: /* --------------------------------- Shallow Water ----------------------------------- */

    149: typedef struct {
    150:   PetscReal gravity;
    151: } ShallowCtx;

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

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

    168:   PetscFunctionBeginUser;
    169:   PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
    170:   cL = PetscSqrtScalar(g * L.h);
    171:   cR = PetscSqrtScalar(g * R.h);
    172:   c  = PetscMax(cL, cR);
    173:   {
    174:     /* Solve for star state */
    175:     const PetscInt maxits = 50;
    176:     PetscScalar    tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */
    177:     h0 = h;
    178:     for (i = 0; i < maxits; i++) {
    179:       PetscScalar fr, fl, dfr, dfl;
    180:       fl  = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */
    181:                       : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h);         /* rarefaction */
    182:       fr  = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */
    183:                       : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h);         /* rarefaction */
    184:       res = R.u - L.u + fr + fl;
    185:       PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
    186:       if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h - h0) < PETSC_SQRT_MACHINE_EPSILON)) {
    187:         star.h = h;
    188:         star.u = L.u - fl;
    189:         goto converged;
    190:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
    191:         h = 0.8 * h0 + 0.2 * h;
    192:         continue;
    193:       }
    194:       /* Accept the last step and take another */
    195:       res0 = res;
    196:       h0   = h;
    197:       dfl  = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h);
    198:       dfr  = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h);
    199:       tmp  = h - res / (dfr + dfl);
    200:       if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
    201:       else h = tmp;
    202:       PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h);
    203:     }
    204:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i);
    205:   }
    206: converged:
    207:   cstar = PetscSqrtScalar(g * star.h);
    208:   if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */
    209:     PetscScalar ufan[2];
    210:     ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL);
    211:     ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0];
    212:     ShallowFlux(phys, ufan, flux);
    213:   } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */
    214:     PetscScalar ufan[2];
    215:     ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR);
    216:     ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0];
    217:     ShallowFlux(phys, ufan, flux);
    218:   } 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)) {
    219:     /* 1-wave is right-travelling shock (supersonic) */
    220:     ShallowFlux(phys, uL, flux);
    221:   } 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)) {
    222:     /* 2-wave is left-travelling shock (supersonic) */
    223:     ShallowFlux(phys, uR, flux);
    224:   } else {
    225:     ustar[0] = star.h;
    226:     ustar[1] = star.h * star.u;
    227:     ShallowFlux(phys, ustar, flux);
    228:   }
    229:   *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR));
    230:   PetscFunctionReturn(PETSC_SUCCESS);
    231: }

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

    242:   PetscFunctionBeginUser;
    243:   /* Positivity preserving modification*/
    244:   if (L.h < tol) L.u = 0.0;
    245:   if (R.h < tol) R.u = 0.0;

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

    251:   ShallowFlux(phys, uL, fL);
    252:   ShallowFlux(phys, uR, fR);

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

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

    265:   PetscFunctionBeginUser;
    266:   for (i = 0; i < m; i++) {
    267:     for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j);
    268:     speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
    269:   }
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

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

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

    282:   if (u[0] < tol) { /*Use conservative variables*/
    283:     X[0 * 2 + 0] = 1;
    284:     X[0 * 2 + 1] = 0;
    285:     X[1 * 2 + 0] = 0;
    286:     X[1 * 2 + 1] = 1;
    287:   } else {
    288:     speeds[0]    = u[1] / u[0] - c;
    289:     speeds[1]    = u[1] / u[0] + c;
    290:     X[0 * 2 + 0] = 1;
    291:     X[0 * 2 + 1] = speeds[0];
    292:     X[1 * 2 + 0] = 1;
    293:     X[1 * 2 + 1] = speeds[1];
    294:   }

    296:   PetscCall(PetscArraycpy(Xi, X, 4));
    297:   PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
    298:   PetscFunctionReturn(PETSC_SUCCESS);
    299: }

    301: static PetscErrorCode PhysicsSample_Shallow(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    302: {
    303:   PetscFunctionBeginUser;
    304:   PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0");
    305:   switch (initial) {
    306:   case 0:
    307:     u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */
    308:     u[1] = (x < 0) ? 0 : 0;
    309:     break;
    310:   case 1:
    311:     u[0] = (x < 10) ? 1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */
    312:     u[1] = (x < 10) ? 2.5 : 0;
    313:     break;
    314:   case 2:
    315:     u[0] = (x < 25) ? 1 : 1;
    316:     u[1] = (x < 25) ? -5 : 5;
    317:     break;
    318:   case 3:
    319:     u[0] = (x < 20) ? 1 : 1e-12;
    320:     u[1] = (x < 20) ? 0 : 0;
    321:     break;
    322:   case 4:
    323:     u[0] = (x < 30) ? 1e-12 : 1;
    324:     u[1] = (x < 30) ? 0 : 0;
    325:     break;
    326:   case 5:
    327:     u[0] = (x < 25) ? 0.1 : 0.1;
    328:     u[1] = (x < 25) ? -0.3 : 0.3;
    329:     break;
    330:   case 6:
    331:     u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x);
    332:     u[1] = 1 * u[0];
    333:     break;
    334:   case 7:
    335:     if (x < -0.1) {
    336:       u[0] = 1e-9;
    337:       u[1] = 0.0;
    338:     } else if (x < 0.1) {
    339:       u[0] = 1.0;
    340:       u[1] = 0.0;
    341:     } else {
    342:       u[0] = 1e-9;
    343:       u[1] = 0.0;
    344:     }
    345:     break;
    346:   case 8:
    347:     if (x < -0.1) {
    348:       u[0] = 2;
    349:       u[1] = 0.0;
    350:     } else if (x < 0.1) {
    351:       u[0] = 3.0;
    352:       u[1] = 0.0;
    353:     } else {
    354:       u[0] = 2;
    355:       u[1] = 0.0;
    356:     }
    357:     break;
    358:   default:
    359:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    360:   }
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }

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

    370:   PetscFunctionBeginUser;
    371:   if (ctx->bctype == FVBC_INFLOW) {
    372:     switch (ctx->initial) {
    373:     case 0:
    374:     case 1:
    375:     case 2:
    376:     case 3:
    377:     case 4:
    378:     case 5:
    379:       u[0] = 0;
    380:       u[1] = 0.0; /* Left boundary conditions */
    381:       u[2] = 0;
    382:       u[3] = 0.0; /* Right boundary conditions */
    383:       break;
    384:     case 6:
    385:       u[0] = 0;
    386:       u[1] = 0.0; /* Left boundary conditions */
    387:       u[2] = 0;
    388:       u[3] = 0.0; /* Right boundary conditions */
    389:       break;
    390:     case 7:
    391:       u[0] = 0;
    392:       u[1] = 0.0; /* Left boundary conditions */
    393:       u[2] = 0;
    394:       u[3] = 0.0; /* Right boundary conditions */
    395:       break;
    396:     case 8:
    397:       u[0] = 0;
    398:       u[1] = 1.0; /* Left boundary conditions */
    399:       u[2] = 0;
    400:       u[3] = -1.0; /* Right boundary conditions */
    401:       break;
    402:     default:
    403:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    404:     }
    405:   }
    406:   PetscFunctionReturn(PETSC_SUCCESS);
    407: }

    409: /* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */
    410: static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx)
    411: {
    412:   PetscFunctionBeginUser;
    413:   switch (ctx->initial) {
    414:   case 0:
    415:   case 1:
    416:   case 2:
    417:   case 3:
    418:   case 4:
    419:   case 5:
    420:   case 6:
    421:   case 7:
    422:   case 8: /* Fix left and right momentum, height is outflow */
    423:     ctx->physics2.bcinflowindex[0] = PETSC_FALSE;
    424:     ctx->physics2.bcinflowindex[1] = PETSC_TRUE;
    425:     ctx->physics2.bcinflowindex[2] = PETSC_FALSE;
    426:     ctx->physics2.bcinflowindex[3] = PETSC_TRUE;
    427:     break;
    428:   default:
    429:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    430:   }
    431:   PetscFunctionReturn(PETSC_SUCCESS);
    432: }

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

    440:   PetscFunctionBeginUser;
    441:   PetscCall(PetscNew(&user));
    442:   ctx->physics2.sample2         = PhysicsSample_Shallow;
    443:   ctx->physics2.inflow          = PhysicsInflow_Shallow;
    444:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
    445:   ctx->physics2.riemann2        = PhysicsRiemann_Shallow_Rusanov;
    446:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow;
    447:   ctx->physics2.user            = user;
    448:   ctx->physics2.dof             = 2;

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

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

    456:   user->gravity = 9.81;

    458:   PetscCall(RiemannListAdd_2WaySplit(&rlist, "exact", PhysicsRiemann_Shallow_Exact));
    459:   PetscCall(RiemannListAdd_2WaySplit(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov));
    460:   PetscCall(ReconstructListAdd_2WaySplit(&rclist, "characteristic", PhysicsCharacteristic_Shallow));
    461:   PetscCall(ReconstructListAdd_2WaySplit(&rclist, "conservative", PhysicsCharacteristic_Conservative));
    462:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", "");
    463:   PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL));
    464:   PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
    465:   PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
    466:   PetscOptionsEnd();
    467:   PetscCall(RiemannListFind_2WaySplit(rlist, rname, &ctx->physics2.riemann2));
    468:   PetscCall(ReconstructListFind_2WaySplit(rclist, rcname, &ctx->physics2.characteristic2));
    469:   PetscCall(PetscFunctionListDestroy(&rlist));
    470:   PetscCall(PetscFunctionListDestroy(&rclist));
    471:   PetscFunctionReturn(PETSC_SUCCESS);
    472: }

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

    481:   PetscFunctionBeginUser;
    482:   PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
    483:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    484:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    485:   PetscCall(DMDAVecGetArray(da, U, &u));
    486:   PetscCall(PetscMalloc1(dof, &uj));
    487:   hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    488:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    489:   for (i = xs; i < xs + xm; i++) {
    490:     if (i < ctx->sf) {
    491:       xi = ctx->xmin + 0.5 * hs + i * hs;
    492:       /* Integrate over cell i using trapezoid rule with N points. */
    493:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    494:       for (j = 0; j < N + 1; j++) {
    495:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    496:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    497:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    498:       }
    499:     } else if (i < ctx->fs) {
    500:       xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf;
    501:       /* Integrate over cell i using trapezoid rule with N points. */
    502:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    503:       for (j = 0; j < N + 1; j++) {
    504:         xj = xi + hf * (j - N / 2) / (PetscReal)N;
    505:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    506:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    507:       }
    508:     } else {
    509:       xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs;
    510:       /* Integrate over cell i using trapezoid rule with N points. */
    511:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    512:       for (j = 0; j < N + 1; j++) {
    513:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    514:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    515:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    516:       }
    517:     }
    518:   }
    519:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    520:   PetscCall(PetscFree(uj));
    521:   PetscFunctionReturn(PETSC_SUCCESS);
    522: }

    524: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
    525: {
    526:   Vec                Y;
    527:   PetscInt           i, Mx;
    528:   const PetscScalar *ptr_X, *ptr_Y;
    529:   PetscReal          hs, hf;

    531:   PetscFunctionBeginUser;
    532:   PetscCall(VecGetSize(X, &Mx));
    533:   PetscCall(VecDuplicate(X, &Y));
    534:   PetscCall(FVSample_2WaySplit(ctx, da, t, Y));
    535:   hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    536:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    537:   PetscCall(VecGetArrayRead(X, &ptr_X));
    538:   PetscCall(VecGetArrayRead(Y, &ptr_Y));
    539:   for (i = 0; i < Mx; i++) {
    540:     if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
    541:     else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
    542:   }
    543:   PetscCall(VecRestoreArrayRead(X, &ptr_X));
    544:   PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
    545:   PetscCall(VecDestroy(&Y));
    546:   PetscFunctionReturn(PETSC_SUCCESS);
    547: }

    549: PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    550: {
    551:   FVCtx       *ctx = (FVCtx *)vctx;
    552:   PetscInt     i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
    553:   PetscReal    hxf, hxs, cfl_idt = 0;
    554:   PetscScalar *x, *f, *slope;
    555:   Vec          Xloc;
    556:   DM           da;

    558:   PetscFunctionBeginUser;
    559:   PetscCall(TSGetDM(ts, &da));
    560:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
    561:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
    562:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    563:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    564:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
    565:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

    567:   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */

    569:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    570:   PetscCall(DMDAVecGetArray(da, F, &f));
    571:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points                                           */
    572:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    574:   if (ctx->bctype == FVBC_OUTFLOW) {
    575:     for (i = xs - 2; i < 0; i++) {
    576:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    577:     }
    578:     for (i = Mx; i < xs + xm + 2; i++) {
    579:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    580:     }
    581:   }

    583:   if (ctx->bctype == FVBC_INFLOW) {
    584:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    585:     pages 137-138 for the scheme. */
    586:     if (xs == 0) { /* Left Boundary */
    587:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
    588:       for (j = 0; j < dof; j++) {
    589:         if (ctx->physics2.bcinflowindex[j]) {
    590:           for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
    591:         } else {
    592:           for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
    593:         }
    594:       }
    595:     }
    596:     if (xs + xm == Mx) { /* Right Boundary */
    597:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
    598:       for (j = 0; j < dof; j++) {
    599:         if (ctx->physics2.bcinflowindex[dof + j]) {
    600:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
    601:         } else {
    602:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
    603:         }
    604:       }
    605:     }
    606:   }

    608:   for (i = xs - 1; i < xs + xm + 1; i++) {
    609:     struct _LimitInfo info;
    610:     PetscScalar      *cjmpL, *cjmpR;
    611:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    612:     PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    613:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    614:     PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    615:     cjmpL = &ctx->cjmpLR[0];
    616:     cjmpR = &ctx->cjmpLR[dof];
    617:     for (j = 0; j < dof; j++) {
    618:       PetscScalar jmpL, jmpR;
    619:       jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    620:       jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    621:       for (k = 0; k < dof; k++) {
    622:         cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    623:         cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    624:       }
    625:     }
    626:     /* Apply limiter to the left and right characteristic jumps */
    627:     info.m   = dof;
    628:     info.hxs = hxs;
    629:     info.hxf = hxf;
    630:     (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    631:     for (j = 0; j < dof; j++) {
    632:       PetscScalar tmp = 0;
    633:       for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    634:       slope[i * dof + j] = tmp;
    635:     }
    636:   }

    638:   for (i = xs; i < xs + xm + 1; i++) {
    639:     PetscReal    maxspeed;
    640:     PetscScalar *uL, *uR;
    641:     uL = &ctx->uLR[0];
    642:     uR = &ctx->uLR[dof];
    643:     if (i < sf) { /* slow region */
    644:       for (j = 0; j < dof; j++) {
    645:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    646:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    647:       }
    648:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    649:       if (i > xs) {
    650:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    651:       }
    652:       if (i < xs + xm) {
    653:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    654:       }
    655:     } else if (i == sf) { /* interface between the slow region and the fast region */
    656:       for (j = 0; j < dof; j++) {
    657:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    658:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    659:       }
    660:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    661:       if (i > xs) {
    662:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    663:       }
    664:       if (i < xs + xm) {
    665:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
    666:       }
    667:     } else if (i > sf && i < fs) { /* fast region */
    668:       for (j = 0; j < dof; j++) {
    669:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    670:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    671:       }
    672:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    673:       if (i > xs) {
    674:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
    675:       }
    676:       if (i < xs + xm) {
    677:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
    678:       }
    679:     } else if (i == fs) { /* interface between the fast region and the slow region */
    680:       for (j = 0; j < dof; j++) {
    681:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    682:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    683:       }
    684:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    685:       if (i > xs) {
    686:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
    687:       }
    688:       if (i < xs + xm) {
    689:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    690:       }
    691:     } else { /* slow region */
    692:       for (j = 0; j < dof; j++) {
    693:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    694:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    695:       }
    696:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    697:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
    698:       if (i > xs) {
    699:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    700:       }
    701:       if (i < xs + xm) {
    702:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    703:       }
    704:     }
    705:   }
    706:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    707:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    708:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    709:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    710:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    711:   if (0) {
    712:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
    713:     PetscReal dt, tnow;
    714:     PetscCall(TSGetTimeStep(ts, &dt));
    715:     PetscCall(TSGetTime(ts, &tnow));
    716:     if (dt > 0.5 / ctx->cfl_idt) {
    717:       if (1) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
    718:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
    719:     }
    720:   }
    721:   PetscFunctionReturn(PETSC_SUCCESS);
    722: }

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

    734:   PetscFunctionBeginUser;
    735:   PetscCall(TSGetDM(ts, &da));
    736:   PetscCall(DMGetLocalVector(da, &Xloc));
    737:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    738:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    739:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    740:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    741:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    742:   PetscCall(VecZeroEntries(F));
    743:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    744:   PetscCall(VecGetArray(F, &f));
    745:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    746:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    748:   if (ctx->bctype == FVBC_OUTFLOW) {
    749:     for (i = xs - 2; i < 0; i++) {
    750:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    751:     }
    752:     for (i = Mx; i < xs + xm + 2; i++) {
    753:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    754:     }
    755:   }
    756:   if (ctx->bctype == FVBC_INFLOW) {
    757:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    758:     pages 137-138 for the scheme. */
    759:     if (xs == 0) { /* Left Boundary */
    760:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
    761:       for (j = 0; j < dof; j++) {
    762:         if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
    763:           for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
    764:         } else {
    765:           for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
    766:         }
    767:       }
    768:     }
    769:     if (xs + xm == Mx) { /* Right Boundary */
    770:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
    771:       for (j = 0; j < dof; j++) {
    772:         if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
    773:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
    774:         } else {
    775:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
    776:         }
    777:       }
    778:     }
    779:   }
    780:   for (i = xs - 1; i < xs + xm + 1; i++) {
    781:     struct _LimitInfo info;
    782:     PetscScalar      *cjmpL, *cjmpR;
    783:     if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
    784:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    785:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    786:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    787:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    788:       cjmpL = &ctx->cjmpLR[0];
    789:       cjmpR = &ctx->cjmpLR[dof];
    790:       for (j = 0; j < dof; j++) {
    791:         PetscScalar jmpL, jmpR;
    792:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    793:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    794:         for (k = 0; k < dof; k++) {
    795:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    796:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    797:         }
    798:       }
    799:       /* Apply limiter to the left and right characteristic jumps */
    800:       info.m   = dof;
    801:       info.hxs = hxs;
    802:       info.hxf = hxf;
    803:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    804:       for (j = 0; j < dof; j++) {
    805:         PetscScalar tmp = 0;
    806:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    807:         slope[i * dof + j] = tmp;
    808:       }
    809:     }
    810:   }

    812:   for (i = xs; i < xs + xm + 1; i++) {
    813:     PetscReal    maxspeed;
    814:     PetscScalar *uL, *uR;
    815:     uL = &ctx->uLR[0];
    816:     uR = &ctx->uLR[dof];
    817:     if (i < sf - lsbwidth) { /* slow region */
    818:       for (j = 0; j < dof; j++) {
    819:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    820:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    821:       }
    822:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    823:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
    824:       if (i > xs) {
    825:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    826:       }
    827:       if (i < xs + xm) {
    828:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    829:         islow++;
    830:       }
    831:     }
    832:     if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
    833:       for (j = 0; j < dof; j++) {
    834:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    835:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    836:       }
    837:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    838:       if (i > xs) {
    839:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    840:       }
    841:     }
    842:     if (i == fs + rsbwidth) { /* slow region */
    843:       for (j = 0; j < dof; j++) {
    844:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    845:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    846:       }
    847:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    848:       if (i < xs + xm) {
    849:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    850:         islow++;
    851:       }
    852:     }
    853:     if (i > fs + rsbwidth) { /* slow region */
    854:       for (j = 0; j < dof; j++) {
    855:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    856:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    857:       }
    858:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    859:       if (i > xs) {
    860:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    861:       }
    862:       if (i < xs + xm) {
    863:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    864:         islow++;
    865:       }
    866:     }
    867:   }
    868:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    869:   PetscCall(VecRestoreArray(F, &f));
    870:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    871:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    872:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    873:   PetscFunctionReturn(PETSC_SUCCESS);
    874: }

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

    885:   PetscFunctionBeginUser;
    886:   PetscCall(TSGetDM(ts, &da));
    887:   PetscCall(DMGetLocalVector(da, &Xloc));
    888:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    889:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    890:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    891:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    892:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    893:   PetscCall(VecZeroEntries(F));
    894:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    895:   PetscCall(VecGetArray(F, &f));
    896:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    897:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    899:   if (ctx->bctype == FVBC_OUTFLOW) {
    900:     for (i = xs - 2; i < 0; i++) {
    901:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    902:     }
    903:     for (i = Mx; i < xs + xm + 2; i++) {
    904:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    905:     }
    906:   }
    907:   if (ctx->bctype == FVBC_INFLOW) {
    908:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    909:     pages 137-138 for the scheme. */
    910:     if (xs == 0) { /* Left Boundary */
    911:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
    912:       for (j = 0; j < dof; j++) {
    913:         if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
    914:           for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
    915:         } else {
    916:           for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
    917:         }
    918:       }
    919:     }
    920:     if (xs + xm == Mx) { /* Right Boundary */
    921:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
    922:       for (j = 0; j < dof; j++) {
    923:         if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
    924:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
    925:         } else {
    926:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
    927:         }
    928:       }
    929:     }
    930:   }
    931:   for (i = xs - 1; i < xs + xm + 1; i++) {
    932:     struct _LimitInfo info;
    933:     PetscScalar      *cjmpL, *cjmpR;
    934:     if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
    935:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    936:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    937:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    938:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    939:       cjmpL = &ctx->cjmpLR[0];
    940:       cjmpR = &ctx->cjmpLR[dof];
    941:       for (j = 0; j < dof; j++) {
    942:         PetscScalar jmpL, jmpR;
    943:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    944:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    945:         for (k = 0; k < dof; k++) {
    946:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    947:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    948:         }
    949:       }
    950:       /* Apply limiter to the left and right characteristic jumps */
    951:       info.m   = dof;
    952:       info.hxs = hxs;
    953:       info.hxf = hxf;
    954:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    955:       for (j = 0; j < dof; j++) {
    956:         PetscScalar tmp = 0;
    957:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    958:         slope[i * dof + j] = tmp;
    959:       }
    960:     }
    961:   }

    963:   for (i = xs; i < xs + xm + 1; i++) {
    964:     PetscReal    maxspeed;
    965:     PetscScalar *uL, *uR;
    966:     uL = &ctx->uLR[0];
    967:     uR = &ctx->uLR[dof];
    968:     if (i == sf - lsbwidth) {
    969:       for (j = 0; j < dof; j++) {
    970:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    971:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    972:       }
    973:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    974:       if (i < xs + xm) {
    975:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    976:         islow++;
    977:       }
    978:     }
    979:     if (i > sf - lsbwidth && i < sf) {
    980:       for (j = 0; j < dof; j++) {
    981:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    982:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    983:       }
    984:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    985:       if (i > xs) {
    986:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    987:       }
    988:       if (i < xs + xm) {
    989:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    990:         islow++;
    991:       }
    992:     }
    993:     if (i == sf) { /* interface between the slow region and the fast region */
    994:       for (j = 0; j < dof; j++) {
    995:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    996:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    997:       }
    998:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    999:       if (i > xs) {
    1000:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    1001:       }
    1002:     }
    1003:     if (i == fs) { /* interface between the fast region and the slow region */
    1004:       for (j = 0; j < dof; j++) {
    1005:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    1006:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    1007:       }
    1008:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1009:       if (i < xs + xm) {
    1010:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    1011:         islow++;
    1012:       }
    1013:     }
    1014:     if (i > fs && i < fs + rsbwidth) {
    1015:       for (j = 0; j < dof; j++) {
    1016:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    1017:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    1018:       }
    1019:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1020:       if (i > xs) {
    1021:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    1022:       }
    1023:       if (i < xs + xm) {
    1024:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    1025:         islow++;
    1026:       }
    1027:     }
    1028:     if (i == fs + rsbwidth) {
    1029:       for (j = 0; j < dof; j++) {
    1030:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    1031:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    1032:       }
    1033:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1034:       if (i > xs) {
    1035:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    1036:       }
    1037:     }
    1038:   }
    1039:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    1040:   PetscCall(VecRestoreArray(F, &f));
    1041:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    1042:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    1043:   PetscFunctionReturn(PETSC_SUCCESS);
    1044: }

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

    1056:   PetscFunctionBeginUser;
    1057:   PetscCall(TSGetDM(ts, &da));
    1058:   PetscCall(DMGetLocalVector(da, &Xloc));
    1059:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    1060:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    1061:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    1062:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    1063:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    1064:   PetscCall(VecZeroEntries(F));
    1065:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    1066:   PetscCall(VecGetArray(F, &f));
    1067:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    1068:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    1070:   if (ctx->bctype == FVBC_OUTFLOW) {
    1071:     for (i = xs - 2; i < 0; i++) {
    1072:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    1073:     }
    1074:     for (i = Mx; i < xs + xm + 2; i++) {
    1075:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    1076:     }
    1077:   }
    1078:   if (ctx->bctype == FVBC_INFLOW) {
    1079:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
    1080:     pages 137-138 for the scheme. */
    1081:     if (xs == 0) { /* Left Boundary */
    1082:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
    1083:       for (j = 0; j < dof; j++) {
    1084:         if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
    1085:           for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
    1086:         } else {
    1087:           for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
    1088:         }
    1089:       }
    1090:     }
    1091:     if (xs + xm == Mx) { /* Right Boundary */
    1092:       PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
    1093:       for (j = 0; j < dof; j++) {
    1094:         if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
    1095:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
    1096:         } else {
    1097:           for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
    1098:         }
    1099:       }
    1100:     }
    1101:   }
    1102:   for (i = xs - 1; i < xs + xm + 1; i++) {
    1103:     struct _LimitInfo info;
    1104:     PetscScalar      *cjmpL, *cjmpR;
    1105:     if (i > sf - 2 && i < fs + 1) {
    1106:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    1107:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    1108:       cjmpL = &ctx->cjmpLR[0];
    1109:       cjmpR = &ctx->cjmpLR[dof];
    1110:       for (j = 0; j < dof; j++) {
    1111:         PetscScalar jmpL, jmpR;
    1112:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    1113:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    1114:         for (k = 0; k < dof; k++) {
    1115:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    1116:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    1117:         }
    1118:       }
    1119:       /* Apply limiter to the left and right characteristic jumps */
    1120:       info.m   = dof;
    1121:       info.hxs = hxs;
    1122:       info.hxf = hxf;
    1123:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    1124:       for (j = 0; j < dof; j++) {
    1125:         PetscScalar tmp = 0;
    1126:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    1127:         slope[i * dof + j] = tmp;
    1128:       }
    1129:     }
    1130:   }

    1132:   for (i = xs; i < xs + xm + 1; i++) {
    1133:     PetscReal    maxspeed;
    1134:     PetscScalar *uL, *uR;
    1135:     uL = &ctx->uLR[0];
    1136:     uR = &ctx->uLR[dof];
    1137:     if (i == sf) { /* interface between the slow region and the fast region */
    1138:       for (j = 0; j < dof; j++) {
    1139:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    1140:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    1141:       }
    1142:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1143:       if (i < xs + xm) {
    1144:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
    1145:         ifast++;
    1146:       }
    1147:     }
    1148:     if (i > sf && i < fs) { /* fast region */
    1149:       for (j = 0; j < dof; j++) {
    1150:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    1151:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    1152:       }
    1153:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1154:       if (i > xs) {
    1155:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
    1156:       }
    1157:       if (i < xs + xm) {
    1158:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
    1159:         ifast++;
    1160:       }
    1161:     }
    1162:     if (i == fs) { /* interface between the fast region and the slow region */
    1163:       for (j = 0; j < dof; j++) {
    1164:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    1165:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    1166:       }
    1167:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    1168:       if (i > xs) {
    1169:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
    1170:       }
    1171:     }
    1172:   }
    1173:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    1174:   PetscCall(VecRestoreArray(F, &f));
    1175:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    1176:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    1177:   PetscFunctionReturn(PETSC_SUCCESS);
    1178: }

    1180: int main(int argc, char *argv[])
    1181: {
    1182:   char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
    1183:   PetscFunctionList limiters = 0, physics = 0;
    1184:   MPI_Comm          comm;
    1185:   TS                ts;
    1186:   DM                da;
    1187:   Vec               X, X0, R;
    1188:   FVCtx             ctx;
    1189:   PetscInt          i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, islowbuffer = 0, *index_slow, *index_fast, *index_slowbuffer;
    1190:   PetscBool         view_final = PETSC_FALSE;
    1191:   PetscReal         ptime, maxtime;

    1193:   PetscFunctionBeginUser;
    1194:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
    1195:   comm = PETSC_COMM_WORLD;
    1196:   PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

    1198:   /* Register limiters to be available on the command line */
    1199:   PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
    1200:   PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
    1201:   PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
    1202:   PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
    1203:   PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
    1204:   PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
    1205:   PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
    1206:   PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));

    1208:   /* Register physical models to be available on the command line */
    1209:   PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow));
    1210:   PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));

    1212:   ctx.comm       = comm;
    1213:   ctx.cfl        = 0.9;
    1214:   ctx.bctype     = FVBC_PERIODIC;
    1215:   ctx.xmin       = -1.0;
    1216:   ctx.xmax       = 1.0;
    1217:   ctx.initial    = 1;
    1218:   ctx.hratio     = 2;
    1219:   maxtime        = 10.0;
    1220:   ctx.simulation = PETSC_FALSE;
    1221:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
    1222:   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
    1223:   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
    1224:   PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
    1225:   PetscCall(PetscOptionsFList("-physics", "Name of physics model to use", "", physics, physname, physname, sizeof(physname), NULL));
    1226:   PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
    1227:   PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
    1228:   PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
    1229:   PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
    1230:   PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
    1231:   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
    1232:   PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
    1233:   PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
    1234:   PetscOptionsEnd();

    1236:   /* Choose the limiter from the list of registered limiters */
    1237:   PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2));
    1238:   PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);

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

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

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

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

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

    1272:   /* create index for slow parts and fast parts,
    1273:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
    1274:   count_slow = Mx * 3 / (3 + ctx.hratio); // compute Mx / (1.0 + ctx.hratio / 3.0);
    1275:   PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hratio/3) is even");
    1276:   count_fast = Mx - count_slow;
    1277:   ctx.sf     = count_slow / 2;
    1278:   ctx.fs     = ctx.sf + count_fast;
    1279:   PetscCall(PetscMalloc1(xm * dof, &index_slow));
    1280:   PetscCall(PetscMalloc1(xm * dof, &index_fast));
    1281:   PetscCall(PetscMalloc1(8 * dof, &index_slowbuffer));
    1282:   ctx.lsbwidth = 4;
    1283:   ctx.rsbwidth = 4;

    1285:   for (i = xs; i < xs + xm; i++) {
    1286:     if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
    1287:       for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    1288:     else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
    1289:       for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
    1290:     else
    1291:       for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
    1292:   }
    1293:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
    1294:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
    1295:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));

    1297:   /* Create a time-stepping object */
    1298:   PetscCall(TSCreate(comm, &ts));
    1299:   PetscCall(TSSetDM(ts, da));
    1300:   PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
    1301:   PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
    1302:   PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
    1303:   PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
    1304:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
    1305:   PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
    1306:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));

    1308:   PetscCall(TSSetType(ts, TSMPRK));
    1309:   PetscCall(TSSetMaxTime(ts, maxtime));
    1310:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    1312:   /* Compute initial conditions and starting time step */
    1313:   PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
    1314:   PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
    1315:   PetscCall(VecCopy(X0, X));                                      /* The function value was not used so we set X=X0 again */
    1316:   PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
    1317:   PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
    1318:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    1319:   {
    1320:     PetscInt           steps;
    1321:     PetscScalar        mass_initial, mass_final, mass_difference, mass_differenceg;
    1322:     const PetscScalar *ptr_X, *ptr_X0;
    1323:     const PetscReal    hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow;
    1324:     const PetscReal    hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;

    1326:     PetscCall(TSSolve(ts, X));
    1327:     PetscCall(TSGetSolveTime(ts, &ptime));
    1328:     PetscCall(TSGetStepNumber(ts, &steps));
    1329:     /* calculate the total mass at initial time and final time */
    1330:     mass_initial = 0.0;
    1331:     mass_final   = 0.0;
    1332:     PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
    1333:     PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
    1334:     for (i = xs; i < xs + xm; i++) {
    1335:       if (i < ctx.sf || i > ctx.fs - 1) {
    1336:         for (k = 0; k < dof; k++) {
    1337:           mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
    1338:           mass_final   = mass_final + hs * ptr_X[i * dof + k];
    1339:         }
    1340:       } else {
    1341:         for (k = 0; k < dof; k++) {
    1342:           mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
    1343:           mass_final   = mass_final + hf * ptr_X[i * dof + k];
    1344:         }
    1345:       }
    1346:     }
    1347:     PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
    1348:     PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
    1349:     mass_difference = mass_final - mass_initial;
    1350:     PetscCallMPI(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
    1351:     PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
    1352:     PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    1353:     PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1.0 / ctx.cfl_idt)));
    1354:     if (ctx.exact) {
    1355:       PetscReal nrm1 = 0;
    1356:       PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1));
    1357:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    1358:     }
    1359:     if (ctx.simulation) {
    1360:       PetscReal          nrm1 = 0;
    1361:       PetscViewer        fd;
    1362:       char               filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
    1363:       Vec                XR;
    1364:       PetscBool          flg;
    1365:       const PetscScalar *ptr_XR;
    1366:       PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
    1367:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
    1368:       PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
    1369:       PetscCall(VecDuplicate(X0, &XR));
    1370:       PetscCall(VecLoad(XR, fd));
    1371:       PetscCall(PetscViewerDestroy(&fd));
    1372:       PetscCall(VecGetArrayRead(X, &ptr_X));
    1373:       PetscCall(VecGetArrayRead(XR, &ptr_XR));
    1374:       for (i = xs; i < xs + xm; i++) {
    1375:         if (i < ctx.sf || i > ctx.fs - 1)
    1376:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
    1377:         else
    1378:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
    1379:       }
    1380:       PetscCall(VecRestoreArrayRead(X, &ptr_X));
    1381:       PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
    1382:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    1383:       PetscCall(VecDestroy(&XR));
    1384:     }
    1385:   }

    1387:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    1388:   if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
    1389:   if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
    1390:   if (draw & 0x4) {
    1391:     Vec Y;
    1392:     PetscCall(VecDuplicate(X, &Y));
    1393:     PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y));
    1394:     PetscCall(VecAYPX(Y, -1, X));
    1395:     PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    1396:     PetscCall(VecDestroy(&Y));
    1397:   }

    1399:   if (view_final) {
    1400:     PetscViewer viewer;
    1401:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    1402:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    1403:     PetscCall(VecView(X, viewer));
    1404:     PetscCall(PetscViewerPopFormat(viewer));
    1405:     PetscCall(PetscViewerDestroy(&viewer));
    1406:   }

    1408:   /* Clean up */
    1409:   PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
    1410:   for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
    1411:   PetscCall(PetscFree(ctx.physics2.bcinflowindex));
    1412:   PetscCall(PetscFree(ctx.ub));
    1413:   PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
    1414:   PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
    1415:   PetscCall(VecDestroy(&X));
    1416:   PetscCall(VecDestroy(&X0));
    1417:   PetscCall(VecDestroy(&R));
    1418:   PetscCall(DMDestroy(&da));
    1419:   PetscCall(TSDestroy(&ts));
    1420:   PetscCall(ISDestroy(&ctx.iss));
    1421:   PetscCall(ISDestroy(&ctx.isf));
    1422:   PetscCall(ISDestroy(&ctx.issb));
    1423:   PetscCall(PetscFree(index_slow));
    1424:   PetscCall(PetscFree(index_fast));
    1425:   PetscCall(PetscFree(index_slowbuffer));
    1426:   PetscCall(PetscFunctionListDestroy(&limiters));
    1427:   PetscCall(PetscFunctionListDestroy(&physics));
    1428:   PetscCall(PetscFinalize());
    1429:   return 0;
    1430: }

    1432: /*TEST

    1434:     build:
    1435:       requires: !complex !single
    1436:       depends: finitevolume1d.c

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

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

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

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

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