Actual source code: ex6.c

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

      7: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
      8:                            "  advection   - Constant coefficient scalar advection\n"
      9:                            "                u_t       + (a*u)_x               = 0\n"
     10:                            "  for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
     11:                            "                hxs  = hratio*hxf \n"
     12:                            "  where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
     13:                            "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
     14:                            "                the states across shocks and rarefactions\n"
     15:                            "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
     16:                            "                also the reference solution should be generated by user and stored in a binary file.\n"
     17:                            "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
     18:                            "Several initial conditions can be chosen with -initial N\n\n"
     19:                            "The problem size should be set with -da_grid_x M\n\n";

     21: #include <petscts.h>
     22: #include <petscdm.h>
     23: #include <petscdmda.h>
     24: #include <petscdraw.h>
     25: #include "finitevolume1d.h"

     27: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
     28: {
     29:   PetscReal range = xmax - xmin;
     30:   return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
     31: }

     33: /* --------------------------------- Advection ----------------------------------- */
     34: typedef struct {
     35:   PetscReal a; /* advective velocity */
     36: } AdvectCtx;

     38: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
     39: {
     40:   AdvectCtx *ctx = (AdvectCtx *)vctx;
     41:   PetscReal  speed;

     43:   PetscFunctionBeginUser;
     44:   speed     = ctx->a;
     45:   flux[0]   = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
     46:   *maxspeed = speed;
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

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

     54:   PetscFunctionBeginUser;
     55:   X[0]      = 1.;
     56:   Xi[0]     = 1.;
     57:   speeds[0] = ctx->a;
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

     61: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
     62: {
     63:   AdvectCtx *ctx = (AdvectCtx *)vctx;
     64:   PetscReal  a   = ctx->a, x0;

     66:   PetscFunctionBeginUser;
     67:   switch (bctype) {
     68:   case FVBC_OUTFLOW:
     69:     x0 = x - a * t;
     70:     break;
     71:   case FVBC_PERIODIC:
     72:     x0 = RangeMod(x - a * t, xmin, xmax);
     73:     break;
     74:   default:
     75:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
     76:   }
     77:   switch (initial) {
     78:   case 0:
     79:     u[0] = (x0 < 0) ? 1 : -1;
     80:     break;
     81:   case 1:
     82:     u[0] = (x0 < 0) ? -1 : 1;
     83:     break;
     84:   case 2:
     85:     u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
     86:     break;
     87:   case 3:
     88:     u[0] = PetscSinReal(2 * PETSC_PI * x0);
     89:     break;
     90:   case 4:
     91:     u[0] = PetscAbs(x0);
     92:     break;
     93:   case 5:
     94:     u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
     95:     break;
     96:   case 6:
     97:     u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
     98:     break;
     99:   case 7:
    100:     u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    101:     break;
    102:   default:
    103:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    104:   }
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
    109: {
    110:   AdvectCtx *user;

    112:   PetscFunctionBeginUser;
    113:   PetscCall(PetscNew(&user));
    114:   ctx->physics2.sample2         = PhysicsSample_Advect;
    115:   ctx->physics2.riemann2        = PhysicsRiemann_Advect;
    116:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
    117:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
    118:   ctx->physics2.user            = user;
    119:   ctx->physics2.dof             = 1;
    120:   PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0]));
    121:   user->a = 1;
    122:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    123:   {
    124:     PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
    125:   }
    126:   PetscOptionsEnd();
    127:   PetscFunctionReturn(PETSC_SUCCESS);
    128: }

    130: PetscErrorCode FVSample_2WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U)
    131: {
    132:   PetscScalar   *u, *uj, xj, xi;
    133:   PetscInt       i, j, k, dof, xs, xm, Mx;
    134:   const PetscInt N = 200;
    135:   PetscReal      hs, hf;

    137:   PetscFunctionBeginUser;
    138:   PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
    139:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    140:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    141:   PetscCall(DMDAVecGetArray(da, U, &u));
    142:   PetscCall(PetscMalloc1(dof, &uj));
    143:   hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    144:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    145:   for (i = xs; i < xs + xm; i++) {
    146:     if (i < ctx->sf) {
    147:       xi = ctx->xmin + 0.5 * hs + i * hs;
    148:       /* Integrate over cell i using trapezoid rule with N points. */
    149:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    150:       for (j = 0; j < N + 1; j++) {
    151:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    152:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    153:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    154:       }
    155:     } else if (i < ctx->fs) {
    156:       xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf;
    157:       /* Integrate over cell i using trapezoid rule with N points. */
    158:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    159:       for (j = 0; j < N + 1; j++) {
    160:         xj = xi + hf * (j - N / 2) / (PetscReal)N;
    161:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    162:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    163:       }
    164:     } else {
    165:       xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs;
    166:       /* Integrate over cell i using trapezoid rule with N points. */
    167:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    168:       for (j = 0; j < N + 1; j++) {
    169:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    170:         PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    171:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    172:       }
    173:     }
    174:   }
    175:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    176:   PetscCall(PetscFree(uj));
    177:   PetscFunctionReturn(PETSC_SUCCESS);
    178: }

    180: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
    181: {
    182:   Vec                Y;
    183:   PetscInt           i, Mx;
    184:   const PetscScalar *ptr_X, *ptr_Y;
    185:   PetscReal          hs, hf;

    187:   PetscFunctionBeginUser;
    188:   PetscCall(VecGetSize(X, &Mx));
    189:   PetscCall(VecDuplicate(X, &Y));
    190:   PetscCall(FVSample_2WaySplit(ctx, da, t, Y));
    191:   hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    192:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    193:   PetscCall(VecGetArrayRead(X, &ptr_X));
    194:   PetscCall(VecGetArrayRead(Y, &ptr_Y));
    195:   for (i = 0; i < Mx; i++) {
    196:     if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
    197:     else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
    198:   }
    199:   PetscCall(VecRestoreArrayRead(X, &ptr_X));
    200:   PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
    201:   PetscCall(VecDestroy(&Y));
    202:   PetscFunctionReturn(PETSC_SUCCESS);
    203: }

    205: PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    206: {
    207:   FVCtx       *ctx = (FVCtx *)vctx;
    208:   PetscInt     i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
    209:   PetscReal    hxf, hxs, cfl_idt = 0;
    210:   PetscScalar *x, *f, *slope;
    211:   Vec          Xloc;
    212:   DM           da;

    214:   PetscFunctionBeginUser;
    215:   PetscCall(TSGetDM(ts, &da));
    216:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
    217:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
    218:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    219:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    220:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
    221:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

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

    225:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    226:   PetscCall(DMDAVecGetArray(da, F, &f));
    227:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points                                           */

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

    231:   if (ctx->bctype == FVBC_OUTFLOW) {
    232:     for (i = xs - 2; i < 0; i++) {
    233:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    234:     }
    235:     for (i = Mx; i < xs + xm + 2; i++) {
    236:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    237:     }
    238:   }
    239:   for (i = xs - 1; i < xs + xm + 1; i++) {
    240:     struct _LimitInfo info;
    241:     PetscScalar      *cjmpL, *cjmpR;
    242:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    243:     PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    244:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    245:     PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    246:     cjmpL = &ctx->cjmpLR[0];
    247:     cjmpR = &ctx->cjmpLR[dof];
    248:     for (j = 0; j < dof; j++) {
    249:       PetscScalar jmpL, jmpR;
    250:       jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    251:       jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    252:       for (k = 0; k < dof; k++) {
    253:         cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    254:         cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    255:       }
    256:     }
    257:     /* Apply limiter to the left and right characteristic jumps */
    258:     info.m   = dof;
    259:     info.hxs = hxs;
    260:     info.hxf = hxf;
    261:     (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    262:     for (j = 0; j < dof; j++) {
    263:       PetscScalar tmp = 0;
    264:       for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    265:       slope[i * dof + j] = tmp;
    266:     }
    267:   }

    269:   for (i = xs; i < xs + xm + 1; i++) {
    270:     PetscReal    maxspeed;
    271:     PetscScalar *uL, *uR;
    272:     uL = &ctx->uLR[0];
    273:     uR = &ctx->uLR[dof];
    274:     if (i < sf) { /* slow region */
    275:       for (j = 0; j < dof; j++) {
    276:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    277:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    278:       }
    279:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    280:       if (i > xs) {
    281:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    282:       }
    283:       if (i < xs + xm) {
    284:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    285:       }
    286:     } else if (i == sf) { /* interface between the slow region and the fast region */
    287:       for (j = 0; j < dof; j++) {
    288:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    289:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    290:       }
    291:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    292:       if (i > xs) {
    293:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    294:       }
    295:       if (i < xs + xm) {
    296:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
    297:       }
    298:     } else if (i > sf && i < fs) { /* fast region */
    299:       for (j = 0; j < dof; j++) {
    300:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    301:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    302:       }
    303:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    304:       if (i > xs) {
    305:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
    306:       }
    307:       if (i < xs + xm) {
    308:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
    309:       }
    310:     } else if (i == fs) { /* interface between the fast region and the slow region */
    311:       for (j = 0; j < dof; j++) {
    312:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    313:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    314:       }
    315:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    316:       if (i > xs) {
    317:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
    318:       }
    319:       if (i < xs + xm) {
    320:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    321:       }
    322:     } else { /* slow region */
    323:       for (j = 0; j < dof; j++) {
    324:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    325:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    326:       }
    327:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    328:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
    329:       if (i > xs) {
    330:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
    331:       }
    332:       if (i < xs + xm) {
    333:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
    334:       }
    335:     }
    336:   }
    337:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    338:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    339:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    340:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    341:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    342:   if (0) {
    343:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
    344:     PetscReal dt, tnow;
    345:     PetscCall(TSGetTimeStep(ts, &dt));
    346:     PetscCall(TSGetTime(ts, &tnow));
    347:     if (dt > 0.5 / ctx->cfl_idt) {
    348:       if (1) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(1 / (2 * ctx->cfl_idt))));
    349:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
    350:     }
    351:   }
    352:   PetscFunctionReturn(PETSC_SUCCESS);
    353: }

    355: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
    356: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    357: {
    358:   FVCtx       *ctx = (FVCtx *)vctx;
    359:   PetscInt     i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
    360:   PetscReal    hxs, hxf, cfl_idt = 0;
    361:   PetscScalar *x, *f, *slope;
    362:   Vec          Xloc;
    363:   DM           da;

    365:   PetscFunctionBeginUser;
    366:   PetscCall(TSGetDM(ts, &da));
    367:   PetscCall(DMGetLocalVector(da, &Xloc));
    368:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    369:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    370:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    371:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    372:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    373:   PetscCall(VecZeroEntries(F));
    374:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    375:   PetscCall(VecGetArray(F, &f));
    376:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    377:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    379:   if (ctx->bctype == FVBC_OUTFLOW) {
    380:     for (i = xs - 2; i < 0; i++) {
    381:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    382:     }
    383:     for (i = Mx; i < xs + xm + 2; i++) {
    384:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    385:     }
    386:   }
    387:   for (i = xs - 1; i < xs + xm + 1; i++) {
    388:     struct _LimitInfo info;
    389:     PetscScalar      *cjmpL, *cjmpR;
    390:     if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
    391:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    392:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    393:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    394:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    395:       cjmpL = &ctx->cjmpLR[0];
    396:       cjmpR = &ctx->cjmpLR[dof];
    397:       for (j = 0; j < dof; j++) {
    398:         PetscScalar jmpL, jmpR;
    399:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    400:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    401:         for (k = 0; k < dof; k++) {
    402:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    403:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    404:         }
    405:       }
    406:       /* Apply limiter to the left and right characteristic jumps */
    407:       info.m   = dof;
    408:       info.hxs = hxs;
    409:       info.hxf = hxf;
    410:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    411:       for (j = 0; j < dof; j++) {
    412:         PetscScalar tmp = 0;
    413:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    414:         slope[i * dof + j] = tmp;
    415:       }
    416:     }
    417:   }

    419:   for (i = xs; i < xs + xm + 1; i++) {
    420:     PetscReal    maxspeed;
    421:     PetscScalar *uL, *uR;
    422:     uL = &ctx->uLR[0];
    423:     uR = &ctx->uLR[dof];
    424:     if (i < sf - lsbwidth) { /* slow region */
    425:       for (j = 0; j < dof; j++) {
    426:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    427:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    428:       }
    429:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    430:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
    431:       if (i > xs) {
    432:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    433:       }
    434:       if (i < xs + xm) {
    435:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    436:         islow++;
    437:       }
    438:     }
    439:     if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
    440:       for (j = 0; j < dof; j++) {
    441:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    442:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    443:       }
    444:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    445:       if (i > xs) {
    446:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    447:       }
    448:     }
    449:     if (i == fs + rsbwidth) { /* slow region */
    450:       for (j = 0; j < dof; j++) {
    451:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    452:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    453:       }
    454:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    455:       if (i < xs + xm) {
    456:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    457:         islow++;
    458:       }
    459:     }
    460:     if (i > fs + rsbwidth) { /* slow region */
    461:       for (j = 0; j < dof; j++) {
    462:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    463:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    464:       }
    465:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    466:       if (i > xs) {
    467:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    468:       }
    469:       if (i < xs + xm) {
    470:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    471:         islow++;
    472:       }
    473:     }
    474:   }
    475:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    476:   PetscCall(VecRestoreArray(F, &f));
    477:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    478:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    479:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    480:   PetscFunctionReturn(PETSC_SUCCESS);
    481: }

    483: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    484: {
    485:   FVCtx       *ctx = (FVCtx *)vctx;
    486:   PetscInt     i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
    487:   PetscReal    hxs, hxf;
    488:   PetscScalar *x, *f, *slope;
    489:   Vec          Xloc;
    490:   DM           da;

    492:   PetscFunctionBeginUser;
    493:   PetscCall(TSGetDM(ts, &da));
    494:   PetscCall(DMGetLocalVector(da, &Xloc));
    495:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    496:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    497:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    498:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    499:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    500:   PetscCall(VecZeroEntries(F));
    501:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    502:   PetscCall(VecGetArray(F, &f));
    503:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    504:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    506:   if (ctx->bctype == FVBC_OUTFLOW) {
    507:     for (i = xs - 2; i < 0; i++) {
    508:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    509:     }
    510:     for (i = Mx; i < xs + xm + 2; i++) {
    511:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    512:     }
    513:   }
    514:   for (i = xs - 1; i < xs + xm + 1; i++) {
    515:     struct _LimitInfo info;
    516:     PetscScalar      *cjmpL, *cjmpR;
    517:     if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
    518:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
    519:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    520:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    521:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    522:       cjmpL = &ctx->cjmpLR[0];
    523:       cjmpR = &ctx->cjmpLR[dof];
    524:       for (j = 0; j < dof; j++) {
    525:         PetscScalar jmpL, jmpR;
    526:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    527:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    528:         for (k = 0; k < dof; k++) {
    529:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    530:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    531:         }
    532:       }
    533:       /* Apply limiter to the left and right characteristic jumps */
    534:       info.m   = dof;
    535:       info.hxs = hxs;
    536:       info.hxf = hxf;
    537:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    538:       for (j = 0; j < dof; j++) {
    539:         PetscScalar tmp = 0;
    540:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    541:         slope[i * dof + j] = tmp;
    542:       }
    543:     }
    544:   }

    546:   for (i = xs; i < xs + xm + 1; i++) {
    547:     PetscReal    maxspeed;
    548:     PetscScalar *uL, *uR;
    549:     uL = &ctx->uLR[0];
    550:     uR = &ctx->uLR[dof];
    551:     if (i == sf - lsbwidth) {
    552:       for (j = 0; j < dof; j++) {
    553:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    554:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    555:       }
    556:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    557:       if (i < xs + xm) {
    558:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    559:         islow++;
    560:       }
    561:     }
    562:     if (i > sf - lsbwidth && i < sf) {
    563:       for (j = 0; j < dof; j++) {
    564:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    565:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    566:       }
    567:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    568:       if (i > xs) {
    569:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    570:       }
    571:       if (i < xs + xm) {
    572:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    573:         islow++;
    574:       }
    575:     }
    576:     if (i == sf) { /* interface between the slow region and the fast region */
    577:       for (j = 0; j < dof; j++) {
    578:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    579:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    580:       }
    581:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    582:       if (i > xs) {
    583:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    584:       }
    585:     }
    586:     if (i == fs) { /* interface between the fast region and the slow region */
    587:       for (j = 0; j < dof; j++) {
    588:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    589:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    590:       }
    591:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    592:       if (i < xs + xm) {
    593:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    594:         islow++;
    595:       }
    596:     }
    597:     if (i > fs && i < fs + rsbwidth) {
    598:       for (j = 0; j < dof; j++) {
    599:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    600:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    601:       }
    602:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    603:       if (i > xs) {
    604:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    605:       }
    606:       if (i < xs + xm) {
    607:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
    608:         islow++;
    609:       }
    610:     }
    611:     if (i == fs + rsbwidth) {
    612:       for (j = 0; j < dof; j++) {
    613:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    614:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    615:       }
    616:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    617:       if (i > xs) {
    618:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
    619:       }
    620:     }
    621:   }
    622:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    623:   PetscCall(VecRestoreArray(F, &f));
    624:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    625:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    626:   PetscFunctionReturn(PETSC_SUCCESS);
    627: }

    629: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
    630: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    631: {
    632:   FVCtx       *ctx = (FVCtx *)vctx;
    633:   PetscInt     i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
    634:   PetscReal    hxs, hxf;
    635:   PetscScalar *x, *f, *slope;
    636:   Vec          Xloc;
    637:   DM           da;

    639:   PetscFunctionBeginUser;
    640:   PetscCall(TSGetDM(ts, &da));
    641:   PetscCall(DMGetLocalVector(da, &Xloc));
    642:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    643:   hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
    644:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
    645:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    646:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    647:   PetscCall(VecZeroEntries(F));
    648:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    649:   PetscCall(VecGetArray(F, &f));
    650:   PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
    651:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

    653:   if (ctx->bctype == FVBC_OUTFLOW) {
    654:     for (i = xs - 2; i < 0; i++) {
    655:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    656:     }
    657:     for (i = Mx; i < xs + xm + 2; i++) {
    658:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    659:     }
    660:   }
    661:   for (i = xs - 1; i < xs + xm + 1; i++) {
    662:     struct _LimitInfo info;
    663:     PetscScalar      *cjmpL, *cjmpR;
    664:     if (i > sf - 2 && i < fs + 1) {
    665:       PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
    666:       PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
    667:       cjmpL = &ctx->cjmpLR[0];
    668:       cjmpR = &ctx->cjmpLR[dof];
    669:       for (j = 0; j < dof; j++) {
    670:         PetscScalar jmpL, jmpR;
    671:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
    672:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
    673:         for (k = 0; k < dof; k++) {
    674:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
    675:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
    676:         }
    677:       }
    678:       /* Apply limiter to the left and right characteristic jumps */
    679:       info.m   = dof;
    680:       info.hxs = hxs;
    681:       info.hxf = hxf;
    682:       (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
    683:       for (j = 0; j < dof; j++) {
    684:         PetscScalar tmp = 0;
    685:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
    686:         slope[i * dof + j] = tmp;
    687:       }
    688:     }
    689:   }

    691:   for (i = xs; i < xs + xm + 1; i++) {
    692:     PetscReal    maxspeed;
    693:     PetscScalar *uL, *uR;
    694:     uL = &ctx->uLR[0];
    695:     uR = &ctx->uLR[dof];
    696:     if (i == sf) { /* interface between the slow region and the fast region */
    697:       for (j = 0; j < dof; j++) {
    698:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
    699:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    700:       }
    701:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    702:       if (i < xs + xm) {
    703:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
    704:         ifast++;
    705:       }
    706:     }
    707:     if (i > sf && i < fs) { /* fast region */
    708:       for (j = 0; j < dof; j++) {
    709:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    710:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
    711:       }
    712:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    713:       if (i > xs) {
    714:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
    715:       }
    716:       if (i < xs + xm) {
    717:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
    718:         ifast++;
    719:       }
    720:     }
    721:     if (i == fs) { /* interface between the fast region and the slow region */
    722:       for (j = 0; j < dof; j++) {
    723:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
    724:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
    725:       }
    726:       PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
    727:       if (i > xs) {
    728:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
    729:       }
    730:     }
    731:   }
    732:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    733:   PetscCall(VecRestoreArray(F, &f));
    734:   PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
    735:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    736:   PetscFunctionReturn(PETSC_SUCCESS);
    737: }

    739: int main(int argc, char *argv[])
    740: {
    741:   char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
    742:   PetscFunctionList limiters = 0, physics = 0;
    743:   MPI_Comm          comm;
    744:   TS                ts;
    745:   DM                da;
    746:   Vec               X, X0, R;
    747:   FVCtx             ctx;
    748:   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;
    749:   PetscBool         view_final = PETSC_FALSE;
    750:   PetscReal         ptime;

    752:   PetscFunctionBeginUser;
    753:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
    754:   comm = PETSC_COMM_WORLD;
    755:   PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

    757:   /* Register limiters to be available on the command line */
    758:   PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
    759:   PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
    760:   PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
    761:   PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
    762:   PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
    763:   PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
    764:   PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
    765:   PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));

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

    770:   ctx.comm   = comm;
    771:   ctx.cfl    = 0.9;
    772:   ctx.bctype = FVBC_PERIODIC;
    773:   ctx.xmin   = -1.0;
    774:   ctx.xmax   = 1.0;
    775:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
    776:   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
    777:   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
    778:   PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
    779:   PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
    780:   PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
    781:   PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
    782:   PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
    783:   PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
    784:   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
    785:   PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
    786:   PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
    787:   PetscOptionsEnd();

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

    793:   /* Choose the physics from the list of registered models */
    794:   {
    795:     PetscErrorCode (*r)(FVCtx *);
    796:     PetscCall(PetscFunctionListFind(physics, physname, &r));
    797:     PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
    798:     /* Create the physics, will set the number of fields and their names */
    799:     PetscCall((*r)(&ctx));
    800:   }

    802:   /* Create a DMDA to manage the parallel grid */
    803:   PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
    804:   PetscCall(DMSetFromOptions(da));
    805:   PetscCall(DMSetUp(da));
    806:   /* Inform the DMDA of the field names provided by the physics. */
    807:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
    808:   for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
    809:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    810:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

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

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

    819:   /* Create a vector to store the solution and to save the initial state */
    820:   PetscCall(DMCreateGlobalVector(da, &X));
    821:   PetscCall(VecDuplicate(X, &X0));
    822:   PetscCall(VecDuplicate(X, &R));

    824:   /* create index for slow parts and fast parts,
    825:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
    826:   count_slow = Mx / (1.0 + ctx.hratio / 3.0);
    827:   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+hartio/3) is even");
    828:   count_fast = Mx - count_slow;
    829:   ctx.sf     = count_slow / 2;
    830:   ctx.fs     = ctx.sf + count_fast;
    831:   PetscCall(PetscMalloc1(xm * dof, &index_slow));
    832:   PetscCall(PetscMalloc1(xm * dof, &index_fast));
    833:   PetscCall(PetscMalloc1(6 * dof, &index_slowbuffer));
    834:   if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
    835:     ctx.lsbwidth = 2;
    836:     ctx.rsbwidth = 4;
    837:   } else {
    838:     ctx.lsbwidth = 4;
    839:     ctx.rsbwidth = 2;
    840:   }
    841:   for (i = xs; i < xs + xm; i++) {
    842:     if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
    843:       for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    844:     else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
    845:       for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
    846:     else
    847:       for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
    848:   }
    849:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
    850:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
    851:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));

    853:   /* Create a time-stepping object */
    854:   PetscCall(TSCreate(comm, &ts));
    855:   PetscCall(TSSetDM(ts, da));
    856:   PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
    857:   PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
    858:   PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
    859:   PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
    860:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
    861:   PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
    862:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));

    864:   PetscCall(TSSetType(ts, TSSSP));
    865:   /*PetscCall(TSSetType(ts,TSMPRK));*/
    866:   PetscCall(TSSetMaxTime(ts, 10));
    867:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    869:   /* Compute initial conditions and starting time step */
    870:   PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
    871:   PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
    872:   PetscCall(VecCopy(X0, X));                                      /* The function value was not used so we set X=X0 again */
    873:   PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
    874:   PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
    875:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    876:   {
    877:     PetscInt           steps;
    878:     PetscScalar        mass_initial, mass_final, mass_difference, mass_differenceg;
    879:     const PetscScalar *ptr_X, *ptr_X0;
    880:     const PetscReal    hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow;
    881:     const PetscReal    hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;

    883:     PetscCall(TSSolve(ts, X));
    884:     PetscCall(TSGetSolveTime(ts, &ptime));
    885:     PetscCall(TSGetStepNumber(ts, &steps));
    886:     /* calculate the total mass at initial time and final time */
    887:     mass_initial = 0.0;
    888:     mass_final   = 0.0;
    889:     PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
    890:     PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
    891:     for (i = xs; i < xs + xm; i++) {
    892:       if (i < ctx.sf || i > ctx.fs - 1) {
    893:         for (k = 0; k < dof; k++) {
    894:           mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
    895:           mass_final   = mass_final + hs * ptr_X[i * dof + k];
    896:         }
    897:       } else {
    898:         for (k = 0; k < dof; k++) {
    899:           mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
    900:           mass_final   = mass_final + hf * ptr_X[i * dof + k];
    901:         }
    902:       }
    903:     }
    904:     PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
    905:     PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
    906:     mass_difference = mass_final - mass_initial;
    907:     PetscCallMPI(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
    908:     PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
    909:     PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    910:     PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1 / ctx.cfl_idt)));
    911:     if (ctx.exact) {
    912:       PetscReal nrm1 = 0;
    913:       PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1));
    914:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    915:     }
    916:     if (ctx.simulation) {
    917:       PetscReal          nrm1 = 0;
    918:       PetscViewer        fd;
    919:       char               filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
    920:       Vec                XR;
    921:       PetscBool          flg;
    922:       const PetscScalar *ptr_XR;
    923:       PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
    924:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
    925:       PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
    926:       PetscCall(VecDuplicate(X0, &XR));
    927:       PetscCall(VecLoad(XR, fd));
    928:       PetscCall(PetscViewerDestroy(&fd));
    929:       PetscCall(VecGetArrayRead(X, &ptr_X));
    930:       PetscCall(VecGetArrayRead(XR, &ptr_XR));
    931:       for (i = xs; i < xs + xm; i++) {
    932:         if (i < ctx.sf || i > ctx.fs - 1)
    933:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
    934:         else
    935:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
    936:       }
    937:       PetscCall(VecRestoreArrayRead(X, &ptr_X));
    938:       PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
    939:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    940:       PetscCall(VecDestroy(&XR));
    941:     }
    942:   }

    944:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    945:   if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
    946:   if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
    947:   if (draw & 0x4) {
    948:     Vec Y;
    949:     PetscCall(VecDuplicate(X, &Y));
    950:     PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y));
    951:     PetscCall(VecAYPX(Y, -1, X));
    952:     PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    953:     PetscCall(VecDestroy(&Y));
    954:   }

    956:   if (view_final) {
    957:     PetscViewer viewer;
    958:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    959:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    960:     PetscCall(VecView(X, viewer));
    961:     PetscCall(PetscViewerPopFormat(viewer));
    962:     PetscCall(PetscViewerDestroy(&viewer));
    963:   }

    965:   /* Clean up */
    966:   PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
    967:   for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
    968:   PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
    969:   PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
    970:   PetscCall(VecDestroy(&X));
    971:   PetscCall(VecDestroy(&X0));
    972:   PetscCall(VecDestroy(&R));
    973:   PetscCall(DMDestroy(&da));
    974:   PetscCall(TSDestroy(&ts));
    975:   PetscCall(ISDestroy(&ctx.iss));
    976:   PetscCall(ISDestroy(&ctx.isf));
    977:   PetscCall(ISDestroy(&ctx.issb));
    978:   PetscCall(PetscFree(index_slow));
    979:   PetscCall(PetscFree(index_fast));
    980:   PetscCall(PetscFree(index_slowbuffer));
    981:   PetscCall(PetscFunctionListDestroy(&limiters));
    982:   PetscCall(PetscFunctionListDestroy(&physics));
    983:   PetscCall(PetscFinalize());
    984:   return 0;
    985: }

    987: /*TEST

    989:     build:
    990:       requires: !complex
    991:       depends: finitevolume1d.c

    993:     test:
    994:       suffix: 1
    995:       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

    997:     test:
    998:       suffix: 2
    999:       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 0
    1000:       output_file: output/ex6_1.out

    1002: TEST*/