Lines Matching +full:- +full:j
3 -hratio is the ratio between mesh size of coarse grids and fine grids
4 -ts_rk_dtratio is the ratio between the large stepsize and the small stepsize
7 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscret…
8 " advection - Constant coefficient scalar advection\n"
10 … this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
13 …" exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect…
15 …" simulation - use reference solution which is generated by smaller time step size to be true so…
17 … " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
18 "Several initial conditions can be chosen with -initial N\n\n"
19 "The problem size should be set with -da_grid_x M\n\n";
29 PetscReal range = xmax - xmin; in RangeMod()
33 /* --------------------------------- Advection ----------------------------------- */
44 speed = ctx->a; in PhysicsRiemann_Advect()
57 speeds[0] = ctx->a; in PhysicsCharacteristic_Advect()
64 PetscReal a = ctx->a, x0; in PhysicsSample_Advect()
69 x0 = x - a * t; in PhysicsSample_Advect()
72 x0 = RangeMod(x - a * t, xmin, xmax); in PhysicsSample_Advect()
79 u[0] = (x0 < 0) ? 1 : -1; in PhysicsSample_Advect()
82 u[0] = (x0 < 0) ? -1 : 1; in PhysicsSample_Advect()
97 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); in PhysicsSample_Advect()
114 ctx->physics2.sample2 = PhysicsSample_Advect; in PhysicsCreate_Advect()
115 ctx->physics2.riemann2 = PhysicsRiemann_Advect; in PhysicsCreate_Advect()
116 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; in PhysicsCreate_Advect()
117 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; in PhysicsCreate_Advect()
118 ctx->physics2.user = user; in PhysicsCreate_Advect()
119 ctx->physics2.dof = 1; in PhysicsCreate_Advect()
120 PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0])); in PhysicsCreate_Advect()
121 user->a = 1; in PhysicsCreate_Advect()
122 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); in PhysicsCreate_Advect()
124 PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL)); in PhysicsCreate_Advect()
133 PetscInt i, j, k, dof, xs, xm, Mx; in FVSample_2WaySplit() local
138 …PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a samp… in FVSample_2WaySplit()
143 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVSample_2WaySplit()
144 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVSample_2WaySplit()
146 if (i < ctx->sf) { in FVSample_2WaySplit()
147 xi = ctx->xmin + 0.5 * hs + i * hs; in FVSample_2WaySplit()
150 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
151 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
152 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
153 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
155 } else if (i < ctx->fs) { in FVSample_2WaySplit()
156 xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf; in FVSample_2WaySplit()
159 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
160 xj = xi + hf * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
161 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
162 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
165 xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs; in FVSample_2WaySplit()
168 for (j = 0; j < N + 1; j++) { in FVSample_2WaySplit()
169 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_2WaySplit()
170 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_2WaySplit()
171 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_2WaySplit()
191 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in SolutionErrorNorms_2WaySplit()
192 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in SolutionErrorNorms_2WaySplit()
196 if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_2WaySplit()
197 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_2WaySplit()
208 PetscInt i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs; in FVRHSFunction_2WaySplit() local
218 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunction_2WaySplit()
219 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunction_2WaySplit()
223 PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */ in FVRHSFunction_2WaySplit()
231 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunction_2WaySplit()
232 for (i = xs - 2; i < 0; i++) { in FVRHSFunction_2WaySplit()
233 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunction_2WaySplit()
236 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunction_2WaySplit()
239 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunction_2WaySplit()
242 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunction_2WaySplit()
243 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunction_2WaySplit()
244 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunction_2WaySplit()
245 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunction_2WaySplit()
246 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunction_2WaySplit()
247 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunction_2WaySplit()
248 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
250 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunction_2WaySplit()
251 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunction_2WaySplit()
253 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunction_2WaySplit()
254 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunction_2WaySplit()
261 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunction_2WaySplit()
262 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
264 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunction_2WaySplit()
265 slope[i * dof + j] = tmp; in FVRHSFunction_2WaySplit()
272 uL = &ctx->uLR[0]; in FVRHSFunction_2WaySplit()
273 uR = &ctx->uLR[dof]; in FVRHSFunction_2WaySplit()
275 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
276 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
277 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
279 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
281 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
284 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
287 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
288 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
289 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
291 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
293 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
296 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
299 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
300 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
301 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
303 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
305 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
308 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
311 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
312 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_2WaySplit()
313 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
315 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
317 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_2WaySplit()
320 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
323 for (j = 0; j < dof; j++) { in FVRHSFunction_2WaySplit()
324 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
325 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_2WaySplit()
327 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_2WaySplit()
330 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
333 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_2WaySplit()
341 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunction_2WaySplit()
347 if (dt > 0.5 / ctx->cfl_idt) { in FVRHSFunction_2WaySplit()
348 …cCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, … in FVRHSFunction_2WaySplit()
349 …OWN_TYPE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt)); in FVRHSFunction_2WaySplit()
355 /* --------------------------------- Finite Volume Solver for slow components ---------------------…
359 …PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbw… in FVRHSFunctionslow_2WaySplit() local
369 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionslow_2WaySplit()
370 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionslow_2WaySplit()
379 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslow_2WaySplit()
380 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslow_2WaySplit()
381 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslow_2WaySplit()
384 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslow_2WaySplit()
387 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslow_2WaySplit()
390 …if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fa… in FVRHSFunctionslow_2WaySplit()
391 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslow_2WaySplit()
392 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslow_2WaySplit()
393 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslow_2WaySplit()
394 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslow_2WaySplit()
395 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslow_2WaySplit()
396 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslow_2WaySplit()
397 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
399 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslow_2WaySplit()
400 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslow_2WaySplit()
402 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslow_2WaySplit()
403 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslow_2WaySplit()
410 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionslow_2WaySplit()
411 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
413 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslow_2WaySplit()
414 slope[i * dof + j] = tmp; in FVRHSFunctionslow_2WaySplit()
422 uL = &ctx->uLR[0]; in FVRHSFunctionslow_2WaySplit()
423 uR = &ctx->uLR[dof]; in FVRHSFunctionslow_2WaySplit()
424 if (i < sf - lsbwidth) { /* slow region */ in FVRHSFunctionslow_2WaySplit()
425 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
426 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
427 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
429 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
432 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
435 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
439 if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */ in FVRHSFunctionslow_2WaySplit()
440 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
441 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
442 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
444 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
446 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
450 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
451 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
452 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
454 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
456 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
461 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_2WaySplit()
462 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
463 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_2WaySplit()
465 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_2WaySplit()
467 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
470 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_2WaySplit()
479 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunctionslow_2WaySplit()
486 …PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbw… in FVRHSFunctionslowbuffer_2WaySplit() local
496 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionslowbuffer_2WaySplit()
497 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionslowbuffer_2WaySplit()
506 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslowbuffer_2WaySplit()
507 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslowbuffer_2WaySplit()
508 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslowbuffer_2WaySplit()
511 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
514 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslowbuffer_2WaySplit()
517 if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) { in FVRHSFunctionslowbuffer_2WaySplit()
518 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslowbuffer_2WaySplit()
519 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslowbuffer_2WaySplit()
520 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslowbuffer_2WaySplit()
521 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslowbuffer_2WaySplit()
522 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslowbuffer_2WaySplit()
523 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslowbuffer_2WaySplit()
524 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
526 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
527 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslowbuffer_2WaySplit()
529 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslowbuffer_2WaySplit()
530 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslowbuffer_2WaySplit()
537 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionslowbuffer_2WaySplit()
538 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
540 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslowbuffer_2WaySplit()
541 slope[i * dof + j] = tmp; in FVRHSFunctionslowbuffer_2WaySplit()
549 uL = &ctx->uLR[0]; in FVRHSFunctionslowbuffer_2WaySplit()
550 uR = &ctx->uLR[dof]; in FVRHSFunctionslowbuffer_2WaySplit()
551 if (i == sf - lsbwidth) { in FVRHSFunctionslowbuffer_2WaySplit()
552 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
553 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
554 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
556 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
558 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
562 if (i > sf - lsbwidth && i < sf) { in FVRHSFunctionslowbuffer_2WaySplit()
563 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
564 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
565 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
567 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
569 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
572 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
577 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
578 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
579 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionslowbuffer_2WaySplit()
581 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
583 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
587 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
588 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionslowbuffer_2WaySplit()
589 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
591 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
593 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
598 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
599 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
600 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
602 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
604 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
607 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
612 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_2WaySplit()
613 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
614 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_2WaySplit()
616 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_2WaySplit()
618 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_2WaySplit()
629 /* --------------------------------- Finite Volume Solver for fast parts -------------------------…
633 PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs; in FVRHSFunctionfast_2WaySplit() local
643 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; in FVRHSFunctionfast_2WaySplit()
644 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); in FVRHSFunctionfast_2WaySplit()
653 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionfast_2WaySplit()
654 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionfast_2WaySplit()
655 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionfast_2WaySplit()
658 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionfast_2WaySplit()
661 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionfast_2WaySplit()
664 if (i > sf - 2 && i < fs + 1) { in FVRHSFunctionfast_2WaySplit()
665 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionfast_2WaySplit()
666 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionfast_2WaySplit()
667 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionfast_2WaySplit()
668 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionfast_2WaySplit()
669 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
671 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionfast_2WaySplit()
672 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionfast_2WaySplit()
674 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionfast_2WaySplit()
675 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionfast_2WaySplit()
682 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); in FVRHSFunctionfast_2WaySplit()
683 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
685 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionfast_2WaySplit()
686 slope[i * dof + j] = tmp; in FVRHSFunctionfast_2WaySplit()
694 uL = &ctx->uLR[0]; in FVRHSFunctionfast_2WaySplit()
695 uR = &ctx->uLR[dof]; in FVRHSFunctionfast_2WaySplit()
697 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
698 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionfast_2WaySplit()
699 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
701 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
703 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
708 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
709 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
710 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
712 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
714 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
717 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
722 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_2WaySplit()
723 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_2WaySplit()
724 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionfast_2WaySplit()
726 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_2WaySplit()
728 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_2WaySplit()
759 PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff)); in main()
760 PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming)); in main()
773 ctx.xmin = -1.0; in main()
776 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); in main()
777 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); in main()
778 …PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, si… in main()
779 …PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final… in main()
780 …PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given … in main()
781 …PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initia… in main()
782 …PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exa… in main()
783 …PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simula… in main()
784 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); in main()
785 …PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype,… in main()
786 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL)); in main()
807 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ in main()
813 …SetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax … in main()
827 … PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) s… in main()
828 count_fast = Mx - count_slow; in main()
834 if (((AdvectCtx *)ctx.physics2.user)->a > 0) { in main()
842 if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1) in main()
844 …else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwid… in main()
853 /* Create a time-stepping object */ in main()
880 const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow; in main()
881 const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast; in main()
892 if (i < ctx.sf || i > ctx.fs - 1) { in main()
906 mass_difference = mass_final - mass_initial; in main()
914 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
923 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); in main()
924 … PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); in main()
932 if (i < ctx.sf || i > ctx.fs - 1) in main()
933 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
935 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
939 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
951 PetscCall(VecAYPX(Y, -1, X)); in main()
995 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…
999 …args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…