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"
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 */

87: static PetscErrorCode PhysicsFlux_Advect(void *vctx, const PetscScalar *u, PetscScalar *flux, PetscReal *maxspeed)
88: {
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: {
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: }

147: {

150:   PetscFunctionBeginUser;
151:   PetscCall(PetscNew(&user));
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:   PetscCall(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));
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:     PetscCall(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da)));
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;
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:   }
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 */

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));
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;
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:     }
794:     mass_difference = mass_final - mass_initial;
795:     PetscCall(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");
813:       PetscCall(VecDuplicate(X0, &XR));
815:       PetscCall(PetscViewerDestroy(&fd));
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:       }
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*/
```