Actual source code: ex7.c

      1: static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter with semidiscrete time stepping.\n"
      2:                            "  advection   - Constant coefficient scalar advection\n"
      3:                            "                u_t       + (a*u)_x               = 0\n"
      4:                            "  for this toy problem, we choose different meshsizes for different sub-domains, say\n"
      5:                            "                hxs  = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n"
      6:                            "                hxf  = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n"
      7:                            "  with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of coarse\n\n"
      8:                            "  grids and fine grids is hratio.\n"
      9:                            "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
     10:                            "                the states across shocks and rarefactions\n"
     11:                            "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
     12:                            "                also the reference solution should be generated by user and stored in a binary file.\n"
     13:                            "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
     14:                            "Several initial conditions can be chosen with -initial N\n\n"
     15:                            "The problem size should be set with -da_grid_x M\n\n"
     16:                            "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n"
     17:                            "                             u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t))                 \n"
     18:                            "                     limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2))))    \n"
     19:                            "                             r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1)))                                   \n"
     20:                            "                             alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1))                     \n"
     21:                            "                             gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1))                 \n";

     23: #include <petscts.h>
     24: #include <petscdm.h>
     25: #include <petscdmda.h>
     26: #include <petscdraw.h>
     27: #include <petscmath.h>

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

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

     37: typedef enum {
     38:   FVBC_PERIODIC,
     39:   FVBC_OUTFLOW
     40: } FVBCType;
     41: static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", 0};

     43: typedef struct {
     44:   PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *);
     45:   PetscErrorCode (*flux)(void *, const PetscScalar *, PetscScalar *, PetscReal *);
     46:   PetscErrorCode (*destroy)(void *);
     47:   void    *user;
     48:   PetscInt dof;
     49:   char    *fieldname[16];
     50: } PhysicsCtx;

     52: typedef struct {
     53:   PhysicsCtx physics;
     54:   MPI_Comm   comm;
     55:   char       prefix[256];

     57:   /* Local work arrays */
     58:   PetscScalar *flux;   /* Flux across interface                                                      */
     59:   PetscReal   *speeds; /* Speeds of each wave                                                        */
     60:   PetscScalar *u;      /* value at face                                                              */

     62:   PetscReal cfl_idt; /* Max allowable value of 1/Delta t                                           */
     63:   PetscReal cfl;
     64:   PetscReal xmin, xmax;
     65:   PetscInt  initial;
     66:   PetscBool exact;
     67:   PetscBool simulation;
     68:   FVBCType  bctype;
     69:   PetscInt  hratio; /* hratio = hslow/hfast */
     70:   IS        isf, iss;
     71:   PetscInt  sf, fs; /* slow-fast and fast-slow interfaces */
     72: } FVCtx;

     74: /* --------------------------------- Physics ----------------------------------- */
     75: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
     76: {
     77:   PetscFunctionBeginUser;
     78:   PetscCall(PetscFree(vctx));
     79:   PetscFunctionReturn(PETSC_SUCCESS);
     80: }

     82: /* --------------------------------- Advection ----------------------------------- */
     83: typedef struct {
     84:   PetscReal a; /* advective velocity */
     85: } AdvectCtx;

     87: static PetscErrorCode PhysicsFlux_Advect(void *vctx, const PetscScalar *u, PetscScalar *flux, PetscReal *maxspeed)
     88: {
     89:   AdvectCtx *ctx = (AdvectCtx *)vctx;
     90:   PetscReal  speed;

     92:   PetscFunctionBeginUser;
     93:   speed     = ctx->a;
     94:   flux[0]   = speed * u[0];
     95:   *maxspeed = speed;
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
    100: {
    101:   AdvectCtx *ctx = (AdvectCtx *)vctx;
    102:   PetscReal  a   = ctx->a, x0;

    104:   PetscFunctionBeginUser;
    105:   switch (bctype) {
    106:   case FVBC_OUTFLOW:
    107:     x0 = x - a * t;
    108:     break;
    109:   case FVBC_PERIODIC:
    110:     x0 = RangeMod(x - a * t, xmin, xmax);
    111:     break;
    112:   default:
    113:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
    114:   }
    115:   switch (initial) {
    116:   case 0:
    117:     u[0] = (x0 < 0) ? 1 : -1;
    118:     break;
    119:   case 1:
    120:     u[0] = (x0 < 0) ? -1 : 1;
    121:     break;
    122:   case 2:
    123:     u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
    124:     break;
    125:   case 3:
    126:     u[0] = PetscSinReal(2 * PETSC_PI * x0);
    127:     break;
    128:   case 4:
    129:     u[0] = PetscAbs(x0);
    130:     break;
    131:   case 5:
    132:     u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
    133:     break;
    134:   case 6:
    135:     u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
    136:     break;
    137:   case 7:
    138:     u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
    139:     break;
    140:   default:
    141:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
    142:   }
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

    146: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
    147: {
    148:   AdvectCtx *user;

    150:   PetscFunctionBeginUser;
    151:   PetscCall(PetscNew(&user));
    152:   ctx->physics.sample  = PhysicsSample_Advect;
    153:   ctx->physics.flux    = PhysicsFlux_Advect;
    154:   ctx->physics.destroy = PhysicsDestroy_SimpleFree;
    155:   ctx->physics.user    = user;
    156:   ctx->physics.dof     = 1;
    157:   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
    158:   user->a = 1;
    159:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
    160:   {
    161:     PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
    162:   }
    163:   PetscOptionsEnd();
    164:   PetscFunctionReturn(PETSC_SUCCESS);
    165: }

    167: /* --------------------------------- Finite Volume Solver ----------------------------------- */

    169: static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    170: {
    171:   FVCtx       *ctx = (FVCtx *)vctx;
    172:   PetscInt     i, j, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
    173:   PetscReal    hf, hs, cfl_idt = 0;
    174:   PetscScalar *x, *f, *r, *min, *alpha, *gamma;
    175:   Vec          Xloc;
    176:   DM           da;

    178:   PetscFunctionBeginUser;
    179:   PetscCall(TSGetDM(ts, &da));
    180:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
    181:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
    182:   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
    183:   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
    184:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
    185:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    186:   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
    187:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    188:   PetscCall(DMDAVecGetArray(da, F, &f));
    189:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    190:   PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

    192:   if (ctx->bctype == FVBC_OUTFLOW) {
    193:     for (i = xs - 2; i < 0; i++) {
    194:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    195:     }
    196:     for (i = Mx; i < xs + xm + 2; i++) {
    197:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    198:     }
    199:   }

    201:   for (i = xs; i < xs + xm + 1; i++) {
    202:     PetscReal    maxspeed;
    203:     PetscScalar *u;
    204:     if (i < sf || i > fs + 1) {
    205:       u        = &ctx->u[0];
    206:       alpha[0] = 1.0 / 6.0;
    207:       gamma[0] = 1.0 / 3.0;
    208:       for (j = 0; j < dof; j++) {
    209:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    210:         min[j] = PetscMin(r[j], 2.0);
    211:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    212:       }
    213:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    214:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hs));
    215:       if (i > xs) {
    216:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
    217:       }
    218:       if (i < xs + xm) {
    219:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
    220:       }
    221:     } else if (i == sf) {
    222:       u        = &ctx->u[0];
    223:       alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
    224:       gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
    225:       for (j = 0; j < dof; j++) {
    226:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    227:         min[j] = PetscMin(r[j], 2.0);
    228:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    229:       }
    230:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    231:       if (i > xs) {
    232:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
    233:       }
    234:       if (i < xs + xm) {
    235:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
    236:       }
    237:     } else if (i == sf + 1) {
    238:       u        = &ctx->u[0];
    239:       alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
    240:       gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
    241:       for (j = 0; j < dof; j++) {
    242:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    243:         min[j] = PetscMin(r[j], 2.0);
    244:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    245:       }
    246:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    247:       if (i > xs) {
    248:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
    249:       }
    250:       if (i < xs + xm) {
    251:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
    252:       }
    253:     } else if (i > sf + 1 && i < fs) {
    254:       u        = &ctx->u[0];
    255:       alpha[0] = 1.0 / 6.0;
    256:       gamma[0] = 1.0 / 3.0;
    257:       for (j = 0; j < dof; j++) {
    258:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    259:         min[j] = PetscMin(r[j], 2.0);
    260:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    261:       }
    262:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    263:       if (i > xs) {
    264:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
    265:       }
    266:       if (i < xs + xm) {
    267:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
    268:       }
    269:     } else if (i == fs) {
    270:       u        = &ctx->u[0];
    271:       alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
    272:       gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
    273:       for (j = 0; j < dof; j++) {
    274:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    275:         min[j] = PetscMin(r[j], 2.0);
    276:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    277:       }
    278:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    279:       if (i > xs) {
    280:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
    281:       }
    282:       if (i < xs + xm) {
    283:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
    284:       }
    285:     } else if (i == fs + 1) {
    286:       u        = &ctx->u[0];
    287:       alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
    288:       gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
    289:       for (j = 0; j < dof; j++) {
    290:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    291:         min[j] = PetscMin(r[j], 2.0);
    292:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    293:       }
    294:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    295:       if (i > xs) {
    296:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
    297:       }
    298:       if (i < xs + xm) {
    299:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
    300:       }
    301:     }
    302:   }
    303:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    304:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    305:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    306:   PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
    307:   if (0) {
    308:     /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
    309:     PetscReal dt, tnow;
    310:     PetscCall(TSGetTimeStep(ts, &dt));
    311:     PetscCall(TSGetTime(ts, &tnow));
    312:     if (dt > 0.5 / ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
    313:   }
    314:   PetscCall(PetscFree4(r, min, alpha, gamma));
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: static PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    319: {
    320:   FVCtx       *ctx = (FVCtx *)vctx;
    321:   PetscInt     i, j, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs;
    322:   PetscReal    hf, hs;
    323:   PetscScalar *x, *f, *r, *min, *alpha, *gamma;
    324:   Vec          Xloc;
    325:   DM           da;

    327:   PetscFunctionBeginUser;
    328:   PetscCall(TSGetDM(ts, &da));
    329:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
    330:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
    331:   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
    332:   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
    333:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
    334:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    335:   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
    336:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    337:   PetscCall(VecGetArray(F, &f));
    338:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    339:   PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

    341:   if (ctx->bctype == FVBC_OUTFLOW) {
    342:     for (i = xs - 2; i < 0; i++) {
    343:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    344:     }
    345:     for (i = Mx; i < xs + xm + 2; i++) {
    346:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    347:     }
    348:   }

    350:   for (i = xs; i < xs + xm + 1; i++) {
    351:     PetscReal    maxspeed;
    352:     PetscScalar *u;
    353:     if (i < sf) {
    354:       u        = &ctx->u[0];
    355:       alpha[0] = 1.0 / 6.0;
    356:       gamma[0] = 1.0 / 3.0;
    357:       for (j = 0; j < dof; j++) {
    358:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    359:         min[j] = PetscMin(r[j], 2.0);
    360:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    361:       }
    362:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    363:       if (i > xs) {
    364:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
    365:       }
    366:       if (i < xs + xm) {
    367:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
    368:         islow++;
    369:       }
    370:     } else if (i == sf) {
    371:       u        = &ctx->u[0];
    372:       alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
    373:       gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
    374:       for (j = 0; j < dof; j++) {
    375:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    376:         min[j] = PetscMin(r[j], 2.0);
    377:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    378:       }
    379:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    380:       if (i > xs) {
    381:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
    382:       }
    383:     } else if (i == fs) {
    384:       u        = &ctx->u[0];
    385:       alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
    386:       gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
    387:       for (j = 0; j < dof; j++) {
    388:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    389:         min[j] = PetscMin(r[j], 2.0);
    390:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    391:       }
    392:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    393:       if (i < xs + xm) {
    394:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
    395:         islow++;
    396:       }
    397:     } else if (i == fs + 1) {
    398:       u        = &ctx->u[0];
    399:       alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
    400:       gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
    401:       for (j = 0; j < dof; j++) {
    402:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    403:         min[j] = PetscMin(r[j], 2.0);
    404:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    405:       }
    406:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    407:       if (i > xs) {
    408:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
    409:       }
    410:       if (i < xs + xm) {
    411:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
    412:         islow++;
    413:       }
    414:     } else if (i > fs + 1) {
    415:       u        = &ctx->u[0];
    416:       alpha[0] = 1.0 / 6.0;
    417:       gamma[0] = 1.0 / 3.0;
    418:       for (j = 0; j < dof; j++) {
    419:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    420:         min[j] = PetscMin(r[j], 2.0);
    421:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    422:       }
    423:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    424:       if (i > xs) {
    425:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
    426:       }
    427:       if (i < xs + xm) {
    428:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
    429:         islow++;
    430:       }
    431:     }
    432:   }
    433:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    434:   PetscCall(VecRestoreArray(F, &f));
    435:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    436:   PetscCall(PetscFree4(r, min, alpha, gamma));
    437:   PetscFunctionReturn(PETSC_SUCCESS);
    438: }

    440: static PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
    441: {
    442:   FVCtx       *ctx = (FVCtx *)vctx;
    443:   PetscInt     i, j, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
    444:   PetscReal    hf, hs;
    445:   PetscScalar *x, *f, *r, *min, *alpha, *gamma;
    446:   Vec          Xloc;
    447:   DM           da;

    449:   PetscFunctionBeginUser;
    450:   PetscCall(TSGetDM(ts, &da));
    451:   PetscCall(DMGetLocalVector(da, &Xloc));                                 /* Xloc contains ghost points                                     */
    452:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points                              */
    453:   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
    454:   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
    455:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points       */
    456:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    457:   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
    458:   PetscCall(DMDAVecGetArray(da, Xloc, &x));
    459:   PetscCall(VecGetArray(F, &f));
    460:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    461:   PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));

    463:   if (ctx->bctype == FVBC_OUTFLOW) {
    464:     for (i = xs - 2; i < 0; i++) {
    465:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
    466:     }
    467:     for (i = Mx; i < xs + xm + 2; i++) {
    468:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
    469:     }
    470:   }

    472:   for (i = xs; i < xs + xm + 1; i++) {
    473:     PetscReal    maxspeed;
    474:     PetscScalar *u;
    475:     if (i == sf) {
    476:       u        = &ctx->u[0];
    477:       alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
    478:       gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
    479:       for (j = 0; j < dof; j++) {
    480:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    481:         min[j] = PetscMin(r[j], 2.0);
    482:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    483:       }
    484:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    485:       if (i < xs + xm) {
    486:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
    487:         ifast++;
    488:       }
    489:     } else if (i == sf + 1) {
    490:       u        = &ctx->u[0];
    491:       alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
    492:       gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
    493:       for (j = 0; j < dof; j++) {
    494:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    495:         min[j] = PetscMin(r[j], 2.0);
    496:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    497:       }
    498:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    499:       if (i > xs) {
    500:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
    501:       }
    502:       if (i < xs + xm) {
    503:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
    504:         ifast++;
    505:       }
    506:     } else if (i > sf + 1 && i < fs) {
    507:       u        = &ctx->u[0];
    508:       alpha[0] = 1.0 / 6.0;
    509:       gamma[0] = 1.0 / 3.0;
    510:       for (j = 0; j < dof; j++) {
    511:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    512:         min[j] = PetscMin(r[j], 2.0);
    513:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    514:       }
    515:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    516:       if (i > xs) {
    517:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
    518:       }
    519:       if (i < xs + xm) {
    520:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
    521:         ifast++;
    522:       }
    523:     } else if (i == fs) {
    524:       u        = &ctx->u[0];
    525:       alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
    526:       gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
    527:       for (j = 0; j < dof; j++) {
    528:         r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    529:         min[j] = PetscMin(r[j], 2.0);
    530:         u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
    531:       }
    532:       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
    533:       if (i > xs) {
    534:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
    535:       }
    536:     }
    537:   }
    538:   PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
    539:   PetscCall(VecRestoreArray(F, &f));
    540:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    541:   PetscCall(PetscFree4(r, min, alpha, gamma));
    542:   PetscFunctionReturn(PETSC_SUCCESS);
    543: }

    545: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */

    547: PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
    548: {
    549:   PetscScalar   *u, *uj, xj, xi;
    550:   PetscInt       i, j, k, dof, xs, xm, Mx, count_slow, count_fast;
    551:   const PetscInt N = 200;

    553:   PetscFunctionBeginUser;
    554:   PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
    555:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    556:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    557:   PetscCall(DMDAVecGetArray(da, U, &u));
    558:   PetscCall(PetscMalloc1(dof, &uj));
    559:   const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
    560:   const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
    561:   count_slow         = Mx / (1 + ctx->hratio);
    562:   count_fast         = Mx - count_slow;
    563:   for (i = xs; i < xs + xm; i++) {
    564:     if (i * hs + 0.5 * hs < (ctx->xmax - ctx->xmin) * 0.25) {
    565:       xi = ctx->xmin + 0.5 * hs + i * hs;
    566:       /* Integrate over cell i using trapezoid rule with N points. */
    567:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    568:       for (j = 0; j < N + 1; j++) {
    569:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    570:         PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    571:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    572:       }
    573:     } else if ((ctx->xmax - ctx->xmin) * 0.25 + (i - count_slow / 2) * hf + 0.5 * hf < (ctx->xmax - ctx->xmin) * 0.75) {
    574:       xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.25 + 0.5 * hf + (i - count_slow / 2) * hf;
    575:       /* Integrate over cell i using trapezoid rule with N points. */
    576:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    577:       for (j = 0; j < N + 1; j++) {
    578:         xj = xi + hf * (j - N / 2) / (PetscReal)N;
    579:         PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    580:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    581:       }
    582:     } else {
    583:       xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.75 + 0.5 * hs + (i - count_slow / 2 - count_fast) * hs;
    584:       /* Integrate over cell i using trapezoid rule with N points. */
    585:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
    586:       for (j = 0; j < N + 1; j++) {
    587:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
    588:         PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
    589:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
    590:       }
    591:     }
    592:   }
    593:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    594:   PetscCall(PetscFree(uj));
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

    598: static PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
    599: {
    600:   PetscReal          xmin, xmax;
    601:   PetscScalar        sum, tvsum, tvgsum;
    602:   const PetscScalar *x;
    603:   PetscInt           imin, imax, Mx, i, j, xs, xm, dof;
    604:   Vec                Xloc;
    605:   PetscBool          iascii;

    607:   PetscFunctionBeginUser;
    608:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    609:   if (iascii) {
    610:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
    611:     PetscCall(DMGetLocalVector(da, &Xloc));
    612:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    613:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
    614:     PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
    615:     PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
    616:     PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    617:     tvsum = 0;
    618:     for (i = xs; i < xs + xm; i++) {
    619:       for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
    620:     }
    621:     PetscCallMPI(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da)));
    622:     PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
    623:     PetscCall(DMRestoreLocalVector(da, &Xloc));

    625:     PetscCall(VecMin(X, &imin, &xmin));
    626:     PetscCall(VecMax(X, &imax, &xmax));
    627:     PetscCall(VecSum(X, &sum));
    628:     PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%g,%g] with minimum at %" PetscInt_FMT ", mean %g, ||x||_TV %g\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvgsum / Mx)));
    629:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
    634: {
    635:   Vec                Y;
    636:   PetscInt           i, Mx, count_slow = 0, count_fast = 0;
    637:   const PetscScalar *ptr_X, *ptr_Y;

    639:   PetscFunctionBeginUser;
    640:   PetscCall(VecGetSize(X, &Mx));
    641:   PetscCall(VecDuplicate(X, &Y));
    642:   PetscCall(FVSample(ctx, da, t, Y));
    643:   const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
    644:   const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
    645:   count_slow         = (PetscReal)Mx / (1.0 + ctx->hratio);
    646:   count_fast         = Mx - count_slow;
    647:   PetscCall(VecGetArrayRead(X, &ptr_X));
    648:   PetscCall(VecGetArrayRead(Y, &ptr_Y));
    649:   for (i = 0; i < Mx; i++) {
    650:     if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
    651:     else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
    652:   }
    653:   PetscCall(VecRestoreArrayRead(X, &ptr_X));
    654:   PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
    655:   PetscCall(VecDestroy(&Y));
    656:   PetscFunctionReturn(PETSC_SUCCESS);
    657: }

    659: int main(int argc, char *argv[])
    660: {
    661:   char              physname[256] = "advect", final_fname[256] = "solution.m";
    662:   PetscFunctionList physics = 0;
    663:   MPI_Comm          comm;
    664:   TS                ts;
    665:   DM                da;
    666:   Vec               X, X0, R;
    667:   FVCtx             ctx;
    668:   PetscInt          i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, *index_slow, *index_fast;
    669:   PetscBool         view_final = PETSC_FALSE;
    670:   PetscReal         ptime;

    672:   PetscFunctionBeginUser;
    673:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
    674:   comm = PETSC_COMM_WORLD;
    675:   PetscCall(PetscMemzero(&ctx, sizeof(ctx)));

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

    680:   ctx.comm   = comm;
    681:   ctx.cfl    = 0.9;
    682:   ctx.bctype = FVBC_PERIODIC;
    683:   ctx.xmin   = -1.0;
    684:   ctx.xmax   = 1.0;
    685:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
    686:   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
    687:   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
    688:   PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
    689:   PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
    690:   PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
    691:   PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
    692:   PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
    693:   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
    694:   PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
    695:   PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
    696:   PetscOptionsEnd();

    698:   /* Choose the physics from the list of registered models */
    699:   {
    700:     PetscErrorCode (*r)(FVCtx *);
    701:     PetscCall(PetscFunctionListFind(physics, physname, &r));
    702:     PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
    703:     /* Create the physics, will set the number of fields and their names */
    704:     PetscCall((*r)(&ctx));
    705:   }

    707:   /* Create a DMDA to manage the parallel grid */
    708:   PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
    709:   PetscCall(DMSetFromOptions(da));
    710:   PetscCall(DMSetUp(da));
    711:   /* Inform the DMDA of the field names provided by the physics. */
    712:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
    713:   for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
    714:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
    715:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

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

    720:   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
    721:   PetscCall(PetscMalloc3(dof, &ctx.u, dof, &ctx.flux, dof, &ctx.speeds));

    723:   /* Create a vector to store the solution and to save the initial state */
    724:   PetscCall(DMCreateGlobalVector(da, &X));
    725:   PetscCall(VecDuplicate(X, &X0));
    726:   PetscCall(VecDuplicate(X, &R));

    728:   /* create index for slow parts and fast parts*/
    729:   count_slow = Mx / (1 + ctx.hratio);
    730:   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) is even");
    731:   count_fast = Mx - count_slow;
    732:   ctx.sf     = count_slow / 2;
    733:   ctx.fs     = ctx.sf + count_fast;
    734:   PetscCall(PetscMalloc1(xm * dof, &index_slow));
    735:   PetscCall(PetscMalloc1(xm * dof, &index_fast));
    736:   for (i = xs; i < xs + xm; i++) {
    737:     if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1)
    738:       for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
    739:     else
    740:       for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
    741:   }
    742:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
    743:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));

    745:   /* Create a time-stepping object */
    746:   PetscCall(TSCreate(comm, &ts));
    747:   PetscCall(TSSetDM(ts, da));
    748:   PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
    749:   PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
    750:   PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
    751:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
    752:   PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));

    754:   PetscCall(TSSetType(ts, TSMPRK));
    755:   PetscCall(TSSetMaxTime(ts, 10));
    756:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    758:   /* Compute initial conditions and starting time step */
    759:   PetscCall(FVSample(&ctx, da, 0, X0));
    760:   PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
    761:   PetscCall(VecCopy(X0, X));                            /* The function value was not used so we set X=X0 again */
    762:   PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
    763:   PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
    764:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    765:   {
    766:     PetscInt           steps;
    767:     PetscScalar        mass_initial, mass_final, mass_difference, mass_differenceg;
    768:     const PetscScalar *ptr_X, *ptr_X0;
    769:     const PetscReal    hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow;
    770:     const PetscReal    hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast;
    771:     PetscCall(TSSolve(ts, X));
    772:     PetscCall(TSGetSolveTime(ts, &ptime));
    773:     PetscCall(TSGetStepNumber(ts, &steps));
    774:     /* calculate the total mass at initial time and final time */
    775:     mass_initial = 0.0;
    776:     mass_final   = 0.0;
    777:     PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
    778:     PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
    779:     for (i = xs; i < xs + xm; i++) {
    780:       if (i < ctx.sf || i > ctx.fs - 1) {
    781:         for (k = 0; k < dof; k++) {
    782:           mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
    783:           mass_final   = mass_final + hs * ptr_X[i * dof + k];
    784:         }
    785:       } else {
    786:         for (k = 0; k < dof; k++) {
    787:           mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
    788:           mass_final   = mass_final + hf * ptr_X[i * dof + k];
    789:         }
    790:       }
    791:     }
    792:     PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
    793:     PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
    794:     mass_difference = mass_final - mass_initial;
    795:     PetscCallMPI(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
    796:     PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
    797:     PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
    798:     if (ctx.exact) {
    799:       PetscReal nrm1 = 0;
    800:       PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1));
    801:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    802:     }
    803:     if (ctx.simulation) {
    804:       PetscReal          nrm1 = 0;
    805:       PetscViewer        fd;
    806:       char               filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
    807:       Vec                XR;
    808:       PetscBool          flg;
    809:       const PetscScalar *ptr_XR;
    810:       PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
    811:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
    812:       PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
    813:       PetscCall(VecDuplicate(X0, &XR));
    814:       PetscCall(VecLoad(XR, fd));
    815:       PetscCall(PetscViewerDestroy(&fd));
    816:       PetscCall(VecGetArrayRead(X, &ptr_X));
    817:       PetscCall(VecGetArrayRead(XR, &ptr_XR));
    818:       for (i = 0; i < Mx; i++) {
    819:         if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] - ptr_XR[i]);
    820:         else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]);
    821:       }
    822:       PetscCall(VecRestoreArrayRead(X, &ptr_X));
    823:       PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
    824:       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
    825:       PetscCall(VecDestroy(&XR));
    826:     }
    827:   }

    829:   PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
    830:   if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
    831:   if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
    832:   if (draw & 0x4) {
    833:     Vec Y;
    834:     PetscCall(VecDuplicate(X, &Y));
    835:     PetscCall(FVSample(&ctx, da, ptime, Y));
    836:     PetscCall(VecAYPX(Y, -1, X));
    837:     PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
    838:     PetscCall(VecDestroy(&Y));
    839:   }

    841:   if (view_final) {
    842:     PetscViewer viewer;
    843:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
    844:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
    845:     PetscCall(VecView(X, viewer));
    846:     PetscCall(PetscViewerPopFormat(viewer));
    847:     PetscCall(PetscViewerDestroy(&viewer));
    848:   }

    850:   /* Clean up */
    851:   PetscCall((*ctx.physics.destroy)(ctx.physics.user));
    852:   for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
    853:   PetscCall(PetscFree3(ctx.u, ctx.flux, ctx.speeds));
    854:   PetscCall(ISDestroy(&ctx.iss));
    855:   PetscCall(ISDestroy(&ctx.isf));
    856:   PetscCall(VecDestroy(&X));
    857:   PetscCall(VecDestroy(&X0));
    858:   PetscCall(VecDestroy(&R));
    859:   PetscCall(DMDestroy(&da));
    860:   PetscCall(TSDestroy(&ts));
    861:   PetscCall(PetscFree(index_slow));
    862:   PetscCall(PetscFree(index_fast));
    863:   PetscCall(PetscFunctionListDestroy(&physics));
    864:   PetscCall(PetscFinalize());
    865:   return 0;
    866: }

    868: /*TEST

    870:     build:
    871:       requires: !complex

    873:     test:
    874:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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

    876:     test:
    877:       suffix: 2
    878:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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
    879:       output_file: output/ex7_1.out

    881:     test:
    882:       suffix: 3
    883:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0

    885:     test:
    886:       suffix: 4
    887:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
    888:       output_file: output/ex7_3.out

    890:     test:
    891:       suffix: 5
    892:       nsize: 2
    893:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
    894:       output_file: output/ex7_3.out
    895: TEST*/