Actual source code: ex8.c

  1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form 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 (slow-medium-fast-medium-slow), \n"
  5:                            "  the meshsize ratio between two adjacient sub-domains is controlled with -hratio,\n"
  6:                            "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
  7:                            "                the states across shocks and rarefactions\n"
  8:                            "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
  9:                            "                also the reference solution should be generated by user and stored in a binary file.\n"
 10:                            "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 11:                            "Several initial conditions can be chosen with -initial N\n\n"
 12:                            "The problem size should be set with -da_grid_x M\n\n";

 14: #include <petscts.h>
 15: #include <petscdm.h>
 16: #include <petscdmda.h>
 17: #include <petscdraw.h>
 18: #include "finitevolume1d.h"

 20: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
 21: {
 22:   PetscReal range = xmax - xmin;
 23:   return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
 24: }

 26: /* --------------------------------- Advection ----------------------------------- */
 27: typedef struct {
 28:   PetscReal a; /* advective velocity */
 29: } AdvectCtx;

 31: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
 32: {
 33:   AdvectCtx *ctx = (AdvectCtx *)vctx;
 34:   PetscReal  speed;

 37:   speed     = ctx->a;
 38:   flux[0]   = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
 39:   *maxspeed = speed;
 40:   return 0;
 41: }

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

 48:   X[0]      = 1.;
 49:   Xi[0]     = 1.;
 50:   speeds[0] = ctx->a;
 51:   return 0;
 52: }

 54: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
 55: {
 56:   AdvectCtx *ctx = (AdvectCtx *)vctx;
 57:   PetscReal  a   = ctx->a, x0;

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

101: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
102: {
103:   AdvectCtx *user;

106:   PetscNew(&user);
107:   ctx->physics2.sample2         = PhysicsSample_Advect;
108:   ctx->physics2.riemann2        = PhysicsRiemann_Advect;
109:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
110:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
111:   ctx->physics2.user            = user;
112:   ctx->physics2.dof             = 1;
113:   PetscStrallocpy("u", &ctx->physics2.fieldname[0]);
114:   user->a = 1;
115:   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
116:   {
117:     PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL);
118:   }
119:   PetscOptionsEnd();
120:   return 0;
121: }

123: PetscErrorCode FVSample_3WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U)
124: {
125:   PetscScalar   *u, *uj, xj, xi;
126:   PetscInt       i, j, k, dof, xs, xm, Mx;
127:   const PetscInt N = 200;
128:   PetscReal      hs, hm, hf;

132:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
133:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
134:   DMDAVecGetArray(da, U, &u);
135:   PetscMalloc1(dof, &uj);

137:   hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
138:   hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
139:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
140:   for (i = xs; i < xs + xm; i++) {
141:     if (i < ctx->sm) {
142:       xi = ctx->xmin + 0.5 * hs + i * hs;
143:       /* Integrate over cell i using trapezoid rule with N points. */
144:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
145:       for (j = 0; j < N + 1; j++) {
146:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
147:         (*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj);
148:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
149:       }
150:     } else if (i < ctx->mf) {
151:       xi = ctx->xmin + ctx->sm * hs + 0.5 * hm + (i - ctx->sm) * hm;
152:       /* Integrate over cell i using trapezoid rule with N points. */
153:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
154:       for (j = 0; j < N + 1; j++) {
155:         xj = xi + hm * (j - N / 2) / (PetscReal)N;
156:         (*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj);
157:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
158:       }
159:     } else if (i < ctx->fm) {
160:       xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + 0.5 * hf + (i - ctx->mf) * hf;
161:       /* Integrate over cell i using trapezoid rule with N points. */
162:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
163:       for (j = 0; j < N + 1; j++) {
164:         xj = xi + hf * (j - N / 2) / (PetscReal)N;
165:         (*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj);
166:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
167:       }
168:     } else if (i < ctx->ms) {
169:       xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + 0.5 * hm + (i - ctx->fm) * hm;
170:       /* Integrate over cell i using trapezoid rule with N points. */
171:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
172:       for (j = 0; j < N + 1; j++) {
173:         xj = xi + hm * (j - N / 2) / (PetscReal)N;
174:         (*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj);
175:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
176:       }
177:     } else {
178:       xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + (ctx->ms - ctx->fm) * hm + 0.5 * hs + (i - ctx->ms) * hs;
179:       /* Integrate over cell i using trapezoid rule with N points. */
180:       for (k = 0; k < dof; k++) u[i * dof + k] = 0;
181:       for (j = 0; j < N + 1; j++) {
182:         xj = xi + hs * (j - N / 2) / (PetscReal)N;
183:         (*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj);
184:         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
185:       }
186:     }
187:   }
188:   DMDAVecRestoreArray(da, U, &u);
189:   PetscFree(uj);
190:   return 0;
191: }

193: static PetscErrorCode SolutionErrorNorms_3WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
194: {
195:   Vec                Y;
196:   PetscInt           i, Mx;
197:   const PetscScalar *ptr_X, *ptr_Y;
198:   PetscReal          hs, hm, hf;

201:   VecGetSize(X, &Mx);
202:   hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
203:   hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
204:   hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
205:   VecDuplicate(X, &Y);
206:   FVSample_3WaySplit(ctx, da, t, Y);
207:   VecGetArrayRead(X, &ptr_X);
208:   VecGetArrayRead(Y, &ptr_Y);
209:   for (i = 0; i < Mx; i++) {
210:     if (i < ctx->sm || i > ctx->ms - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
211:     else if (i < ctx->mf || i > ctx->fm - 1) *nrm1 += hm * PetscAbs(ptr_X[i] - ptr_Y[i]);
212:     else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
213:   }
214:   VecRestoreArrayRead(X, &ptr_X);
215:   VecRestoreArrayRead(Y, &ptr_Y);
216:   VecDestroy(&Y);
217:   return 0;
218: }

220: PetscErrorCode FVRHSFunction_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
221: {
222:   FVCtx       *ctx = (FVCtx *)vctx;
223:   PetscInt     i, j, k, Mx, dof, xs, xm, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->ms;
224:   PetscReal    hxf, hxm, hxs, cfl_idt = 0;
225:   PetscScalar *x, *f, *slope;
226:   Vec          Xloc;
227:   DM           da;

230:   TSGetDM(ts, &da);
231:   DMGetLocalVector(da, &Xloc);                                 /* Xloc contains ghost points                                     */
232:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0); /* Mx is the number of center points                              */
233:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
234:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
235:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
236:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc); /* X is solution vector which does not contain ghost points       */
237:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);

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

241:   DMDAVecGetArray(da, Xloc, &x);
242:   DMDAVecGetArray(da, F, &f);
243:   DMDAGetArray(da, PETSC_TRUE, &slope); /* contains ghost points                                           */

245:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

247:   if (ctx->bctype == FVBC_OUTFLOW) {
248:     for (i = xs - 2; i < 0; i++) {
249:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
250:     }
251:     for (i = Mx; i < xs + xm + 2; i++) {
252:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
253:     }
254:   }
255:   for (i = xs - 1; i < xs + xm + 1; i++) {
256:     struct _LimitInfo info;
257:     PetscScalar      *cjmpL, *cjmpR;
258:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
259:     (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
260:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
261:     PetscArrayzero(ctx->cjmpLR, 2 * dof);
262:     cjmpL = &ctx->cjmpLR[0];
263:     cjmpR = &ctx->cjmpLR[dof];
264:     for (j = 0; j < dof; j++) {
265:       PetscScalar jmpL, jmpR;
266:       jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
267:       jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
268:       for (k = 0; k < dof; k++) {
269:         cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
270:         cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
271:       }
272:     }
273:     /* Apply limiter to the left and right characteristic jumps */
274:     info.m   = dof;
275:     info.hxs = hxs;
276:     info.hxm = hxm;
277:     info.hxf = hxf;
278:     (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
279:     for (j = 0; j < dof; j++) {
280:       PetscScalar tmp = 0;
281:       for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
282:       slope[i * dof + j] = tmp;
283:     }
284:   }

286:   for (i = xs; i < xs + xm + 1; i++) {
287:     PetscReal    maxspeed;
288:     PetscScalar *uL, *uR;
289:     uL = &ctx->uLR[0];
290:     uR = &ctx->uLR[dof];
291:     if (i < sm || i > ms) { /* slow region */
292:       for (j = 0; j < dof; j++) {
293:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
294:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
295:       }
296:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
297:       if (i > xs) {
298:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
299:       }
300:       if (i < xs + xm) {
301:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
302:       }
303:     } else if (i == sm) { /* interface between slow and medium component */
304:       for (j = 0; j < dof; j++) {
305:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
306:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
307:       }
308:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
309:       if (i > xs) {
310:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
311:       }
312:       if (i < xs + xm) {
313:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
314:       }
315:     } else if (i == ms) { /* interface between medium and slow regions */
316:       for (j = 0; j < dof; j++) {
317:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
318:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
319:       }
320:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
321:       if (i > xs) {
322:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
323:       }
324:       if (i < xs + xm) {
325:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
326:       }
327:     } else if (i < mf || i > fm) { /* medium region */
328:       for (j = 0; j < dof; j++) {
329:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
330:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
331:       }
332:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
333:       if (i > xs) {
334:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
335:       }
336:       if (i < xs + xm) {
337:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
338:       }
339:     } else if (i == mf) { /* interface between medium and fast regions */
340:       for (j = 0; j < dof; j++) {
341:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
342:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
343:       }
344:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
345:       if (i > xs) {
346:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
347:       }
348:       if (i < xs + xm) {
349:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
350:       }
351:     } else if (i == fm) { /* interface between fast and medium regions */
352:       for (j = 0; j < dof; j++) {
353:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
354:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
355:       }
356:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
357:       if (i > xs) {
358:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
359:       }
360:       if (i < xs + xm) {
361:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
362:       }
363:     } else { /* fast region */
364:       for (j = 0; j < dof; j++) {
365:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
366:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
367:       }
368:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
369:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
370:       if (i > xs) {
371:         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
372:       }
373:       if (i < xs + xm) {
374:         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
375:       }
376:     }
377:   }
378:   DMDAVecRestoreArray(da, Xloc, &x);
379:   DMDAVecRestoreArray(da, F, &f);
380:   DMDARestoreArray(da, PETSC_TRUE, &slope);
381:   DMRestoreLocalVector(da, &Xloc);
382:   MPI_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da));
383:   if (0) {
384:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
385:     PetscReal dt, tnow;
386:     TSGetTimeStep(ts, &dt);
387:     TSGetTime(ts, &tnow);
388:     if (dt > 0.5 / ctx->cfl_idt) PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt));
389:   }
390:   return 0;
391: }

393: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
394: PetscErrorCode FVRHSFunctionslow_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
395: {
396:   FVCtx       *ctx = (FVCtx *)vctx;
397:   PetscInt     i, j, k, Mx, dof, xs, xm, islow = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
398:   PetscReal    hxs, hxm, hxf, cfl_idt = 0;
399:   PetscScalar *x, *f, *slope;
400:   Vec          Xloc;
401:   DM           da;

404:   TSGetDM(ts, &da);
405:   DMGetLocalVector(da, &Xloc);
406:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
407:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
408:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
409:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
410:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
411:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
412:   VecZeroEntries(F);
413:   DMDAVecGetArray(da, Xloc, &x);
414:   VecGetArray(F, &f);
415:   DMDAGetArray(da, PETSC_TRUE, &slope);
416:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

418:   if (ctx->bctype == FVBC_OUTFLOW) {
419:     for (i = xs - 2; i < 0; i++) {
420:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
421:     }
422:     for (i = Mx; i < xs + xm + 2; i++) {
423:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
424:     }
425:   }
426:   for (i = xs - 1; i < xs + xm + 1; i++) {
427:     struct _LimitInfo info;
428:     PetscScalar      *cjmpL, *cjmpR;
429:     if (i < sm - lsbwidth + 1 || i > ms + rsbwidth - 2) { /* slow components and the first and last fast components */
430:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
431:       (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
432:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
433:       PetscArrayzero(ctx->cjmpLR, 2 * dof);
434:       cjmpL = &ctx->cjmpLR[0];
435:       cjmpR = &ctx->cjmpLR[dof];
436:       for (j = 0; j < dof; j++) {
437:         PetscScalar jmpL, jmpR;
438:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
439:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
440:         for (k = 0; k < dof; k++) {
441:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
442:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
443:         }
444:       }
445:       /* Apply limiter to the left and right characteristic jumps */
446:       info.m   = dof;
447:       info.hxs = hxs;
448:       info.hxm = hxm;
449:       info.hxf = hxf;
450:       (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
451:       for (j = 0; j < dof; j++) {
452:         PetscScalar tmp = 0;
453:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
454:         slope[i * dof + j] = tmp;
455:       }
456:     }
457:   }

459:   for (i = xs; i < xs + xm + 1; i++) {
460:     PetscReal    maxspeed;
461:     PetscScalar *uL, *uR;
462:     uL = &ctx->uLR[0];
463:     uR = &ctx->uLR[dof];
464:     if (i < sm - lsbwidth) { /* slow region */
465:       for (j = 0; j < dof; j++) {
466:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
467:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
468:       }
469:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
470:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
471:       if (i > xs) {
472:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
473:       }
474:       if (i < xs + xm) {
475:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
476:         islow++;
477:       }
478:     }
479:     if (i == sm - lsbwidth) { /* interface between slow and medium regions */
480:       for (j = 0; j < dof; j++) {
481:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
482:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
483:       }
484:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
485:       if (i > xs) {
486:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
487:       }
488:     }
489:     if (i == ms + rsbwidth) { /* interface between medium and slow regions */
490:       for (j = 0; j < dof; j++) {
491:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
492:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
493:       }
494:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
495:       if (i < xs + xm) {
496:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
497:         islow++;
498:       }
499:     }
500:     if (i > ms + rsbwidth) { /* slow region */
501:       for (j = 0; j < dof; j++) {
502:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
503:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
504:       }
505:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
506:       cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
507:       if (i > xs) {
508:         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
509:       }
510:       if (i < xs + xm) {
511:         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
512:         islow++;
513:       }
514:     }
515:   }
516:   DMDAVecRestoreArray(da, Xloc, &x);
517:   VecRestoreArray(F, &f);
518:   DMDARestoreArray(da, PETSC_TRUE, &slope);
519:   DMRestoreLocalVector(da, &Xloc);
520:   MPI_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da));
521:   return 0;
522: }

524: PetscErrorCode FVRHSFunctionslowbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
525: {
526:   FVCtx       *ctx = (FVCtx *)vctx;
527:   PetscInt     i, j, k, Mx, dof, xs, xm, islowbuffer = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
528:   PetscReal    hxs, hxm, hxf;
529:   PetscScalar *x, *f, *slope;
530:   Vec          Xloc;
531:   DM           da;

534:   TSGetDM(ts, &da);
535:   DMGetLocalVector(da, &Xloc);
536:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
537:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
538:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
539:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
540:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
541:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
542:   VecZeroEntries(F);
543:   DMDAVecGetArray(da, Xloc, &x);
544:   VecGetArray(F, &f);
545:   DMDAGetArray(da, PETSC_TRUE, &slope);
546:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

548:   if (ctx->bctype == FVBC_OUTFLOW) {
549:     for (i = xs - 2; i < 0; i++) {
550:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
551:     }
552:     for (i = Mx; i < xs + xm + 2; i++) {
553:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
554:     }
555:   }
556:   for (i = xs - 1; i < xs + xm + 1; i++) {
557:     struct _LimitInfo info;
558:     PetscScalar      *cjmpL, *cjmpR;
559:     if ((i > sm - lsbwidth - 2 && i < sm + 1) || (i > ms - 2 && i < ms + rsbwidth + 1)) {
560:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
561:       (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
562:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
563:       PetscArrayzero(ctx->cjmpLR, 2 * dof);
564:       cjmpL = &ctx->cjmpLR[0];
565:       cjmpR = &ctx->cjmpLR[dof];
566:       for (j = 0; j < dof; j++) {
567:         PetscScalar jmpL, jmpR;
568:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
569:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
570:         for (k = 0; k < dof; k++) {
571:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
572:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
573:         }
574:       }
575:       /* Apply limiter to the left and right characteristic jumps */
576:       info.m   = dof;
577:       info.hxs = hxs;
578:       info.hxm = hxm;
579:       info.hxf = hxf;
580:       (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
581:       for (j = 0; j < dof; j++) {
582:         PetscScalar tmp = 0;
583:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
584:         slope[i * dof + j] = tmp;
585:       }
586:     }
587:   }

589:   for (i = xs; i < xs + xm + 1; i++) {
590:     PetscReal    maxspeed;
591:     PetscScalar *uL, *uR;
592:     uL = &ctx->uLR[0];
593:     uR = &ctx->uLR[dof];
594:     if (i == sm - lsbwidth) {
595:       for (j = 0; j < dof; j++) {
596:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
597:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
598:       }
599:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
600:       if (i < xs + xm) {
601:         for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
602:         islowbuffer++;
603:       }
604:     }
605:     if (i > sm - lsbwidth && i < sm) {
606:       for (j = 0; j < dof; j++) {
607:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
608:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
609:       }
610:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
611:       if (i > xs) {
612:         for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
613:       }
614:       if (i < xs + xm) {
615:         for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
616:         islowbuffer++;
617:       }
618:     }
619:     if (i == sm) { /* interface between the slow region and the medium region */
620:       for (j = 0; j < dof; j++) {
621:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
622:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
623:       }
624:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
625:       if (i > xs) {
626:         for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
627:       }
628:     }
629:     if (i == ms) { /* interface between the medium region and the slow region */
630:       for (j = 0; j < dof; j++) {
631:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
632:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
633:       }
634:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
635:       if (i < xs + xm) {
636:         for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
637:         islowbuffer++;
638:       }
639:     }
640:     if (i > ms && i < ms + rsbwidth) {
641:       for (j = 0; j < dof; j++) {
642:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
643:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
644:       }
645:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
646:       if (i > xs) {
647:         for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
648:       }
649:       if (i < xs + xm) {
650:         for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
651:         islowbuffer++;
652:       }
653:     }
654:     if (i == ms + rsbwidth) {
655:       for (j = 0; j < dof; j++) {
656:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
657:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
658:       }
659:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
660:       if (i > xs) {
661:         for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
662:       }
663:     }
664:   }
665:   DMDAVecRestoreArray(da, Xloc, &x);
666:   VecRestoreArray(F, &f);
667:   DMDARestoreArray(da, PETSC_TRUE, &slope);
668:   DMRestoreLocalVector(da, &Xloc);
669:   return 0;
670: }

672: /* --------------------------------- Finite Volume Solver for medium components ----------------------------------- */
673: PetscErrorCode FVRHSFunctionmedium_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
674: {
675:   FVCtx       *ctx = (FVCtx *)vctx;
676:   PetscInt     i, j, k, Mx, dof, xs, xm, imedium = 0, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->ms, lmbwidth = ctx->lmbwidth, rmbwidth = ctx->rmbwidth;
677:   PetscReal    hxs, hxm, hxf;
678:   PetscScalar *x, *f, *slope;
679:   Vec          Xloc;
680:   DM           da;

683:   TSGetDM(ts, &da);
684:   DMGetLocalVector(da, &Xloc);
685:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
686:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
687:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
688:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
689:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
690:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
691:   VecZeroEntries(F);
692:   DMDAVecGetArray(da, Xloc, &x);
693:   VecGetArray(F, &f);
694:   DMDAGetArray(da, PETSC_TRUE, &slope);
695:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

697:   if (ctx->bctype == FVBC_OUTFLOW) {
698:     for (i = xs - 2; i < 0; i++) {
699:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
700:     }
701:     for (i = Mx; i < xs + xm + 2; i++) {
702:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
703:     }
704:   }
705:   for (i = xs - 1; i < xs + xm + 1; i++) {
706:     struct _LimitInfo info;
707:     PetscScalar      *cjmpL, *cjmpR;
708:     if ((i > sm - 2 && i < mf - lmbwidth + 1) || (i > fm + rmbwidth - 2 && i < ms + 1)) { /* slow components and the first and last fast components */
709:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
710:       (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
711:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
712:       PetscArrayzero(ctx->cjmpLR, 2 * dof);
713:       cjmpL = &ctx->cjmpLR[0];
714:       cjmpR = &ctx->cjmpLR[dof];
715:       for (j = 0; j < dof; j++) {
716:         PetscScalar jmpL, jmpR;
717:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
718:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
719:         for (k = 0; k < dof; k++) {
720:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
721:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
722:         }
723:       }
724:       /* Apply limiter to the left and right characteristic jumps */
725:       info.m   = dof;
726:       info.hxs = hxs;
727:       info.hxm = hxm;
728:       info.hxf = hxf;
729:       (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
730:       for (j = 0; j < dof; j++) {
731:         PetscScalar tmp = 0;
732:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
733:         slope[i * dof + j] = tmp;
734:       }
735:     }
736:   }

738:   for (i = xs; i < xs + xm + 1; i++) {
739:     PetscReal    maxspeed;
740:     PetscScalar *uL, *uR;
741:     uL = &ctx->uLR[0];
742:     uR = &ctx->uLR[dof];
743:     if (i == sm) { /* interface between slow and medium regions */
744:       for (j = 0; j < dof; j++) {
745:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
746:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
747:       }
748:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
749:       if (i < xs + xm) {
750:         for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
751:         imedium++;
752:       }
753:     }
754:     if (i > sm && i < mf - lmbwidth) { /* medium region */
755:       for (j = 0; j < dof; j++) {
756:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
757:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
758:       }
759:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
760:       if (i > xs) {
761:         for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
762:       }
763:       if (i < xs + xm) {
764:         for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
765:         imedium++;
766:       }
767:     }
768:     if (i == mf - lmbwidth) { /* interface between medium and fast regions */
769:       for (j = 0; j < dof; j++) {
770:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
771:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
772:       }
773:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
774:       if (i > xs) {
775:         for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
776:       }
777:     }
778:     if (i == fm + rmbwidth) { /* interface between fast and medium regions */
779:       for (j = 0; j < dof; j++) {
780:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
781:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
782:       }
783:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
784:       if (i < xs + xm) {
785:         for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
786:         imedium++;
787:       }
788:     }
789:     if (i > fm + rmbwidth && i < ms) { /* medium region */
790:       for (j = 0; j < dof; j++) {
791:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
792:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
793:       }
794:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
795:       if (i > xs) {
796:         for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
797:       }
798:       if (i < xs + xm) {
799:         for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
800:         imedium++;
801:       }
802:     }
803:     if (i == ms) { /* interface between medium and slow regions */
804:       for (j = 0; j < dof; j++) {
805:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
806:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
807:       }
808:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
809:       if (i > xs) {
810:         for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
811:       }
812:     }
813:   }
814:   DMDAVecRestoreArray(da, Xloc, &x);
815:   VecRestoreArray(F, &f);
816:   DMDARestoreArray(da, PETSC_TRUE, &slope);
817:   DMRestoreLocalVector(da, &Xloc);
818:   return 0;
819: }

821: PetscErrorCode FVRHSFunctionmediumbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
822: {
823:   FVCtx       *ctx = (FVCtx *)vctx;
824:   PetscInt     i, j, k, Mx, dof, xs, xm, imediumbuffer = 0, mf = ctx->mf, fm = ctx->fm, lmbwidth = ctx->lmbwidth, rmbwidth = ctx->rmbwidth;
825:   PetscReal    hxs, hxm, hxf;
826:   PetscScalar *x, *f, *slope;
827:   Vec          Xloc;
828:   DM           da;

831:   TSGetDM(ts, &da);
832:   DMGetLocalVector(da, &Xloc);
833:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
834:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
835:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
836:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
837:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
838:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
839:   VecZeroEntries(F);
840:   DMDAVecGetArray(da, Xloc, &x);
841:   VecGetArray(F, &f);
842:   DMDAGetArray(da, PETSC_TRUE, &slope);
843:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

845:   if (ctx->bctype == FVBC_OUTFLOW) {
846:     for (i = xs - 2; i < 0; i++) {
847:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
848:     }
849:     for (i = Mx; i < xs + xm + 2; i++) {
850:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
851:     }
852:   }
853:   for (i = xs - 1; i < xs + xm + 1; i++) {
854:     struct _LimitInfo info;
855:     PetscScalar      *cjmpL, *cjmpR;
856:     if ((i > mf - lmbwidth - 2 && i < mf + 1) || (i > fm - 2 && i < fm + rmbwidth + 1)) {
857:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
858:       (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
859:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
860:       PetscArrayzero(ctx->cjmpLR, 2 * dof);
861:       cjmpL = &ctx->cjmpLR[0];
862:       cjmpR = &ctx->cjmpLR[dof];
863:       for (j = 0; j < dof; j++) {
864:         PetscScalar jmpL, jmpR;
865:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
866:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
867:         for (k = 0; k < dof; k++) {
868:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
869:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
870:         }
871:       }
872:       /* Apply limiter to the left and right characteristic jumps */
873:       info.m   = dof;
874:       info.hxs = hxs;
875:       info.hxm = hxm;
876:       info.hxf = hxf;
877:       (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
878:       for (j = 0; j < dof; j++) {
879:         PetscScalar tmp = 0;
880:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
881:         slope[i * dof + j] = tmp;
882:       }
883:     }
884:   }

886:   for (i = xs; i < xs + xm + 1; i++) {
887:     PetscReal    maxspeed;
888:     PetscScalar *uL, *uR;
889:     uL = &ctx->uLR[0];
890:     uR = &ctx->uLR[dof];
891:     if (i == mf - lmbwidth) {
892:       for (j = 0; j < dof; j++) {
893:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
894:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
895:       }
896:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
897:       if (i < xs + xm) {
898:         for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
899:         imediumbuffer++;
900:       }
901:     }
902:     if (i > mf - lmbwidth && i < mf) {
903:       for (j = 0; j < dof; j++) {
904:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
905:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
906:       }
907:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
908:       if (i > xs) {
909:         for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
910:       }
911:       if (i < xs + xm) {
912:         for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
913:         imediumbuffer++;
914:       }
915:     }
916:     if (i == mf) { /* interface between the medium region and the fast region */
917:       for (j = 0; j < dof; j++) {
918:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
919:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
920:       }
921:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
922:       if (i > xs) {
923:         for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
924:       }
925:     }
926:     if (i == fm) { /* interface between the fast region and the medium region */
927:       for (j = 0; j < dof; j++) {
928:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
929:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
930:       }
931:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
932:       if (i < xs + xm) {
933:         for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
934:         imediumbuffer++;
935:       }
936:     }
937:     if (i > fm && i < fm + rmbwidth) {
938:       for (j = 0; j < dof; j++) {
939:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
940:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
941:       }
942:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
943:       if (i > xs) {
944:         for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
945:       }
946:       if (i < xs + xm) {
947:         for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
948:         imediumbuffer++;
949:       }
950:     }
951:     if (i == fm + rmbwidth) {
952:       for (j = 0; j < dof; j++) {
953:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
954:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
955:       }
956:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
957:       if (i > xs) {
958:         for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
959:       }
960:     }
961:   }
962:   DMDAVecRestoreArray(da, Xloc, &x);
963:   VecRestoreArray(F, &f);
964:   DMDARestoreArray(da, PETSC_TRUE, &slope);
965:   DMRestoreLocalVector(da, &Xloc);
966:   return 0;
967: }

969: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
970: PetscErrorCode FVRHSFunctionfast_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
971: {
972:   FVCtx       *ctx = (FVCtx *)vctx;
973:   PetscInt     i, j, k, Mx, dof, xs, xm, ifast = 0, mf = ctx->mf, fm = ctx->fm;
974:   PetscReal    hxs, hxm, hxf;
975:   PetscScalar *x, *f, *slope;
976:   Vec          Xloc;
977:   DM           da;

980:   TSGetDM(ts, &da);
981:   DMGetLocalVector(da, &Xloc);
982:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
983:   hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
984:   hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
985:   hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
986:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
987:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
988:   VecZeroEntries(F);
989:   DMDAVecGetArray(da, Xloc, &x);
990:   VecGetArray(F, &f);
991:   DMDAGetArray(da, PETSC_TRUE, &slope);
992:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

994:   if (ctx->bctype == FVBC_OUTFLOW) {
995:     for (i = xs - 2; i < 0; i++) {
996:       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
997:     }
998:     for (i = Mx; i < xs + xm + 2; i++) {
999:       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
1000:     }
1001:   }
1002:   for (i = xs - 1; i < xs + xm + 1; i++) { /* fast components and the last slow components before fast components and the first slow component after fast components */
1003:     struct _LimitInfo info;
1004:     PetscScalar      *cjmpL, *cjmpR;
1005:     if (i > mf - 2 && i < fm + 1) {
1006:       (*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds);
1007:       PetscArrayzero(ctx->cjmpLR, 2 * dof);
1008:       cjmpL = &ctx->cjmpLR[0];
1009:       cjmpR = &ctx->cjmpLR[dof];
1010:       for (j = 0; j < dof; j++) {
1011:         PetscScalar jmpL, jmpR;
1012:         jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
1013:         jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
1014:         for (k = 0; k < dof; k++) {
1015:           cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
1016:           cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
1017:         }
1018:       }
1019:       /* Apply limiter to the left and right characteristic jumps */
1020:       info.m   = dof;
1021:       info.hxs = hxs;
1022:       info.hxm = hxm;
1023:       info.hxf = hxf;
1024:       (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
1025:       for (j = 0; j < dof; j++) {
1026:         PetscScalar tmp = 0;
1027:         for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
1028:         slope[i * dof + j] = tmp;
1029:       }
1030:     }
1031:   }

1033:   for (i = xs; i < xs + xm + 1; i++) {
1034:     PetscReal    maxspeed;
1035:     PetscScalar *uL, *uR;
1036:     uL = &ctx->uLR[0];
1037:     uR = &ctx->uLR[dof];
1038:     if (i == mf) { /* interface between medium and fast regions */
1039:       for (j = 0; j < dof; j++) {
1040:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
1041:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1042:       }
1043:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
1044:       if (i < xs + xm) {
1045:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1046:         ifast++;
1047:       }
1048:     }
1049:     if (i > mf && i < fm) { /* fast region */
1050:       for (j = 0; j < dof; j++) {
1051:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1052:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1053:       }
1054:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
1055:       if (i > xs) {
1056:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1057:       }
1058:       if (i < xs + xm) {
1059:         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1060:         ifast++;
1061:       }
1062:     }
1063:     if (i == fm) { /* interface between fast and medium regions */
1064:       for (j = 0; j < dof; j++) {
1065:         uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1066:         uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
1067:       }
1068:       (*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed);
1069:       if (i > xs) {
1070:         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1071:       }
1072:     }
1073:   }
1074:   DMDAVecRestoreArray(da, Xloc, &x);
1075:   VecRestoreArray(F, &f);
1076:   DMDARestoreArray(da, PETSC_TRUE, &slope);
1077:   DMRestoreLocalVector(da, &Xloc);
1078:   return 0;
1079: }

1081: int main(int argc, char *argv[])
1082: {
1083:   char              lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
1084:   PetscFunctionList limiters = 0, physics = 0;
1085:   MPI_Comm          comm;
1086:   TS                ts;
1087:   DM                da;
1088:   Vec               X, X0, R;
1089:   FVCtx             ctx;
1090:   PetscInt          i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_medium, count_fast, islow = 0, imedium = 0, ifast = 0, *index_slow, *index_medium, *index_fast, islowbuffer = 0, *index_slowbuffer, imediumbuffer = 0, *index_mediumbuffer;
1091:   PetscBool         view_final = PETSC_FALSE;
1092:   PetscReal         ptime;

1095:   PetscInitialize(&argc, &argv, 0, help);
1096:   comm = PETSC_COMM_WORLD;
1097:   PetscMemzero(&ctx, sizeof(ctx));

1099:   /* Register limiters to be available on the command line */
1100:   PetscFunctionListAdd(&limiters, "upwind", Limit3_Upwind);
1101:   PetscFunctionListAdd(&limiters, "lax-wendroff", Limit3_LaxWendroff);
1102:   PetscFunctionListAdd(&limiters, "beam-warming", Limit3_BeamWarming);
1103:   PetscFunctionListAdd(&limiters, "fromm", Limit3_Fromm);
1104:   PetscFunctionListAdd(&limiters, "minmod", Limit3_Minmod);
1105:   PetscFunctionListAdd(&limiters, "superbee", Limit3_Superbee);
1106:   PetscFunctionListAdd(&limiters, "mc", Limit3_MC);
1107:   PetscFunctionListAdd(&limiters, "koren3", Limit3_Koren3);

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

1112:   ctx.comm   = comm;
1113:   ctx.cfl    = 0.9;
1114:   ctx.bctype = FVBC_PERIODIC;
1115:   ctx.xmin   = -1.0;
1116:   ctx.xmax   = 1.0;
1117:   PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
1118:   PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL);
1119:   PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL);
1120:   PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL);
1121:   PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL);
1122:   PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final);
1123:   PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL);
1124:   PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL);
1125:   PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL);
1126:   PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL);
1127:   PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL);
1128:   PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL);
1129:   PetscOptionsEnd();

1131:   /* Choose the limiter from the list of registered limiters */
1132:   PetscFunctionListFind(limiters, lname, &ctx.limit3);

1135:   /* Choose the physics from the list of registered models */
1136:   {
1137:     PetscErrorCode (*r)(FVCtx *);
1138:     PetscFunctionListFind(physics, physname, &r);
1140:     /* Create the physics, will set the number of fields and their names */
1141:     (*r)(&ctx);
1142:   }

1144:   /* Create a DMDA to manage the parallel grid */
1145:   DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da);
1146:   DMSetFromOptions(da);
1147:   DMSetUp(da);
1148:   /* Inform the DMDA of the field names provided by the physics. */
1149:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1150:   for (i = 0; i < ctx.physics2.dof; i++) DMDASetFieldName(da, i, ctx.physics2.fieldname[i]);
1151:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
1152:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

1154:   /* Set coordinates of cell centers */
1155:   DMDASetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax - ctx.xmin) / Mx, 0, 0, 0, 0);

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

1161:   /* Create a vector to store the solution and to save the initial state */
1162:   DMCreateGlobalVector(da, &X);
1163:   VecDuplicate(X, &X0);
1164:   VecDuplicate(X, &R);

1166:   /* create index for slow parts and fast parts,
1167:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
1168:   count_slow   = Mx / (1 + ctx.hratio) / (1 + ctx.hratio);
1169:   count_medium = 2 * ctx.hratio * count_slow;
1171:   count_fast = ctx.hratio * ctx.hratio * count_slow;
1172:   ctx.sm     = count_slow / 2;
1173:   ctx.mf     = ctx.sm + count_medium / 2;
1174:   ctx.fm     = ctx.mf + count_fast;
1175:   ctx.ms     = ctx.fm + count_medium / 2;
1176:   PetscMalloc1(xm * dof, &index_slow);
1177:   PetscMalloc1(xm * dof, &index_medium);
1178:   PetscMalloc1(xm * dof, &index_fast);
1179:   PetscMalloc1(6 * dof, &index_slowbuffer);
1180:   PetscMalloc1(6 * dof, &index_mediumbuffer);
1181:   if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
1182:     ctx.lsbwidth = 2;
1183:     ctx.rsbwidth = 4;
1184:     ctx.lmbwidth = 2;
1185:     ctx.rmbwidth = 4;
1186:   } else {
1187:     ctx.lsbwidth = 4;
1188:     ctx.rsbwidth = 2;
1189:     ctx.lmbwidth = 4;
1190:     ctx.rmbwidth = 2;
1191:   }

1193:   for (i = xs; i < xs + xm; i++) {
1194:     if (i < ctx.sm - ctx.lsbwidth || i > ctx.ms + ctx.rsbwidth - 1)
1195:       for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
1196:     else if ((i >= ctx.sm - ctx.lsbwidth && i < ctx.sm) || (i > ctx.ms - 1 && i <= ctx.ms + ctx.rsbwidth - 1))
1197:       for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
1198:     else if (i < ctx.mf - ctx.lmbwidth || i > ctx.fm + ctx.rmbwidth - 1)
1199:       for (k = 0; k < dof; k++) index_medium[imedium++] = i * dof + k;
1200:     else if ((i >= ctx.mf - ctx.lmbwidth && i < ctx.mf) || (i > ctx.fm - 1 && i <= ctx.fm + ctx.rmbwidth - 1))
1201:       for (k = 0; k < dof; k++) index_mediumbuffer[imediumbuffer++] = i * dof + k;
1202:     else
1203:       for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
1204:   }
1205:   ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss);
1206:   ISCreateGeneral(PETSC_COMM_WORLD, imedium, index_medium, PETSC_COPY_VALUES, &ctx.ism);
1207:   ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf);
1208:   ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb);
1209:   ISCreateGeneral(PETSC_COMM_WORLD, imediumbuffer, index_mediumbuffer, PETSC_COPY_VALUES, &ctx.ismb);

1211:   /* Create a time-stepping object */
1212:   TSCreate(comm, &ts);
1213:   TSSetDM(ts, da);
1214:   TSSetRHSFunction(ts, R, FVRHSFunction_3WaySplit, &ctx);
1215:   TSRHSSplitSetIS(ts, "slow", ctx.iss);
1216:   TSRHSSplitSetIS(ts, "medium", ctx.ism);
1217:   TSRHSSplitSetIS(ts, "fast", ctx.isf);
1218:   TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb);
1219:   TSRHSSplitSetIS(ts, "mediumbuffer", ctx.ismb);
1220:   TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_3WaySplit, &ctx);
1221:   TSRHSSplitSetRHSFunction(ts, "medium", NULL, FVRHSFunctionmedium_3WaySplit, &ctx);
1222:   TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_3WaySplit, &ctx);
1223:   TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_3WaySplit, &ctx);
1224:   TSRHSSplitSetRHSFunction(ts, "mediumbuffer", NULL, FVRHSFunctionmediumbuffer_3WaySplit, &ctx);

1226:   TSSetType(ts, TSSSP);
1227:   /*TSSetType(ts,TSMPRK);*/
1228:   TSSetMaxTime(ts, 10);
1229:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

1231:   /* Compute initial conditions and starting time step */
1232:   FVSample_3WaySplit(&ctx, da, 0, X0);
1233:   FVRHSFunction_3WaySplit(ts, 0, X0, X, (void *)&ctx); /* Initial function evaluation, only used to determine max speed */
1234:   VecCopy(X0, X);                                      /* The function value was not used so we set X=X0 again */
1235:   TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt);
1236:   TSSetFromOptions(ts); /* Take runtime options */
1237:   SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD);
1238:   {
1239:     PetscInt           steps;
1240:     PetscScalar        mass_initial, mass_final, mass_difference, mass_differenceg;
1241:     const PetscScalar *ptr_X, *ptr_X0;
1242:     const PetscReal    hs = (ctx.xmax - ctx.xmin) / 4.0 / count_slow;
1243:     const PetscReal    hm = (ctx.xmax - ctx.xmin) / 2.0 / count_medium;
1244:     const PetscReal    hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;

1246:     TSSolve(ts, X);
1247:     TSGetSolveTime(ts, &ptime);
1248:     TSGetStepNumber(ts, &steps);
1249:     /* calculate the total mass at initial time and final time */
1250:     mass_initial = 0.0;
1251:     mass_final   = 0.0;
1252:     DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0);
1253:     DMDAVecGetArrayRead(da, X, (void *)&ptr_X);
1254:     for (i = xs; i < xs + xm; i++) {
1255:       if (i < ctx.sm || i > ctx.ms - 1)
1256:         for (k = 0; k < dof; k++) {
1257:           mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
1258:           mass_final   = mass_final + hs * ptr_X[i * dof + k];
1259:         }
1260:       else if (i < ctx.mf || i > ctx.fm - 1)
1261:         for (k = 0; k < dof; k++) {
1262:           mass_initial = mass_initial + hm * ptr_X0[i * dof + k];
1263:           mass_final   = mass_final + hm * ptr_X[i * dof + k];
1264:         }
1265:       else {
1266:         for (k = 0; k < dof; k++) {
1267:           mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
1268:           mass_final   = mass_final + hf * ptr_X[i * dof + k];
1269:         }
1270:       }
1271:     }
1272:     DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0);
1273:     DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X);
1274:     mass_difference = mass_final - mass_initial;
1275:     MPI_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm);
1276:     PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg);
1277:     PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps);
1278:     PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1.0 / ctx.cfl_idt));
1279:     if (ctx.exact) {
1280:       PetscReal nrm1 = 0;
1281:       SolutionErrorNorms_3WaySplit(&ctx, da, ptime, X, &nrm1);
1282:       PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1);
1283:     }
1284:     if (ctx.simulation) {
1285:       PetscReal          nrm1 = 0;
1286:       PetscViewer        fd;
1287:       char               filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
1288:       Vec                XR;
1289:       PetscBool          flg;
1290:       const PetscScalar *ptr_XR;
1291:       PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg);
1293:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd);
1294:       VecDuplicate(X0, &XR);
1295:       VecLoad(XR, fd);
1296:       PetscViewerDestroy(&fd);
1297:       VecGetArrayRead(X, &ptr_X);
1298:       VecGetArrayRead(XR, &ptr_XR);
1299:       for (i = xs; i < xs + xm; i++) {
1300:         if (i < ctx.sm || i > ctx.ms - 1)
1301:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1302:         else if (i < ctx.mf || i < ctx.fm - 1)
1303:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hm * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1304:         else
1305:           for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1306:       }
1307:       VecRestoreArrayRead(X, &ptr_X);
1308:       VecRestoreArrayRead(XR, &ptr_XR);
1309:       PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1);
1310:       VecDestroy(&XR);
1311:     }
1312:   }

1314:   SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD);
1315:   if (draw & 0x1) VecView(X0, PETSC_VIEWER_DRAW_WORLD);
1316:   if (draw & 0x2) VecView(X, PETSC_VIEWER_DRAW_WORLD);
1317:   if (draw & 0x4) {
1318:     Vec Y;
1319:     VecDuplicate(X, &Y);
1320:     FVSample_3WaySplit(&ctx, da, ptime, Y);
1321:     VecAYPX(Y, -1, X);
1322:     VecView(Y, PETSC_VIEWER_DRAW_WORLD);
1323:     VecDestroy(&Y);
1324:   }

1326:   if (view_final) {
1327:     PetscViewer viewer;
1328:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer);
1329:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1330:     VecView(X, viewer);
1331:     PetscViewerPopFormat(viewer);
1332:     PetscViewerDestroy(&viewer);
1333:   }

1335:   /* Clean up */
1336:   (*ctx.physics2.destroy)(ctx.physics2.user);
1337:   for (i = 0; i < ctx.physics2.dof; i++) PetscFree(ctx.physics2.fieldname[i]);
1338:   PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope);
1339:   PetscFree3(ctx.uLR, ctx.flux, ctx.speeds);
1340:   VecDestroy(&X);
1341:   VecDestroy(&X0);
1342:   VecDestroy(&R);
1343:   DMDestroy(&da);
1344:   TSDestroy(&ts);
1345:   ISDestroy(&ctx.iss);
1346:   ISDestroy(&ctx.ism);
1347:   ISDestroy(&ctx.isf);
1348:   ISDestroy(&ctx.issb);
1349:   ISDestroy(&ctx.ismb);
1350:   PetscFree(index_slow);
1351:   PetscFree(index_medium);
1352:   PetscFree(index_fast);
1353:   PetscFree(index_slowbuffer);
1354:   PetscFree(index_mediumbuffer);
1355:   PetscFunctionListDestroy(&limiters);
1356:   PetscFunctionListDestroy(&physics);
1357:   PetscFinalize();
1358:   return 0;
1359: }

1361: /*TEST

1363:     build:
1364:       requires: !complex
1365:       depends: finitevolume1d.c

1367:     test:
1368:       suffix: 1
1369:       args: -da_grid_x 90 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a23 -ts_use_splitrhsfunction 0

1371:     test:
1372:       suffix: 2
1373:       args: -da_grid_x 90 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a23 -ts_use_splitrhsfunction 1

1375: TEST*/