Lines Matching +full:- +full:j
3 -hratio is the ratio between mesh size of coarse grids and fine grids
6 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscret…
7 " advect - Constant coefficient scalar advection\n"
9 " shallow - 1D Shallow water equations (Saint Venant System)\n"
13 … this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
16 …" exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect…
18 …" simulation - use reference solution which is generated by smaller time step size to be true so…
20 … " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
21 … " bc_type - Boundary condition for the problem, options are: periodic, outflow, inflow "
22 …iptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n"
23 "The problem size should be set with -da_grid_x M\n\n";
27 …-da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 7.0 -ts_type mprk -ts…
28 …-da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 2.5 -ts_type mprk -ts…
29 …-da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 4.0 -ts_type mprk -ts…
30 …-da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_time_step 0.01 -ts_max_time 4.0 -ts_type mprk…
31 …-da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 5.0 -ts_type mprk -ts…
45 PetscReal range = xmax - xmin; in RangeMod()
52 /* --------------------------------- Advection ----------------------------------- */
63 speed = ctx->a; in PhysicsRiemann_Advect()
76 speeds[0] = ctx->a; in PhysicsCharacteristic_Advect()
83 PetscReal a = ctx->a, x0; in PhysicsSample_Advect()
88 x0 = x - a * t; in PhysicsSample_Advect()
91 x0 = RangeMod(x - a * t, xmin, xmax); in PhysicsSample_Advect()
98 u[0] = (x0 < 0) ? 1 : -1; in PhysicsSample_Advect()
101 u[0] = (x0 < 0) ? -1 : 1; in PhysicsSample_Advect()
116 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); in PhysicsSample_Advect()
133 ctx->physics2.sample2 = PhysicsSample_Advect; in PhysicsCreate_Advect()
134 ctx->physics2.riemann2 = PhysicsRiemann_Advect; in PhysicsCreate_Advect()
135 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; in PhysicsCreate_Advect()
136 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; in PhysicsCreate_Advect()
137 ctx->physics2.user = user; in PhysicsCreate_Advect()
138 ctx->physics2.dof = 1; in PhysicsCreate_Advect()
140 PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0])); in PhysicsCreate_Advect()
141 user->a = 1; in PhysicsCreate_Advect()
142 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); in PhysicsCreate_Advect()
143 PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL)); in PhysicsCreate_Advect()
147 /* --------------------------------- Shallow Water ----------------------------------- */
156 f[1] = PetscSqr(u[1]) / u[0] + 0.5 * phys->gravity * PetscSqr(u[0]); in ShallowFlux()
162 PetscScalar g = phys->gravity, ustar[2], cL, cR, c, cstar; in PhysicsRiemann_Shallow_Exact()
180 … fl = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */ in PhysicsRiemann_Shallow_Exact()
181 … : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h); /* rarefaction */ in PhysicsRiemann_Shallow_Exact()
182 … fr = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */ in PhysicsRiemann_Shallow_Exact()
183 … : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h); /* rarefaction */ in PhysicsRiemann_Shallow_Exact()
184 res = R.u - L.u + fr + fl; in PhysicsRiemann_Shallow_Exact()
185 …PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number g… in PhysicsRiemann_Shallow_Exact()
186 …if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h - h0) < PETSC_S… in PhysicsRiemann_Shallow_Exact()
188 star.u = L.u - fl; in PhysicsRiemann_Shallow_Exact()
197 …dfl = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar… in PhysicsRiemann_Shallow_Exact()
198 …dfr = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar… in PhysicsRiemann_Shallow_Exact()
199 tmp = h - res / (dfr + dfl); in PhysicsRiemann_Shallow_Exact()
202 …PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h… in PhysicsRiemann_Shallow_Exact()
208 if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */ in PhysicsRiemann_Shallow_Exact()
213 } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */ in PhysicsRiemann_Shallow_Exact()
215 ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR); in PhysicsRiemann_Shallow_Exact()
216 ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0]; in PhysicsRiemann_Shallow_Exact()
218 …} else if ((L.h >= star.h && L.u - c >= 0) || (L.h < star.h && (star.h * star.u - L.h * L.u) / (st… in PhysicsRiemann_Shallow_Exact()
219 /* 1-wave is right-travelling shock (supersonic) */ in PhysicsRiemann_Shallow_Exact()
221 …(star.h <= R.h && R.u + c <= 0) || (star.h > R.h && (R.h * R.u - star.h * star.h) / (R.h - star.h)… in PhysicsRiemann_Shallow_Exact()
222 /* 2-wave is left-travelling shock (supersonic) */ in PhysicsRiemann_Shallow_Exact()
229 *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR)); in PhysicsRiemann_Shallow_Exact()
236 PetscScalar g = phys->gravity, fL[2], fR[2], s; in PhysicsRiemann_Shallow_Rusanov()
240 PetscReal tol = 1e-6; in PhysicsRiemann_Shallow_Rusanov()
255 flux[0] = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (L.h - R.h); in PhysicsRiemann_Shallow_Rusanov()
256 flux[1] = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]); in PhysicsRiemann_Shallow_Rusanov()
263 PetscInt i, j; in PhysicsCharacteristic_Conservative() local
267 for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j); in PhysicsCharacteristic_Conservative()
277 PetscReal tol = 1e-6; in PhysicsCharacteristic_Shallow()
280 c = PetscSqrtScalar(u[0] * phys->gravity); in PhysicsCharacteristic_Shallow()
288 speeds[0] = u[1] / u[0] - c; in PhysicsCharacteristic_Shallow()
316 u[1] = (x < 25) ? -5 : 5; in PhysicsSample_Shallow()
319 u[0] = (x < 20) ? 1 : 1e-12; in PhysicsSample_Shallow()
323 u[0] = (x < 30) ? 1e-12 : 1; in PhysicsSample_Shallow()
328 u[1] = (x < 25) ? -0.3 : 0.3; in PhysicsSample_Shallow()
335 if (x < -0.1) { in PhysicsSample_Shallow()
336 u[0] = 1e-9; in PhysicsSample_Shallow()
342 u[0] = 1e-9; in PhysicsSample_Shallow()
347 if (x < -0.1) { in PhysicsSample_Shallow()
371 if (ctx->bctype == FVBC_INFLOW) { in PhysicsInflow_Shallow()
372 switch (ctx->initial) { in PhysicsInflow_Shallow()
400 u[3] = -1.0; /* Right boundary conditions */ in PhysicsInflow_Shallow()
413 switch (ctx->initial) { in PhysicsSetInflowType_Shallow()
423 ctx->physics2.bcinflowindex[0] = PETSC_FALSE; in PhysicsSetInflowType_Shallow()
424 ctx->physics2.bcinflowindex[1] = PETSC_TRUE; in PhysicsSetInflowType_Shallow()
425 ctx->physics2.bcinflowindex[2] = PETSC_FALSE; in PhysicsSetInflowType_Shallow()
426 ctx->physics2.bcinflowindex[3] = PETSC_TRUE; in PhysicsSetInflowType_Shallow()
442 ctx->physics2.sample2 = PhysicsSample_Shallow; in PhysicsCreate_Shallow()
443 ctx->physics2.inflow = PhysicsInflow_Shallow; in PhysicsCreate_Shallow()
444 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; in PhysicsCreate_Shallow()
445 ctx->physics2.riemann2 = PhysicsRiemann_Shallow_Rusanov; in PhysicsCreate_Shallow()
446 ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow; in PhysicsCreate_Shallow()
447 ctx->physics2.user = user; in PhysicsCreate_Shallow()
448 ctx->physics2.dof = 2; in PhysicsCreate_Shallow()
450 PetscCall(PetscMalloc1(2 * (ctx->physics2.dof), &ctx->physics2.bcinflowindex)); in PhysicsCreate_Shallow()
453 PetscCall(PetscStrallocpy("height", &ctx->physics2.fieldname[0])); in PhysicsCreate_Shallow()
454 PetscCall(PetscStrallocpy("momentum", &ctx->physics2.fieldname[1])); in PhysicsCreate_Shallow()
456 user->gravity = 9.81; in PhysicsCreate_Shallow()
462 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", ""); in PhysicsCreate_Shallow()
463 …PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravit… in PhysicsCreate_Shallow()
464 …PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname,… in PhysicsCreate_Shallow()
465 …PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, … in PhysicsCreate_Shallow()
467 PetscCall(RiemannListFind_2WaySplit(rlist, rname, &ctx->physics2.riemann2)); in PhysicsCreate_Shallow()
468 PetscCall(ReconstructListFind_2WaySplit(rclist, rcname, &ctx->physics2.characteristic2)); in PhysicsCreate_Shallow()
477 PetscInt i, j, k, dof, xs, xm, Mx; in FVSample_2WaySplit() local
482 …PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a samp… in FVSample_2WaySplit()
487 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVSample_2WaySplit()
488 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVSample_2WaySplit()
490 if (i < ctx->sf) { in FVSample_2WaySplit()
491 xi = ctx->xmin + 0.5 * hs + i * hs; in FVSample_2WaySplit()
494 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
495 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
496 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
497 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
499 } else if (i < ctx->fs) { in FVSample_2WaySplit()
500 xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf; in FVSample_2WaySplit()
503 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
504 xj = xi + hf * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
505 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
506 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
509 xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs; in FVSample_2WaySplit()
512 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
513 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
514 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
515 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
535 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in SolutionErrorNorms_2WaySplit()
536 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in SolutionErrorNorms_2WaySplit()
540 if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_2WaySplit()
541 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_2WaySplit()
552 PetscInt i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs; in FVRHSFunction_2WaySplit() local
562 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunction_2WaySplit()
563 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunction_2WaySplit()
567 PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */ in FVRHSFunction_2WaySplit()
574 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunction_2WaySplit()
575 for (i = xs - 2; i < 0; i++) { in FVRHSFunction_2WaySplit()
576 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunction_2WaySplit()
579 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunction_2WaySplit()
583 if (ctx->bctype == FVBC_INFLOW) { in FVRHSFunction_2WaySplit()
585 pages 137-138 for the scheme. */ in FVRHSFunction_2WaySplit()
587 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); in FVRHSFunction_2WaySplit()
588 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
589 if (ctx->physics2.bcinflowindex[j]) { in FVRHSFunction_2WaySplit()
590 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; in FVRHSFunction_2WaySplit()
592 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ in FVRHSFunction_2WaySplit()
597 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); in FVRHSFunction_2WaySplit()
598 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
599 if (ctx->physics2.bcinflowindex[dof + j]) { in FVRHSFunction_2WaySplit()
600 …for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof… in FVRHSFunction_2WaySplit()
602 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ in FVRHSFunction_2WaySplit()
608 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunction_2WaySplit()
611 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunction_2WaySplit()
612 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunction_2WaySplit()
613 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunction_2WaySplit()
614 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunction_2WaySplit()
615 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunction_2WaySplit()
616 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunction_2WaySplit()
617 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
619 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunction_2WaySplit()
620 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunction_2WaySplit()
622 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunction_2WaySplit()
623 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunction_2WaySplit()
630 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunction_2WaySplit()
631 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
633 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunction_2WaySplit()
634 slope[i * dof + j] = tmp; in FVRHSFunction_2WaySplit()
641 uL = &ctx->uLR[0]; in FVRHSFunction_2WaySplit()
642 uR = &ctx->uLR[dof]; in FVRHSFunction_2WaySplit()
644 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
645 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
646 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
648 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
650 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
653 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
656 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
657 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
658 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
660 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
662 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
665 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
668 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
669 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
670 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
672 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
674 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
677 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
680 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
681 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
682 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
684 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
686 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
689 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
692 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
693 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
694 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
696 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
699 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
702 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
710 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunction_2WaySplit()
716 if (dt > 0.5 / ctx->cfl_idt) { in FVRHSFunction_2WaySplit()
717 …scCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow,… in FVRHSFunction_2WaySplit()
718 …ERR_PLIB, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt)); in FVRHSFunction_2WaySplit()
724 /* --------------------------------- Finite Volume Solver for slow components ---------------------…
728 …PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbw… in FVRHSFunctionslow_2WaySplit() local
738 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionslow_2WaySplit()
739 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionslow_2WaySplit()
748 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslow_2WaySplit()
749 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslow_2WaySplit()
750 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslow_2WaySplit()
753 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslow_2WaySplit()
756 if (ctx->bctype == FVBC_INFLOW) { in FVRHSFunctionslow_2WaySplit()
758 pages 137-138 for the scheme. */ in FVRHSFunctionslow_2WaySplit()
760 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); in FVRHSFunctionslow_2WaySplit()
761 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
762 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { in FVRHSFunctionslow_2WaySplit()
763 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; in FVRHSFunctionslow_2WaySplit()
765 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ in FVRHSFunctionslow_2WaySplit()
770 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); in FVRHSFunctionslow_2WaySplit()
771 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
772 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { in FVRHSFunctionslow_2WaySplit()
773 …for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof… in FVRHSFunctionslow_2WaySplit()
775 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ in FVRHSFunctionslow_2WaySplit()
780 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslow_2WaySplit()
783 …if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fa… in FVRHSFunctionslow_2WaySplit()
784 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslow_2WaySplit()
785 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslow_2WaySplit()
786 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslow_2WaySplit()
787 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslow_2WaySplit()
788 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslow_2WaySplit()
789 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslow_2WaySplit()
790 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
792 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslow_2WaySplit()
793 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslow_2WaySplit()
795 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslow_2WaySplit()
796 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslow_2WaySplit()
803 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionslow_2WaySplit()
804 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
806 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslow_2WaySplit()
807 slope[i * dof + j] = tmp; in FVRHSFunctionslow_2WaySplit()
815 uL = &ctx->uLR[0]; in FVRHSFunctionslow_2WaySplit()
816 uR = &ctx->uLR[dof]; in FVRHSFunctionslow_2WaySplit()
817 if (i < sf - lsbwidth) { /* slow region */ in FVRHSFunctionslow_2WaySplit()
818 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
819 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
820 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
822 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
825 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
828 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
832 if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */ in FVRHSFunctionslow_2WaySplit()
833 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
834 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
835 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
837 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
839 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
843 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
844 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
845 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
847 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
849 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
854 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
855 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
856 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
858 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
860 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
863 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
872 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunctionslow_2WaySplit()
879 …PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbw… in FVRHSFunctionslowbuffer_2WaySplit() local
889 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionslowbuffer_2WaySplit()
890 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionslowbuffer_2WaySplit()
899 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslowbuffer_2WaySplit()
900 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslowbuffer_2WaySplit()
901 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslowbuffer_2WaySplit()
904 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
907 if (ctx->bctype == FVBC_INFLOW) { in FVRHSFunctionslowbuffer_2WaySplit()
909 pages 137-138 for the scheme. */ in FVRHSFunctionslowbuffer_2WaySplit()
911 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); in FVRHSFunctionslowbuffer_2WaySplit()
912 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
913 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { in FVRHSFunctionslowbuffer_2WaySplit()
914 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
916 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ in FVRHSFunctionslowbuffer_2WaySplit()
921 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); in FVRHSFunctionslowbuffer_2WaySplit()
922 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
923 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { in FVRHSFunctionslowbuffer_2WaySplit()
924 …for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof… in FVRHSFunctionslowbuffer_2WaySplit()
926 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ in FVRHSFunctionslowbuffer_2WaySplit()
931 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslowbuffer_2WaySplit()
934 if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) { in FVRHSFunctionslowbuffer_2WaySplit()
935 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslowbuffer_2WaySplit()
936 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslowbuffer_2WaySplit()
937 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslowbuffer_2WaySplit()
938 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslowbuffer_2WaySplit()
939 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslowbuffer_2WaySplit()
940 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslowbuffer_2WaySplit()
941 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
943 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
944 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
946 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslowbuffer_2WaySplit()
947 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslowbuffer_2WaySplit()
954 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionslowbuffer_2WaySplit()
955 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
957 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslowbuffer_2WaySplit()
958 slope[i * dof + j] = tmp; in FVRHSFunctionslowbuffer_2WaySplit()
966 uL = &ctx->uLR[0]; in FVRHSFunctionslowbuffer_2WaySplit()
967 uR = &ctx->uLR[dof]; in FVRHSFunctionslowbuffer_2WaySplit()
968 if (i == sf - lsbwidth) { in FVRHSFunctionslowbuffer_2WaySplit()
969 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
970 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
971 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
973 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
975 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
979 if (i > sf - lsbwidth && i < sf) { in FVRHSFunctionslowbuffer_2WaySplit()
980 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
981 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
982 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
984 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
986 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
989 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
994 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
995 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
996 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionslowbuffer_2WaySplit()
998 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
1000 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
1004 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
1005 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1006 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1008 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
1010 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
1015 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
1016 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1017 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1019 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
1021 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
1024 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
1029 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
1030 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1031 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
1033 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
1035 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
1046 /* --------------------------------- Finite Volume Solver for fast parts -------------------------…
1050 PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs; in FVRHSFunctionfast_2WaySplit() local
1060 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionfast_2WaySplit()
1061 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionfast_2WaySplit()
1070 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionfast_2WaySplit()
1071 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionfast_2WaySplit()
1072 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionfast_2WaySplit()
1075 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionfast_2WaySplit()
1078 if (ctx->bctype == FVBC_INFLOW) { in FVRHSFunctionfast_2WaySplit()
1080 pages 137-138 for the scheme. */ in FVRHSFunctionfast_2WaySplit()
1082 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); in FVRHSFunctionfast_2WaySplit()
1083 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1084 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { in FVRHSFunctionfast_2WaySplit()
1085 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; in FVRHSFunctionfast_2WaySplit()
1087 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ in FVRHSFunctionfast_2WaySplit()
1092 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); in FVRHSFunctionfast_2WaySplit()
1093 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1094 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { in FVRHSFunctionfast_2WaySplit()
1095 …for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof… in FVRHSFunctionfast_2WaySplit()
1097 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ in FVRHSFunctionfast_2WaySplit()
1102 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionfast_2WaySplit()
1105 if (i > sf - 2 && i < fs + 1) { in FVRHSFunctionfast_2WaySplit()
1106 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionfast_2WaySplit()
1107 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionfast_2WaySplit()
1108 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionfast_2WaySplit()
1109 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionfast_2WaySplit()
1110 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1112 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionfast_2WaySplit()
1113 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionfast_2WaySplit()
1115 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionfast_2WaySplit()
1116 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionfast_2WaySplit()
1123 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionfast_2WaySplit()
1124 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1126 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionfast_2WaySplit()
1127 slope[i * dof + j] = tmp; in FVRHSFunctionfast_2WaySplit()
1135 uL = &ctx->uLR[0]; in FVRHSFunctionfast_2WaySplit()
1136 uR = &ctx->uLR[dof]; in FVRHSFunctionfast_2WaySplit()
1138 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1139 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionfast_2WaySplit()
1140 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
1142 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
1144 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
1149 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1150 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
1151 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
1153 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
1155 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
1158 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
1163 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
1164 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
1165 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionfast_2WaySplit()
1167 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
1169 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
1200 PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff)); in main()
1201 PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming)); in main()
1215 ctx.xmin = -1.0; in main()
1222 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); in main()
1223 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); in main()
1224 …PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, si… in main()
1225 …PetscCall(PetscOptionsFList("-physics", "Name of physics model to use", "", physics, physname, phy… in main()
1226 …PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final… in main()
1227 …PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given … in main()
1228 …PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initia… in main()
1229 …PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exa… in main()
1230 …PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simula… in main()
1231 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); in main()
1232 …PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype,… in main()
1233 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL)); in main()
1254 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ in main()
1260 …SetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax … in main()
1275 … PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) s… in main()
1276 count_fast = Mx - count_slow; in main()
1286 if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1) in main()
1288 …else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwid… in main()
1297 /* Create a time-stepping object */ in main()
1323 const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow; in main()
1324 const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast; in main()
1335 if (i < ctx.sf || i > ctx.fs - 1) { in main()
1349 mass_difference = mass_final - mass_initial; in main()
1357 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
1366 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); in main()
1367 … PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); in main()
1375 if (i < ctx.sf || i > ctx.fs - 1) in main()
1376 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
1378 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
1382 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
1394 PetscCall(VecAYPX(Y, -1, X)); in main()
1440 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…
1445 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…
1450 …-da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -…
1456 …-da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -…
1462 …-da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -…