Lines Matching +full:- +full:j
1 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscret…
2 " advection - Constant coefficient scalar advection\n"
4 …oy problem, we choose different meshsizes for different sub-domains (slow-medium-fast-medium-slow)…
5 … " the meshsize ratio between two adjacient sub-domains is controlled with -hratio,\n"
6 …" exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect…
8 …" simulation - use reference solution which is generated by smaller time step size to be true so…
10 … " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
11 "Several initial conditions can be chosen with -initial N\n\n"
12 "The problem size should be set with -da_grid_x M\n\n";
22 PetscReal range = xmax - xmin; in RangeMod()
26 /* --------------------------------- Advection ----------------------------------- */
37 speed = ctx->a; in PhysicsRiemann_Advect()
50 speeds[0] = ctx->a; in PhysicsCharacteristic_Advect()
57 PetscReal a = ctx->a, x0; in PhysicsSample_Advect()
62 x0 = x - a * t; in PhysicsSample_Advect()
65 x0 = RangeMod(x - a * t, xmin, xmax); in PhysicsSample_Advect()
72 u[0] = (x0 < 0) ? 1 : -1; in PhysicsSample_Advect()
75 u[0] = (x0 < 0) ? -1 : 1; in PhysicsSample_Advect()
90 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); in PhysicsSample_Advect()
107 ctx->physics2.sample2 = PhysicsSample_Advect; in PhysicsCreate_Advect()
108 ctx->physics2.riemann2 = PhysicsRiemann_Advect; in PhysicsCreate_Advect()
109 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; in PhysicsCreate_Advect()
110 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; in PhysicsCreate_Advect()
111 ctx->physics2.user = user; in PhysicsCreate_Advect()
112 ctx->physics2.dof = 1; in PhysicsCreate_Advect()
113 PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0])); in PhysicsCreate_Advect()
114 user->a = 1; in PhysicsCreate_Advect()
115 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); in PhysicsCreate_Advect()
117 PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL)); in PhysicsCreate_Advect()
126 PetscInt i, j, k, dof, xs, xm, Mx; in FVSample_3WaySplit() local
131 …PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a samp… in FVSample_3WaySplit()
137 hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVSample_3WaySplit()
138 hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVSample_3WaySplit()
139 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVSample_3WaySplit()
141 if (i < ctx->sm) { in FVSample_3WaySplit()
142 xi = ctx->xmin + 0.5 * hs + i * hs; in FVSample_3WaySplit()
145 for (j = 0; j < N + 1; j++) { in FVSample_3WaySplit()
146 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_3WaySplit()
147 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_3WaySplit()
148 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_3WaySplit()
150 } else if (i < ctx->mf) { in FVSample_3WaySplit()
151 xi = ctx->xmin + ctx->sm * hs + 0.5 * hm + (i - ctx->sm) * hm; in FVSample_3WaySplit()
154 for (j = 0; j < N + 1; j++) { in FVSample_3WaySplit()
155 xj = xi + hm * (j - N / 2) / (PetscReal)N; in FVSample_3WaySplit()
156 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_3WaySplit()
157 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_3WaySplit()
159 } else if (i < ctx->fm) { in FVSample_3WaySplit()
160 xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + 0.5 * hf + (i - ctx->mf) * hf; in FVSample_3WaySplit()
163 for (j = 0; j < N + 1; j++) { in FVSample_3WaySplit()
164 xj = xi + hf * (j - N / 2) / (PetscReal)N; in FVSample_3WaySplit()
165 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_3WaySplit()
166 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_3WaySplit()
168 } else if (i < ctx->ms) { in FVSample_3WaySplit()
169 …xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + 0.5 * hm + (… in FVSample_3WaySplit()
172 for (j = 0; j < N + 1; j++) { in FVSample_3WaySplit()
173 xj = xi + hm * (j - N / 2) / (PetscReal)N; in FVSample_3WaySplit()
174 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_3WaySplit()
175 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_3WaySplit()
178 …xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + (ctx->ms - c… in FVSample_3WaySplit()
181 for (j = 0; j < N + 1; j++) { in FVSample_3WaySplit()
182 xj = xi + hs * (j - N / 2) / (PetscReal)N; in FVSample_3WaySplit()
183 …PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->… in FVSample_3WaySplit()
184 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; in FVSample_3WaySplit()
202 hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in SolutionErrorNorms_3WaySplit()
203 hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in SolutionErrorNorms_3WaySplit()
204 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in SolutionErrorNorms_3WaySplit()
210 if (i < ctx->sm || i > ctx->ms - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_3WaySplit()
211 else if (i < ctx->mf || i > ctx->fm - 1) *nrm1 += hm * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_3WaySplit()
212 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]); in SolutionErrorNorms_3WaySplit()
223 PetscInt i, j, k, Mx, dof, xs, xm, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->ms; in FVRHSFunction_3WaySplit() local
233 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunction_3WaySplit()
234 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunction_3WaySplit()
235 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunction_3WaySplit()
239 PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */ in FVRHSFunction_3WaySplit()
247 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunction_3WaySplit()
248 for (i = xs - 2; i < 0; i++) { in FVRHSFunction_3WaySplit()
249 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunction_3WaySplit()
252 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunction_3WaySplit()
255 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunction_3WaySplit()
258 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunction_3WaySplit()
259 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunction_3WaySplit()
260 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunction_3WaySplit()
261 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunction_3WaySplit()
262 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunction_3WaySplit()
263 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunction_3WaySplit()
264 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
266 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunction_3WaySplit()
267 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunction_3WaySplit()
269 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunction_3WaySplit()
270 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunction_3WaySplit()
278 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunction_3WaySplit()
279 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
281 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunction_3WaySplit()
282 slope[i * dof + j] = tmp; in FVRHSFunction_3WaySplit()
289 uL = &ctx->uLR[0]; in FVRHSFunction_3WaySplit()
290 uR = &ctx->uLR[dof]; in FVRHSFunction_3WaySplit()
292 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
293 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_3WaySplit()
294 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_3WaySplit()
296 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
298 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_3WaySplit()
301 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_3WaySplit()
304 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
305 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunction_3WaySplit()
306 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
308 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
310 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunction_3WaySplit()
313 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
316 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
317 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
318 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunction_3WaySplit()
320 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
322 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
325 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; in FVRHSFunction_3WaySplit()
328 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
329 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
330 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
332 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
334 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
337 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
340 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
341 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
342 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_3WaySplit()
344 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
346 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
349 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_3WaySplit()
352 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
353 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_3WaySplit()
354 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunction_3WaySplit()
356 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
358 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_3WaySplit()
361 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm; in FVRHSFunction_3WaySplit()
364 for (j = 0; j < dof; j++) { in FVRHSFunction_3WaySplit()
365 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunction_3WaySplit()
366 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunction_3WaySplit()
368 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunction_3WaySplit()
371 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunction_3WaySplit()
374 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; in FVRHSFunction_3WaySplit()
382 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunction_3WaySplit()
388 …ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\… in FVRHSFunction_3WaySplit()
393 /* --------------------------------- Finite Volume Solver for slow components ---------------------…
397 …PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbw… in FVRHSFunctionslow_3WaySplit() local
407 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunctionslow_3WaySplit()
408 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunctionslow_3WaySplit()
409 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunctionslow_3WaySplit()
418 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslow_3WaySplit()
419 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslow_3WaySplit()
420 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslow_3WaySplit()
423 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslow_3WaySplit()
426 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslow_3WaySplit()
429 …if (i < sm - lsbwidth + 1 || i > ms + rsbwidth - 2) { /* slow components and the first and last fa… in FVRHSFunctionslow_3WaySplit()
430 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslow_3WaySplit()
431 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslow_3WaySplit()
432 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslow_3WaySplit()
433 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslow_3WaySplit()
434 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslow_3WaySplit()
435 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslow_3WaySplit()
436 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
438 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslow_3WaySplit()
439 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslow_3WaySplit()
441 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslow_3WaySplit()
442 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslow_3WaySplit()
450 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunctionslow_3WaySplit()
451 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
453 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslow_3WaySplit()
454 slope[i * dof + j] = tmp; in FVRHSFunctionslow_3WaySplit()
462 uL = &ctx->uLR[0]; in FVRHSFunctionslow_3WaySplit()
463 uR = &ctx->uLR[dof]; in FVRHSFunctionslow_3WaySplit()
464 if (i < sm - lsbwidth) { /* slow region */ in FVRHSFunctionslow_3WaySplit()
465 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
466 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
467 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
469 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_3WaySplit()
472 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
475 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
479 if (i == sm - lsbwidth) { /* interface between slow and medium regions */ in FVRHSFunctionslow_3WaySplit()
480 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
481 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
482 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
484 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_3WaySplit()
486 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
490 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
491 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
492 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
494 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_3WaySplit()
496 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
501 for (j = 0; j < dof; j++) { in FVRHSFunctionslow_3WaySplit()
502 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
503 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslow_3WaySplit()
505 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslow_3WaySplit()
508 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
511 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslow_3WaySplit()
520 …PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((Pe… in FVRHSFunctionslow_3WaySplit()
527 …PetscInt i, j, k, Mx, dof, xs, xm, islowbuffer = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx… in FVRHSFunctionslowbuffer_3WaySplit() local
537 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunctionslowbuffer_3WaySplit()
538 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunctionslowbuffer_3WaySplit()
539 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunctionslowbuffer_3WaySplit()
548 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionslowbuffer_3WaySplit()
549 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionslowbuffer_3WaySplit()
550 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionslowbuffer_3WaySplit()
553 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionslowbuffer_3WaySplit()
556 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionslowbuffer_3WaySplit()
559 if ((i > sm - lsbwidth - 2 && i < sm + 1) || (i > ms - 2 && i < ms + rsbwidth + 1)) { in FVRHSFunctionslowbuffer_3WaySplit()
560 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionslowbuffer_3WaySplit()
561 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionslowbuffer_3WaySplit()
562 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionslowbuffer_3WaySplit()
563 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionslowbuffer_3WaySplit()
564 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionslowbuffer_3WaySplit()
565 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionslowbuffer_3WaySplit()
566 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
568 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionslowbuffer_3WaySplit()
569 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionslowbuffer_3WaySplit()
571 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionslowbuffer_3WaySplit()
572 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionslowbuffer_3WaySplit()
580 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunctionslowbuffer_3WaySplit()
581 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
583 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionslowbuffer_3WaySplit()
584 slope[i * dof + j] = tmp; in FVRHSFunctionslowbuffer_3WaySplit()
592 uL = &ctx->uLR[0]; in FVRHSFunctionslowbuffer_3WaySplit()
593 uR = &ctx->uLR[dof]; in FVRHSFunctionslowbuffer_3WaySplit()
594 if (i == sm - lsbwidth) { in FVRHSFunctionslowbuffer_3WaySplit()
595 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
596 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
597 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
599 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
601 for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
605 if (i > sm - lsbwidth && i < sm) { in FVRHSFunctionslowbuffer_3WaySplit()
606 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
607 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
608 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
610 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
612 for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
615 for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
620 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
621 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
622 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionslowbuffer_3WaySplit()
624 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
626 for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
630 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
631 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionslowbuffer_3WaySplit()
632 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
634 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
636 for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
641 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
642 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
643 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
645 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
647 for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
650 for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
655 for (j = 0; j < dof; j++) { in FVRHSFunctionslowbuffer_3WaySplit()
656 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
657 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionslowbuffer_3WaySplit()
659 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionslowbuffer_3WaySplit()
661 for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs; in FVRHSFunctionslowbuffer_3WaySplit()
672 /* --------------------------------- Finite Volume Solver for medium components -------------------…
676 …nt i, j, k, Mx, dof, xs, xm, imedium = 0, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->… in FVRHSFunctionmedium_3WaySplit() local
686 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunctionmedium_3WaySplit()
687 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunctionmedium_3WaySplit()
688 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunctionmedium_3WaySplit()
697 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionmedium_3WaySplit()
698 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionmedium_3WaySplit()
699 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionmedium_3WaySplit()
702 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionmedium_3WaySplit()
705 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionmedium_3WaySplit()
708 …if ((i > sm - 2 && i < mf - lmbwidth + 1) || (i > fm + rmbwidth - 2 && i < ms + 1)) { /* slow comp… in FVRHSFunctionmedium_3WaySplit()
709 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionmedium_3WaySplit()
710 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionmedium_3WaySplit()
711 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionmedium_3WaySplit()
712 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionmedium_3WaySplit()
713 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionmedium_3WaySplit()
714 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionmedium_3WaySplit()
715 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
717 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionmedium_3WaySplit()
718 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionmedium_3WaySplit()
720 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionmedium_3WaySplit()
721 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionmedium_3WaySplit()
729 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunctionmedium_3WaySplit()
730 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
732 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionmedium_3WaySplit()
733 slope[i * dof + j] = tmp; in FVRHSFunctionmedium_3WaySplit()
741 uL = &ctx->uLR[0]; in FVRHSFunctionmedium_3WaySplit()
742 uR = &ctx->uLR[dof]; in FVRHSFunctionmedium_3WaySplit()
744 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
745 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; in FVRHSFunctionmedium_3WaySplit()
746 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
748 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
750 for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
754 if (i > sm && i < mf - lmbwidth) { /* medium region */ in FVRHSFunctionmedium_3WaySplit()
755 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
756 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
757 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
759 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
761 for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
764 for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
768 if (i == mf - lmbwidth) { /* interface between medium and fast regions */ in FVRHSFunctionmedium_3WaySplit()
769 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
770 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
771 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
773 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
775 for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
779 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
780 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
781 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
783 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
785 for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
790 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
791 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
792 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
794 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
796 for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
799 for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
804 for (j = 0; j < dof; j++) { in FVRHSFunctionmedium_3WaySplit()
805 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmedium_3WaySplit()
806 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; in FVRHSFunctionmedium_3WaySplit()
808 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmedium_3WaySplit()
810 for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmedium_3WaySplit()
824 …PetscInt i, j, k, Mx, dof, xs, xm, imediumbuffer = 0, mf = ctx->mf, fm = ctx->fm, lmbwidth = c… in FVRHSFunctionmediumbuffer_3WaySplit() local
834 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunctionmediumbuffer_3WaySplit()
835 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunctionmediumbuffer_3WaySplit()
836 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunctionmediumbuffer_3WaySplit()
845 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionmediumbuffer_3WaySplit()
846 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionmediumbuffer_3WaySplit()
847 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionmediumbuffer_3WaySplit()
850 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionmediumbuffer_3WaySplit()
853 for (i = xs - 1; i < xs + xm + 1; i++) { in FVRHSFunctionmediumbuffer_3WaySplit()
856 if ((i > mf - lmbwidth - 2 && i < mf + 1) || (i > fm - 2 && i < fm + rmbwidth + 1)) { in FVRHSFunctionmediumbuffer_3WaySplit()
857 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ in FVRHSFunctionmediumbuffer_3WaySplit()
858 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionmediumbuffer_3WaySplit()
859 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ in FVRHSFunctionmediumbuffer_3WaySplit()
860 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionmediumbuffer_3WaySplit()
861 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionmediumbuffer_3WaySplit()
862 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionmediumbuffer_3WaySplit()
863 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
865 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionmediumbuffer_3WaySplit()
866 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionmediumbuffer_3WaySplit()
868 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionmediumbuffer_3WaySplit()
869 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionmediumbuffer_3WaySplit()
877 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunctionmediumbuffer_3WaySplit()
878 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
880 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionmediumbuffer_3WaySplit()
881 slope[i * dof + j] = tmp; in FVRHSFunctionmediumbuffer_3WaySplit()
889 uL = &ctx->uLR[0]; in FVRHSFunctionmediumbuffer_3WaySplit()
890 uR = &ctx->uLR[dof]; in FVRHSFunctionmediumbuffer_3WaySplit()
891 if (i == mf - lmbwidth) { in FVRHSFunctionmediumbuffer_3WaySplit()
892 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
893 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
894 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
896 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
898 for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
902 if (i > mf - lmbwidth && i < mf) { in FVRHSFunctionmediumbuffer_3WaySplit()
903 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
904 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
905 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
907 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
909 for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
912 for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
917 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
918 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
919 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
921 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
923 for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
927 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
928 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
929 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
931 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
933 for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
938 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
939 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
940 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
942 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
944 for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
947 for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
952 for (j = 0; j < dof; j++) { in FVRHSFunctionmediumbuffer_3WaySplit()
953 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
954 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionmediumbuffer_3WaySplit()
956 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionmediumbuffer_3WaySplit()
958 for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm; in FVRHSFunctionmediumbuffer_3WaySplit()
969 /* --------------------------------- Finite Volume Solver for fast parts -------------------------…
973 PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, mf = ctx->mf, fm = ctx->fm; in FVRHSFunctionfast_3WaySplit() local
983 hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm; in FVRHSFunctionfast_3WaySplit()
984 hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm); in FVRHSFunctionfast_3WaySplit()
985 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf); in FVRHSFunctionfast_3WaySplit()
994 if (ctx->bctype == FVBC_OUTFLOW) { in FVRHSFunctionfast_3WaySplit()
995 for (i = xs - 2; i < 0; i++) { in FVRHSFunctionfast_3WaySplit()
996 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; in FVRHSFunctionfast_3WaySplit()
999 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; in FVRHSFunctionfast_3WaySplit()
1002 …for (i = xs - 1; i < xs + xm + 1; i++) { /* fast components and the last slow components before fa… in FVRHSFunctionfast_3WaySplit()
1005 if (i > mf - 2 && i < fm + 1) { in FVRHSFunctionfast_3WaySplit()
1006 …PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv… in FVRHSFunctionfast_3WaySplit()
1007 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); in FVRHSFunctionfast_3WaySplit()
1008 cjmpL = &ctx->cjmpLR[0]; in FVRHSFunctionfast_3WaySplit()
1009 cjmpR = &ctx->cjmpLR[dof]; in FVRHSFunctionfast_3WaySplit()
1010 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_3WaySplit()
1012 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; in FVRHSFunctionfast_3WaySplit()
1013 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; in FVRHSFunctionfast_3WaySplit()
1015 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; in FVRHSFunctionfast_3WaySplit()
1016 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; in FVRHSFunctionfast_3WaySplit()
1024 (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope); in FVRHSFunctionfast_3WaySplit()
1025 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_3WaySplit()
1027 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; in FVRHSFunctionfast_3WaySplit()
1028 slope[i * dof + j] = tmp; in FVRHSFunctionfast_3WaySplit()
1036 uL = &ctx->uLR[0]; in FVRHSFunctionfast_3WaySplit()
1037 uR = &ctx->uLR[dof]; in FVRHSFunctionfast_3WaySplit()
1039 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_3WaySplit()
1040 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2; in FVRHSFunctionfast_3WaySplit()
1041 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_3WaySplit()
1043 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_3WaySplit()
1045 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_3WaySplit()
1050 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_3WaySplit()
1051 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_3WaySplit()
1052 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; in FVRHSFunctionfast_3WaySplit()
1054 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_3WaySplit()
1056 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_3WaySplit()
1059 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; in FVRHSFunctionfast_3WaySplit()
1064 for (j = 0; j < dof; j++) { in FVRHSFunctionfast_3WaySplit()
1065 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; in FVRHSFunctionfast_3WaySplit()
1066 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2; in FVRHSFunctionfast_3WaySplit()
1068 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); in FVRHSFunctionfast_3WaySplit()
1070 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; in FVRHSFunctionfast_3WaySplit()
1101 PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit3_LaxWendroff)); in main()
1102 PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit3_BeamWarming)); in main()
1115 ctx.xmin = -1.0; in main()
1118 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); in main()
1119 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); in main()
1120 …PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, si… in main()
1121 …PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final… in main()
1122 …PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given … in main()
1123 …PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initia… in main()
1124 …PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exa… in main()
1125 …PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simula… in main()
1126 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); in main()
1127 …PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype,… in main()
1128 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL)); in main()
1149 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ in main()
1155 …SetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax … in main()
1170 … PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) s… in main()
1181 if (((AdvectCtx *)ctx.physics2.user)->a > 0) { in main()
1194 if (i < ctx.sm - ctx.lsbwidth || i > ctx.ms + ctx.rsbwidth - 1) in main()
1196 …else if ((i >= ctx.sm - ctx.lsbwidth && i < ctx.sm) || (i > ctx.ms - 1 && i <= ctx.ms + ctx.rsbwid… in main()
1198 else if (i < ctx.mf - ctx.lmbwidth || i > ctx.fm + ctx.rmbwidth - 1) in main()
1200 …else if ((i >= ctx.mf - ctx.lmbwidth && i < ctx.mf) || (i > ctx.fm - 1 && i <= ctx.fm + ctx.rmbwid… in main()
1211 /* Create a time-stepping object */ in main()
1242 const PetscReal hs = (ctx.xmax - ctx.xmin) / 4.0 / count_slow; in main()
1243 const PetscReal hm = (ctx.xmax - ctx.xmin) / 2.0 / count_medium; in main()
1244 const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast; in main()
1255 if (i < ctx.sm || i > ctx.ms - 1) in main()
1260 else if (i < ctx.mf || i > ctx.fm - 1) in main()
1274 mass_difference = mass_final - mass_initial; in main()
1282 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
1291 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); in main()
1292 … PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); in main()
1300 if (i < ctx.sm || i > ctx.ms - 1) in main()
1301 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
1302 else if (i < ctx.mf || i < ctx.fm - 1) in main()
1303 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hm * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
1305 … for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); in main()
1309 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); in main()
1321 PetscCall(VecAYPX(Y, -1, X)); in main()
1369 …args: -da_grid_x 90 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…
1373 …args: -da_grid_x 90 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_st…