2                            "  advection   - Constant coefficient scalar advection\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"
9 …"  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect…
11 …"  simulation  - use reference solution which is generated by smaller time step size to be true so…
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…
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))            …
19 …"                             r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1)))                  …
20 …"                             alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1))    …
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))…
31   PetscReal range = xmax - xmin;  in RangeMod()
35 /* --------------------------------- Finite Volume data structures --------------------------------…
71   PetscInt  sf, fs; /* slow-fast and fast-slow interfaces */
74 /* --------------------------------- Physics ----------------------------------- */
82 /* --------------------------------- Advection ----------------------------------- */
93   speed     = ctx->a;  in PhysicsFlux_Advect()
102   PetscReal  a   = ctx->a, x0;  in PhysicsSample_Advect()
107     x0 = x - a * t;  in PhysicsSample_Advect()
110     x0 = RangeMod(x - a * t, xmin, xmax);  in PhysicsSample_Advect()
117     u[0] = (x0 < 0) ? 1 : -1;  in PhysicsSample_Advect()
120     u[0] = (x0 < 0) ? -1 : 1;  in PhysicsSample_Advect()
135     u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));  in PhysicsSample_Advect()
152   ctx->physics.sample  = PhysicsSample_Advect;  in PhysicsCreate_Advect()
153   ctx->physics.flux    = PhysicsFlux_Advect;  in PhysicsCreate_Advect()
154   ctx->physics.destroy = PhysicsDestroy_SimpleFree;  in PhysicsCreate_Advect()
155   ctx->physics.user    = user;  in PhysicsCreate_Advect()
156   ctx->physics.dof     = 1;  in PhysicsCreate_Advect()
157   PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));  in PhysicsCreate_Advect()
158   user->a = 1;  in PhysicsCreate_Advect()
159   PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");  in PhysicsCreate_Advect()
161     PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));  in PhysicsCreate_Advect()
167 /* --------------------------------- Finite Volume Solver ----------------------------------- */
172   PetscInt     i, j, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;  in FVRHSFunction()  local
182   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;  in FVRHSFunction()
183   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;  in FVRHSFunction()
186   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */  in FVRHSFunction()
192   if (ctx->bctype == FVBC_OUTFLOW) {  in FVRHSFunction()
193     for (i = xs - 2; i < 0; i++) {  in FVRHSFunction()
194       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];  in FVRHSFunction()
197       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];  in FVRHSFunction()
205       u        = &ctx->u[0];  in FVRHSFunction()
208       for (j = 0; j < dof; j++) {  in FVRHSFunction()
209 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
210         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
211 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
213       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
216         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunction()
219         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;  in FVRHSFunction()
222       u        = &ctx->u[0];  in FVRHSFunction()
225       for (j = 0; j < dof; j++) {  in FVRHSFunction()
226 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
227         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
228 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
230       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
232         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunction()
235         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;  in FVRHSFunction()
238       u        = &ctx->u[0];  in FVRHSFunction()
241       for (j = 0; j < dof; j++) {  in FVRHSFunction()
242 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
243         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
244 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
246       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
248         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunction()
251         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;  in FVRHSFunction()
254       u        = &ctx->u[0];  in FVRHSFunction()
257       for (j = 0; j < dof; j++) {  in FVRHSFunction()
258 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
259         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
260 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
262       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
264         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunction()
267         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;  in FVRHSFunction()
270       u        = &ctx->u[0];  in FVRHSFunction()
273       for (j = 0; j < dof; j++) {  in FVRHSFunction()
274 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
275         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
276 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
278       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
280         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunction()
283         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;  in FVRHSFunction()
286       u        = &ctx->u[0];  in FVRHSFunction()
289       for (j = 0; j < dof; j++) {  in FVRHSFunction()
290 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunction()
291         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunction()
292 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunction()
294       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunction()
296         for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunction()
299         for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;  in FVRHSFunction()
306 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe…  in FVRHSFunction()
312 …ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\…  in FVRHSFunction()
321   PetscInt     i, j, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs;  in FVRHSFunctionslow()  local
331   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;  in FVRHSFunctionslow()
332   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;  in FVRHSFunctionslow()
335   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */  in FVRHSFunctionslow()
341   if (ctx->bctype == FVBC_OUTFLOW) {  in FVRHSFunctionslow()
342     for (i = xs - 2; i < 0; i++) {  in FVRHSFunctionslow()
343       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];  in FVRHSFunctionslow()
346       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];  in FVRHSFunctionslow()
354       u        = &ctx->u[0];  in FVRHSFunctionslow()
357       for (j = 0; j < dof; j++) {  in FVRHSFunctionslow()
358 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionslow()
359         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionslow()
360 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionslow()
362       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionslow()
364         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunctionslow()
367         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;  in FVRHSFunctionslow()
371       u        = &ctx->u[0];  in FVRHSFunctionslow()
374       for (j = 0; j < dof; j++) {  in FVRHSFunctionslow()
375 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionslow()
376         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionslow()
377 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionslow()
379       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionslow()
381         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunctionslow()
384       u        = &ctx->u[0];  in FVRHSFunctionslow()
387       for (j = 0; j < dof; j++) {  in FVRHSFunctionslow()
388 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionslow()
389         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionslow()
390 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionslow()
392       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionslow()
394         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;  in FVRHSFunctionslow()
398       u        = &ctx->u[0];  in FVRHSFunctionslow()
401       for (j = 0; j < dof; j++) {  in FVRHSFunctionslow()
402 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionslow()
403         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionslow()
404 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionslow()
406       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionslow()
408         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunctionslow()
411         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;  in FVRHSFunctionslow()
415       u        = &ctx->u[0];  in FVRHSFunctionslow()
418       for (j = 0; j < dof; j++) {  in FVRHSFunctionslow()
419 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionslow()
420         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionslow()
421 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionslow()
423       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionslow()
425         for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;  in FVRHSFunctionslow()
428         for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;  in FVRHSFunctionslow()
443   PetscInt     i, j, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;  in FVRHSFunctionfast()  local
453   hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;  in FVRHSFunctionfast()
454   hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;  in FVRHSFunctionfast()
457   PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */  in FVRHSFunctionfast()
463   if (ctx->bctype == FVBC_OUTFLOW) {  in FVRHSFunctionfast()
464     for (i = xs - 2; i < 0; i++) {  in FVRHSFunctionfast()
465       for (j = 0; j < dof; j++) x[i * dof + j] = x[j];  in FVRHSFunctionfast()
468       for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];  in FVRHSFunctionfast()
476       u        = &ctx->u[0];  in FVRHSFunctionfast()
479       for (j = 0; j < dof; j++) {  in FVRHSFunctionfast()
480 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionfast()
481         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionfast()
482 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionfast()
484       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionfast()
486         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;  in FVRHSFunctionfast()
490       u        = &ctx->u[0];  in FVRHSFunctionfast()
493       for (j = 0; j < dof; j++) {  in FVRHSFunctionfast()
494 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionfast()
495         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionfast()
496 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionfast()
498       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionfast()
500         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunctionfast()
503         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;  in FVRHSFunctionfast()
507       u        = &ctx->u[0];  in FVRHSFunctionfast()
510       for (j = 0; j < dof; j++) {  in FVRHSFunctionfast()
511 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionfast()
512         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionfast()
513 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionfast()
515       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionfast()
517         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunctionfast()
520         for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;  in FVRHSFunctionfast()
524       u        = &ctx->u[0];  in FVRHSFunctionfast()
527       for (j = 0; j < dof; j++) {  in FVRHSFunctionfast()
528 …  r[j]   = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);  in FVRHSFunctionfast()
529         min[j] = PetscMin(r[j], 2.0);  in FVRHSFunctionfast()
530 …u[j]   = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i …  in FVRHSFunctionfast()
532       PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));  in FVRHSFunctionfast()
534         for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;  in FVRHSFunctionfast()
545 /* --------------------------------- Finite Volume Solver for slow components ---------------------…
550   PetscInt       i, j, k, dof, xs, xm, Mx, count_slow, count_fast;  in FVSample()  local
554 …PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampli…  in FVSample()
559   const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;  in FVSample()
560   const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;  in FVSample()
561   count_slow         = Mx / (1 + ctx->hratio);  in FVSample()
562   count_fast         = Mx - count_slow;  in FVSample()
564     if (i * hs + 0.5 * hs < (ctx->xmax - ctx->xmin) * 0.25) {  in FVSample()
565       xi = ctx->xmin + 0.5 * hs + i * hs;  in FVSample()
568       for (j = 0; j < N + 1; j++) {  in FVSample()
569         xj = xi + hs * (j - N / 2) / (PetscReal)N;  in FVSample()
570 …PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xma…  in FVSample()
571         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;  in FVSample()
573 …} else if ((ctx->xmax - ctx->xmin) * 0.25 + (i - count_slow / 2) * hf + 0.5 * hf < (ctx->xmax - ct…  in FVSample()
574       xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.25 + 0.5 * hf + (i - count_slow / 2) * hf;  in FVSample()
577       for (j = 0; j < N + 1; j++) {  in FVSample()
578         xj = xi + hf * (j - N / 2) / (PetscReal)N;  in FVSample()
579 …PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xma…  in FVSample()
580         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;  in FVSample()
583 …xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.75 + 0.5 * hs + (i - count_slow / 2 - count_fast) * h…  in FVSample()
586       for (j = 0; j < N + 1; j++) {  in FVSample()
587         xj = xi + hs * (j - N / 2) / (PetscReal)N;  in FVSample()
588 …PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xma…  in FVSample()
589         for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;  in FVSample()
603   PetscInt           imin, imax, Mx, i, j, xs, xm, dof;  in SolutionStatsView()  local
619     for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);  in SolutionStatsView()
642   const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;  in SolutionErrorNorms()
643   const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;  in SolutionErrorNorms()
644   count_slow         = (PetscReal)Mx / (1.0 + ctx->hratio);  in SolutionErrorNorms()
645   count_fast         = Mx - count_slow;  in SolutionErrorNorms()
649 …f (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - pt…  in SolutionErrorNorms()
650     else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);  in SolutionErrorNorms()
682   ctx.xmin   = -1.0;  in main()
685   PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));  in main()
686   PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));  in main()
687 …PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final…  in main()
688 …PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given …  in main()
689 …PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initia…  in main()
690 …PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exa…  in main()
691 …PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simula…  in main()
692   PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));  in main()
693 …PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype,…  in main()
694   PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));  in main()
711   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */  in main()
717 …SetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax …  in main()
729 … PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) s…  in main()
730   count_fast = Mx - count_slow;  in main()
736     if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1)  in main()
744   /* Create a time-stepping object */  in main()
768     const PetscReal    hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow;  in main()
769     const PetscReal    hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast;  in main()
779       if (i < ctx.sf || i > ctx.fs - 1) {  in main()
793     mass_difference = mass_final - mass_initial;  in main()
800       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));  in main()
809       PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));  in main()
810 … PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");  in main()
818 …(i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] -…  in main()
819         else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]);  in main()
823       PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));  in main()
835     PetscCall(VecAYPX(Y, -1, X));  in main()
873 … -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type…
877 … -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type…
882 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts…
886 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts…
892 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts…