Actual source code: ex5.c

      1: /*
      2:   Note:
      3:   To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_dt 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin
      4:   Errors can be computed in the following runs with -simulation -f reference.bin

      6:   Multirate RK options:
      7:   -ts_rk_dtratio is the ratio between larger time step size and small time step size
      8:   -ts_rk_multirate_type has three choices: none (for single step RK)
      9:                                            combined (for multirate method and user just need to provide the same RHS with the single step RK)
     10:                                            partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components
     11: */

     13: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
     14:                            " advection   - Variable coefficient scalar advection\n"
     15:                            "                u_t       + (a*u)_x               = 0\n"
     16:                            " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n"
     17:                            " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n"
     18:                            " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n"
     19:                            " more precisely, you need firstly generate a reference to a binary file say file.bin, then on command line,\n"
     20:                            " you should type -simulation -f file.bin.\n"
     21:                            " you can choose the number of grids by -da_grid_x.\n"
     22:                            " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n";

     24: #include <petscts.h>
     25: #include <petscdm.h>
     26: #include <petscdmda.h>
     27: #include <petscdraw.h>
     28: #include <petsc/private/tsimpl.h>

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

     32: #include "finitevolume1d.h"

     34: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
     35: {
     36:   PetscReal range = xmax - xmin;
     37:   return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
     38: }

     40: /* --------------------------------- Advection ----------------------------------- */

     42: typedef struct {
     43:   PetscReal a[2]; /* advective velocity */
     44: } AdvectCtx;

     46: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed, PetscReal x, PetscReal xmin, PetscReal xmax)
     47: {
     48:   AdvectCtx *ctx = (AdvectCtx *)vctx;
     49:   PetscReal *speed;

     51:   PetscFunctionBeginUser;
     52:   speed = ctx->a;
     53:   if (x == 0 || x == xmin || x == xmax)
     54:     flux[0] = PetscMax(0, (speed[0] + speed[1]) / 2.0) * uL[0] + PetscMin(0, (speed[0] + speed[1]) / 2.0) * uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */
     55:   else if (x < 0) flux[0] = PetscMax(0, speed[0]) * uL[0] + PetscMin(0, speed[0]) * uR[0];                         /* else if condition need to be changed based on different problem, 'x = 0' is discontinuous point of a */
     56:   else flux[0] = PetscMax(0, speed[1]) * uL[0] + PetscMin(0, speed[1]) * uR[0];
     57:   *maxspeed = *speed;
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

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

     65:   PetscFunctionBeginUser;
     66:   X[0]  = 1.;
     67:   Xi[0] = 1.;
     68:   if (x < 0) speeds[0] = ctx->a[0]; /* x is discontinuous point of a */
     69:   else speeds[0] = ctx->a[1];
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
     74: {
     75:   AdvectCtx *ctx = (AdvectCtx *)vctx;
     76:   PetscReal *a   = ctx->a, x0;

     78:   PetscFunctionBeginUser;
     79:   if (x < 0) { /* x is cell center */
     80:     switch (bctype) {
     81:     case FVBC_OUTFLOW:
     82:       x0 = x - a[0] * t;
     83:       break;
     84:     case FVBC_PERIODIC:
     85:       x0 = RangeMod(x - a[0] * t, xmin, xmax);
     86:       break;
     87:     default:
     88:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
     89:     }
     90:     switch (initial) {
     91:     case 0:
     92:       u[0] = (x0 < 0) ? 1 : -1;
     93:       break;
     94:     case 1:
     95:       u[0] = (x0 < 0) ? -1 : 1;
     96:       break;
     97:     case 2:
     98:       u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
     99:       break;
    100:     case 3:
    101:       u[0] = PetscSinReal(2 * PETSC_PI * x0);
    102:       break;
    103:     case 4:
    104:       u[0] = PetscAbs(x0);
    105:       break;
    106:     case 5:
    107:       u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    108:       break;
    109:     case 6:
    110:       u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    111:       break;
    112:     case 7:
    113:       u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    114:       break;
    115:     default:
    116:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    117:     }
    118:   } else {
    119:     switch (bctype) {
    120:     case FVBC_OUTFLOW:
    121:       x0 = x - a[1] * t;
    122:       break;
    123:     case FVBC_PERIODIC:
    124:       x0 = RangeMod(x - a[1] * t, xmin, xmax);
    125:       break;
    126:     default:
    127:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
    128:     }
    129:     switch (initial) {
    130:     case 0:
    131:       u[0] = (x0 < 0) ? 1 : -1;
    132:       break;
    133:     case 1:
    134:       u[0] = (x0 < 0) ? -1 : 1;
    135:       break;
    136:     case 2:
    137:       u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
    138:       break;
    139:     case 3:
    140:       u[0] = PetscSinReal(2 * PETSC_PI * x0);
    141:       break;
    142:     case 4:
    143:       u[0] = PetscAbs(x0);
    144:       break;
    145:     case 5:
    146:       u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    147:       break;
    148:     case 6:
    149:       u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    150:       break;
    151:     case 7:
    152:       u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    153:       break;
    154:     default:
    155:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    156:     }
    157:   }
    158:   PetscFunctionReturn(PETSC_SUCCESS);
    159: }

    161: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
    162: {
    163:   AdvectCtx *user;

    165:   PetscFunctionBeginUser;
    166:   PetscCall(PetscNew(&user));
    167:   ctx->physics.sample         = PhysicsSample_Advect;
    168:   ctx->physics.riemann        = PhysicsRiemann_Advect;
    169:   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
    170:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
    171:   ctx->physics.user           = user;
    172:   ctx->physics.dof            = 1;
    173:   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
    174:   user->a[0] = 1;
    175:   user->a[1] = 1;
    176:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    177:   {
    178:     PetscCall(PetscOptionsReal("-physics_advect_a1", "Speed1", "", user->a[0], &user->a[0], NULL));
    179:     PetscCall(PetscOptionsReal("-physics_advect_a2", "Speed2", "", user->a[1], &user->a[1], NULL));
    180:   }
    181:   PetscOptionsEnd();
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */

    187: PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    188: {
    189:   FVCtx       *ctx = (FVCtx *)vctx;
    190:   PetscInt     i, j, k, Mx, dof, xs, xm, len_slow;
    191:   PetscReal    hx, cfl_idt = 0;
    192:   PetscScalar *x, *f, *slope;
    193:   Vec          Xloc;
    194:   DM           da;

    196:   PetscFunctionBeginUser;
    197:   PetscCall(TSGetDM(ts, &da));
    198:   PetscCall(DMGetLocalVector(da, &Xloc));
    199:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    200:   hx = (ctx->xmax - ctx->xmin) / Mx;
    201:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    202:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    203:   PetscCall(VecZeroEntries(F));
    204:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    205:   PetscCall(VecGetArray(F, &f));
    206:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    207:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    208:   PetscCall(ISGetSize(ctx->iss, &len_slow));

    210:   if (ctx->bctype == FVBC_OUTFLOW) {
    211:     for (i = xs - 2; i < 0; i++) {
    212:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    213:     }
    214:     for (i = Mx; i < xs + xm + 2; i++) {
    215:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    216:     }
    217:   }
    218:   for (i = xs - 1; i < xs + xm + 1; i++) {
    219:     struct _LimitInfo info;
    220:     PetscScalar      *cjmpL, *cjmpR;
    221:     if (i < len_slow + 1) {
    222:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    223:       PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
    224:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    225:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    226:       cjmpL = &ctx->cjmpLR[0];
    227:       cjmpR = &ctx->cjmpLR[dof];
    228:       for (j = 0; j < dof; j++) {
    229:         PetscScalar jmpL, jmpR;
    230:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    231:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    232:         for (k = 0; k < dof; k++) {
    233:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    234:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    235:         }
    236:       }
    237:       /* Apply limiter to the left and right characteristic jumps */
    238:       info.m  = dof;
    239:       info.hx = hx;
    240:       (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    241:       for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    242:       for (j = 0; j < dof; j++) {
    243:         PetscScalar tmp = 0;
    244:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    245:         slope[i * dof + j] = tmp;
    246:       }
    247:     }
    248:   }

    250:   for (i = xs; i < xs + xm + 1; i++) {
    251:     PetscReal    maxspeed;
    252:     PetscScalar *uL, *uR;
    253:     if (i < len_slow) { /* slow parts can be changed based on a */
    254:       uL = &ctx->uLR[0];
    255:       uR = &ctx->uLR[dof];
    256:       for (j = 0; j < dof; j++) {
    257:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    258:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    259:       }
    260:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    261:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    262:       if (i > xs) {
    263:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
    264:       }
    265:       if (i < xs + xm) {
    266:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
    267:       }
    268:     }
    269:     if (i == len_slow) { /* slow parts can be changed based on a */
    270:       uL = &ctx->uLR[0];
    271:       uR = &ctx->uLR[dof];
    272:       for (j = 0; j < dof; j++) {
    273:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    274:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    275:       }
    276:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    277:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    278:       if (i > xs) {
    279:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
    280:       }
    281:     }
    282:   }
    283:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    284:   PetscCall(VecRestoreArray(F, &f));
    285:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    286:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    287:   PetscFunctionReturn(PETSC_SUCCESS);
    288: }

    290: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */

    292: PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    293: {
    294:   FVCtx       *ctx = (FVCtx *)vctx;
    295:   PetscInt     i, j, k, Mx, dof, xs, xm, len_slow;
    296:   PetscReal    hx, cfl_idt = 0;
    297:   PetscScalar *x, *f, *slope;
    298:   Vec          Xloc;
    299:   DM           da;

    301:   PetscFunctionBeginUser;
    302:   PetscCall(TSGetDM(ts, &da));
    303:   PetscCall(DMGetLocalVector(da, &Xloc));
    304:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    305:   hx = (ctx->xmax - ctx->xmin) / Mx;
    306:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    307:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    308:   PetscCall(VecZeroEntries(F));
    309:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    310:   PetscCall(VecGetArray(F, &f));
    311:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    312:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    313:   PetscCall(ISGetSize(ctx->iss, &len_slow));

    315:   if (ctx->bctype == FVBC_OUTFLOW) {
    316:     for (i = xs - 2; i < 0; i++) {
    317:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    318:     }
    319:     for (i = Mx; i < xs + xm + 2; i++) {
    320:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    321:     }
    322:   }
    323:   for (i = xs - 1; i < xs + xm + 1; i++) {
    324:     struct _LimitInfo info;
    325:     PetscScalar      *cjmpL, *cjmpR;
    326:     if (i > len_slow - 2) {
    327:       PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
    328:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    329:       cjmpL = &ctx->cjmpLR[0];
    330:       cjmpR = &ctx->cjmpLR[dof];
    331:       for (j = 0; j < dof; j++) {
    332:         PetscScalar jmpL, jmpR;
    333:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    334:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    335:         for (k = 0; k < dof; k++) {
    336:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    337:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    338:         }
    339:       }
    340:       /* Apply limiter to the left and right characteristic jumps */
    341:       info.m  = dof;
    342:       info.hx = hx;
    343:       (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    344:       for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    345:       for (j = 0; j < dof; j++) {
    346:         PetscScalar tmp = 0;
    347:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    348:         slope[i * dof + j] = tmp;
    349:       }
    350:     }
    351:   }

    353:   for (i = xs; i < xs + xm + 1; i++) {
    354:     PetscReal    maxspeed;
    355:     PetscScalar *uL, *uR;
    356:     if (i > len_slow) {
    357:       uL = &ctx->uLR[0];
    358:       uR = &ctx->uLR[dof];
    359:       for (j = 0; j < dof; j++) {
    360:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    361:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    362:       }
    363:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    364:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    365:       if (i > xs) {
    366:         for (j = 0; j < dof; j++) f[(i - len_slow - 1) * dof + j] -= ctx->flux[j] / hx;
    367:       }
    368:       if (i < xs + xm) {
    369:         for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx;
    370:       }
    371:     }
    372:     if (i == len_slow) {
    373:       uL = &ctx->uLR[0];
    374:       uR = &ctx->uLR[dof];
    375:       for (j = 0; j < dof; j++) {
    376:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    377:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    378:       }
    379:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    380:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    381:       if (i < xs + xm) {
    382:         for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx;
    383:       }
    384:     }
    385:   }
    386:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    387:   PetscCall(VecRestoreArray(F, &f));
    388:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    389:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }

    393: PetscErrorCode FVRHSFunctionslow2(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    394: {
    395:   FVCtx       *ctx = (FVCtx *)vctx;
    396:   PetscInt     i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2;
    397:   PetscReal    hx, cfl_idt = 0;
    398:   PetscScalar *x, *f, *slope;
    399:   Vec          Xloc;
    400:   DM           da;

    402:   PetscFunctionBeginUser;
    403:   PetscCall(TSGetDM(ts, &da));
    404:   PetscCall(DMGetLocalVector(da, &Xloc));
    405:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    406:   hx = (ctx->xmax - ctx->xmin) / Mx;
    407:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    408:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    409:   PetscCall(VecZeroEntries(F));
    410:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    411:   PetscCall(VecGetArray(F, &f));
    412:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    413:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    414:   PetscCall(ISGetSize(ctx->iss, &len_slow1));
    415:   PetscCall(ISGetSize(ctx->iss2, &len_slow2));
    416:   if (ctx->bctype == FVBC_OUTFLOW) {
    417:     for (i = xs - 2; i < 0; i++) {
    418:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    419:     }
    420:     for (i = Mx; i < xs + xm + 2; i++) {
    421:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    422:     }
    423:   }
    424:   for (i = xs - 1; i < xs + xm + 1; i++) {
    425:     struct _LimitInfo info;
    426:     PetscScalar      *cjmpL, *cjmpR;
    427:     if (i < len_slow1 + len_slow2 + 1 && i > len_slow1 - 2) {
    428:       PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
    429:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    430:       cjmpL = &ctx->cjmpLR[0];
    431:       cjmpR = &ctx->cjmpLR[dof];
    432:       for (j = 0; j < dof; j++) {
    433:         PetscScalar jmpL, jmpR;
    434:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    435:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    436:         for (k = 0; k < dof; k++) {
    437:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    438:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    439:         }
    440:       }
    441:       /* Apply limiter to the left and right characteristic jumps */
    442:       info.m  = dof;
    443:       info.hx = hx;
    444:       (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    445:       for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    446:       for (j = 0; j < dof; j++) {
    447:         PetscScalar tmp = 0;
    448:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    449:         slope[i * dof + j] = tmp;
    450:       }
    451:     }
    452:   }

    454:   for (i = xs; i < xs + xm + 1; i++) {
    455:     PetscReal    maxspeed;
    456:     PetscScalar *uL, *uR;
    457:     if (i < len_slow1 + len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */
    458:       uL = &ctx->uLR[0];
    459:       uR = &ctx->uLR[dof];
    460:       for (j = 0; j < dof; j++) {
    461:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    462:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    463:       }
    464:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    465:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    466:       if (i > xs) {
    467:         for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx;
    468:       }
    469:       if (i < xs + xm) {
    470:         for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx;
    471:       }
    472:     }
    473:     if (i == len_slow1 + len_slow2) { /* slow parts can be changed based on a */
    474:       uL = &ctx->uLR[0];
    475:       uR = &ctx->uLR[dof];
    476:       for (j = 0; j < dof; j++) {
    477:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    478:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    479:       }
    480:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    481:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    482:       if (i > xs) {
    483:         for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx;
    484:       }
    485:     }
    486:     if (i == len_slow1) { /* slow parts can be changed based on a */
    487:       uL = &ctx->uLR[0];
    488:       uR = &ctx->uLR[dof];
    489:       for (j = 0; j < dof; j++) {
    490:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    491:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    492:       }
    493:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    494:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    495:       if (i < xs + xm) {
    496:         for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx;
    497:       }
    498:     }
    499:   }

    501:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    502:   PetscCall(VecRestoreArray(F, &f));
    503:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    504:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    505:   PetscFunctionReturn(PETSC_SUCCESS);
    506: }

    508: PetscErrorCode FVRHSFunctionfast2(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    509: {
    510:   FVCtx       *ctx = (FVCtx *)vctx;
    511:   PetscInt     i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2;
    512:   PetscReal    hx, cfl_idt = 0;
    513:   PetscScalar *x, *f, *slope;
    514:   Vec          Xloc;
    515:   DM           da;

    517:   PetscFunctionBeginUser;
    518:   PetscCall(TSGetDM(ts, &da));
    519:   PetscCall(DMGetLocalVector(da, &Xloc));
    520:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    521:   hx = (ctx->xmax - ctx->xmin) / Mx;
    522:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    523:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    524:   PetscCall(VecZeroEntries(F));
    525:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    526:   PetscCall(VecGetArray(F, &f));
    527:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    528:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    529:   PetscCall(ISGetSize(ctx->iss, &len_slow1));
    530:   PetscCall(ISGetSize(ctx->iss2, &len_slow2));

    532:   if (ctx->bctype == FVBC_OUTFLOW) {
    533:     for (i = xs - 2; i < 0; i++) {
    534:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    535:     }
    536:     for (i = Mx; i < xs + xm + 2; i++) {
    537:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    538:     }
    539:   }
    540:   for (i = xs - 1; i < xs + xm + 1; i++) {
    541:     struct _LimitInfo info;
    542:     PetscScalar      *cjmpL, *cjmpR;
    543:     if (i > len_slow1 + len_slow2 - 2) {
    544:       PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
    545:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    546:       cjmpL = &ctx->cjmpLR[0];
    547:       cjmpR = &ctx->cjmpLR[dof];
    548:       for (j = 0; j < dof; j++) {
    549:         PetscScalar jmpL, jmpR;
    550:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    551:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    552:         for (k = 0; k < dof; k++) {
    553:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    554:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    555:         }
    556:       }
    557:       /* Apply limiter to the left and right characteristic jumps */
    558:       info.m  = dof;
    559:       info.hx = hx;
    560:       (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
    561:       for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
    562:       for (j = 0; j < dof; j++) {
    563:         PetscScalar tmp = 0;
    564:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    565:         slope[i * dof + j] = tmp;
    566:       }
    567:     }
    568:   }

    570:   for (i = xs; i < xs + xm + 1; i++) {
    571:     PetscReal    maxspeed;
    572:     PetscScalar *uL, *uR;
    573:     if (i > len_slow1 + len_slow2) {
    574:       uL = &ctx->uLR[0];
    575:       uR = &ctx->uLR[dof];
    576:       for (j = 0; j < dof; j++) {
    577:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    578:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    579:       }
    580:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    581:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    582:       if (i > xs) {
    583:         for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2 - 1) * dof + j] -= ctx->flux[j] / hx;
    584:       }
    585:       if (i < xs + xm) {
    586:         for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx;
    587:       }
    588:     }
    589:     if (i == len_slow1 + len_slow2) {
    590:       uL = &ctx->uLR[0];
    591:       uR = &ctx->uLR[dof];
    592:       for (j = 0; j < dof; j++) {
    593:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
    594:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
    595:       }
    596:       PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
    597:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
    598:       if (i < xs + xm) {
    599:         for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx;
    600:       }
    601:     }
    602:   }
    603:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    604:   PetscCall(VecRestoreArray(F, &f));
    605:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    606:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    607:   PetscFunctionReturn(PETSC_SUCCESS);
    608: }

    610: int main(int argc, char *argv[])
    611: {
    612:   char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
    613:   PetscFunctionList limiters = 0, physics = 0;
    614:   MPI_Comm          comm;
    615:   TS                ts;
    616:   DM                da;
    617:   Vec               X, X0, R;
    618:   FVCtx             ctx;
    619:   PetscInt          i, k, dof, xs, xm, Mx, draw = 0, *index_slow, *index_fast, islow = 0, ifast = 0;
    620:   PetscBool         view_final = PETSC_FALSE;
    621:   PetscReal         ptime;

    623:   PetscFunctionBeginUser;
    624:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
    625:   comm = PETSC_COMM_WORLD;
    626:   PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

    628:   /* Register limiters to be available on the command line */
    629:   PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind));
    630:   PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff));
    631:   PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming));
    632:   PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm));
    633:   PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod));
    634:   PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee));
    635:   PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit_MC));
    636:   PetscCall(PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer));
    637:   PetscCall(PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada));
    638:   PetscCall(PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD));
    639:   PetscCall(PetscFunctionListAdd(&limiters, "koren", Limit_Koren));
    640:   PetscCall(PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym));
    641:   PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3));
    642:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2));
    643:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1));
    644:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1));
    645:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10));
    646:   PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100));

    648:   /* Register physical models to be available on the command line */
    649:   PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));

    651:   ctx.comm   = comm;
    652:   ctx.cfl    = 0.9;
    653:   ctx.bctype = FVBC_PERIODIC;
    654:   ctx.xmin   = -1.0;
    655:   ctx.xmax   = 1.0;
    656:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
    657:   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
    658:   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
    659:   PetscCall(PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
    660:   PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
    661:   PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
    662:   PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
    663:   PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
    664:   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
    665:   PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
    666:   PetscCall(PetscOptionsBool("-recursive_split", "Split the domain recursively", "", ctx.recursive, &ctx.recursive, NULL));
    667:   PetscOptionsEnd();

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

    673:   /* Choose the physics from the list of registered models */
    674:   {
    675:     PetscErrorCode (*r)(FVCtx *);
    676:     PetscCall(PetscFunctionListFind(physics, physname, &r));
    677:     PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
    678:     /* Create the physics, will set the number of fields and their names */
    679:     PetscCall((*r)(&ctx));
    680:   }

    682:   /* Create a DMDA to manage the parallel grid */
    683:   PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
    684:   PetscCall(DMSetFromOptions(da));
    685:   PetscCall(DMSetUp(da));
    686:   /* Inform the DMDA of the field names provided by the physics. */
    687:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
    688:   for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
    689:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    690:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

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

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

    699:   /* Create a vector to store the solution and to save the initial state */
    700:   PetscCall(DMCreateGlobalVector(da, &X));
    701:   PetscCall(VecDuplicate(X, &X0));
    702:   PetscCall(VecDuplicate(X, &R));

    704:   /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
    705:   PetscCall(PetscMalloc1(xm * dof, &index_slow));
    706:   PetscCall(PetscMalloc1(xm * dof, &index_fast));
    707:   for (i = xs; i < xs + xm; i++) {
    708:     if (ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx < 0)
    709:       for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    710:     else
    711:       for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
    712:   } /* this step need to be changed based on discontinuous point of a */
    713:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
    714:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));

    716:   /* Create a time-stepping object */
    717:   PetscCall(TSCreate(comm, &ts));
    718:   PetscCall(TSSetDM(ts, da));
    719:   PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
    720:   PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
    721:   PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
    722:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
    723:   PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));

    725:   if (ctx.recursive) {
    726:     TS subts;
    727:     islow = 0;
    728:     ifast = 0;
    729:     for (i = xs; i < xs + xm; i++) {
    730:       PetscReal coord = ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx;
    731:       if (coord >= 0 && coord < ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.)
    732:         for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    733:       if (coord >= ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.)
    734:         for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
    735:     } /* this step need to be changed based on discontinuous point of a */
    736:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss2));
    737:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf2));

    739:     PetscCall(TSRHSSplitGetSubTS(ts, "fast", &subts));
    740:     PetscCall(TSRHSSplitSetIS(subts, "slow", ctx.iss2));
    741:     PetscCall(TSRHSSplitSetIS(subts, "fast", ctx.isf2));
    742:     PetscCall(TSRHSSplitSetRHSFunction(subts, "slow", NULL, FVRHSFunctionslow2, &ctx));
    743:     PetscCall(TSRHSSplitSetRHSFunction(subts, "fast", NULL, FVRHSFunctionfast2, &ctx));
    744:   }

    746:   /*PetscCall(TSSetType(ts,TSSSP));*/
    747:   PetscCall(TSSetType(ts, TSMPRK));
    748:   PetscCall(TSSetMaxTime(ts, 10));
    749:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    751:   /* Compute initial conditions and starting time step */
    752:   PetscCall(FVSample(&ctx, da, 0, X0));
    753:   PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
    754:   PetscCall(VecCopy(X0, X));                            /* The function value was not used so we set X=X0 again */
    755:   PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
    756:   PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
    757:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    758:   {
    759:     PetscInt    steps;
    760:     PetscScalar mass_initial, mass_final, mass_difference;

    762:     PetscCall(TSSolve(ts, X));
    763:     PetscCall(TSGetSolveTime(ts, &ptime));
    764:     PetscCall(TSGetStepNumber(ts, &steps));
    765:     PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    766:     /* calculate the total mass at initial time and final time */
    767:     mass_initial = 0.0;
    768:     mass_final   = 0.0;
    769:     PetscCall(VecSum(X0, &mass_initial));
    770:     PetscCall(VecSum(X, &mass_final));
    771:     mass_difference = (ctx.xmax - ctx.xmin) / (PetscScalar)Mx * (mass_final - mass_initial);
    772:     PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
    773:     if (ctx.simulation) {
    774:       PetscViewer fd;
    775:       char        filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
    776:       Vec         XR;
    777:       PetscReal   nrm1, nrmsup;
    778:       PetscBool   flg;

    780:       PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
    781:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
    782:       PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
    783:       PetscCall(VecDuplicate(X0, &XR));
    784:       PetscCall(VecLoad(XR, fd));
    785:       PetscCall(PetscViewerDestroy(&fd));
    786:       PetscCall(VecAYPX(XR, -1, X));
    787:       PetscCall(VecNorm(XR, NORM_1, &nrm1));
    788:       PetscCall(VecNorm(XR, NORM_INFINITY, &nrmsup));
    789:       nrm1 /= Mx;
    790:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g  ||x-x_e||_sup %g\n", (double)nrm1, (double)nrmsup));
    791:       PetscCall(VecDestroy(&XR));
    792:     }
    793:   }

    795:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    796:   if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
    797:   if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
    798:   if (draw & 0x4) {
    799:     Vec Y;
    800:     PetscCall(VecDuplicate(X, &Y));
    801:     PetscCall(FVSample(&ctx, da, ptime, Y));
    802:     PetscCall(VecAYPX(Y, -1, X));
    803:     PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    804:     PetscCall(VecDestroy(&Y));
    805:   }

    807:   if (view_final) {
    808:     PetscViewer viewer;
    809:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    810:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    811:     PetscCall(VecView(X, viewer));
    812:     PetscCall(PetscViewerPopFormat(viewer));
    813:     PetscCall(PetscViewerDestroy(&viewer));
    814:   }

    816:   /* Clean up */
    817:   PetscCall((*ctx.physics.destroy)(ctx.physics.user));
    818:   for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
    819:   PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
    820:   PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
    821:   PetscCall(ISDestroy(&ctx.iss));
    822:   PetscCall(ISDestroy(&ctx.isf));
    823:   PetscCall(ISDestroy(&ctx.iss2));
    824:   PetscCall(ISDestroy(&ctx.isf2));
    825:   PetscCall(VecDestroy(&X));
    826:   PetscCall(VecDestroy(&X0));
    827:   PetscCall(VecDestroy(&R));
    828:   PetscCall(DMDestroy(&da));
    829:   PetscCall(TSDestroy(&ts));
    830:   PetscCall(PetscFree(index_slow));
    831:   PetscCall(PetscFree(index_fast));
    832:   PetscCall(PetscFunctionListDestroy(&limiters));
    833:   PetscCall(PetscFunctionListDestroy(&physics));
    834:   PetscCall(PetscFinalize());
    835:   return 0;
    836: }

    838: /*TEST
    839:     build:
    840:       requires: !complex
    841:       depends: finitevolume1d.c

    843:     test:
    844:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0

    846:     test:
    847:       suffix: 2
    848:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
    849:       output_file: output/ex5_1.out

    851:     test:
    852:       suffix: 3
    853:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0

    855:     test:
    856:       suffix: 4
    857:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
    858:       output_file: output/ex5_3.out

    860:     test:
    861:       suffix: 5
    862:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split

    864:     test:
    865:       suffix: 6
    866:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split
    867: TEST*/