You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
The positivity preserving feature PositivityPreservingLimiterZhangShu does not seem to work for GaussSBP() with the time stepper SSPRK43() or CarpenterKennedy2N54
1D Compressible Euler example with GaussSBP():
using Trixi
using OrdinaryDiffEq
gamma_gas = 1.4
equations = CompressibleEulerEquations1D(gamma_gas)
###############################################################################
# setup the GSBP DG discretization that uses the Gauss operators from
# Chan, Del Rey Fernandez, Carpenter (2019).
# [https://doi.org/10.1137/18M1209234](https://doi.org/10.1137/18M1209234)
# Shu-Osher initial condition for 1D compressible Euler equations
# Example 8, https://www.sciencedirect.com/science/article/pii/0021999189902222
function initial_condition_shu_osher(x, t, equations::CompressibleEulerEquations1D)
x0 = -4
rho_left = 27/7
v_left = 4*sqrt(35)/9
p_left = 31/3
# Replaced v_right = 0 to v_right = 0.1 to avoid positivity issues.
v_right = 0.1
p_right = 1.0
rho = ifelse(x[1] > x0, 1 + 1/5*sin(5*x[1]), rho_left)
v = ifelse(x[1] > x0, v_right, v_left)
p = ifelse(x[1] > x0, p_right, p_left)
return prim2cons(SVector(rho + 0.1*rand(), v + 0.1*rand(), p + 0.1*rand()), equations)
end
initial_condition = initial_condition_shu_osher
surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
polydeg = 3
basis = DGMultiBasis(Line(), polydeg, approximation_type = GaussSBP())
indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)
dg = DGMulti(basis,
surface_integral = SurfaceIntegralWeakForm(surface_flux),
volume_integral = volume_integral)
boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (; :entire_boundary => boundary_condition)
###############################################################################
# setup the 1D mesh
cells_per_dimension = (64,)
mesh = DGMultiMesh(dg, cells_per_dimension,
coordinates_min = (-5.0,), coordinates_max = (5.0,),
periodicity = false)
###############################################################################
# setup the semidiscretization and ODE problem
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
dg, boundary_conditions = boundary_conditions)
tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)
###############################################################################
# setup the callbacks
# prints a summary of the simulation setup and resets the timers
summary_callback = SummaryCallback()
# analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval = 100, uEltype = real(dg))
# handles the re-calculation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl = 0.1)
# collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
# For future work, add PositivityPreservingLimiterZhangShu for DGMulti with Gauss nodes
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
variables=(Trixi.density, pressure))
# ###############################################################################
# # run the simulation
sol = solve(ode, SSPRK43(stage_limiter!), adaptive = true, callback = callbacks, save_everystep = true)
The text was updated successfully, but these errors were encountered:
As I understand, the arguments in Zhang and Shu's positivity limiter only apply when solution values at Gauss-Legendre-Lobatto (GLL) points (or any solution points that include $-1,1$) are made admissibility. One way to apply it to a DG scheme using any other solution points is to (1) Extrapolate/Interpolate the solution polynomial to GLL points, and then (2) Correct those GLL point values.
Yes - the Zhang-Shu limiter involves two steps, not all of which generalize to Gauss-Legendre solvers. These are:
a CFL condition to ensure that the cell average of the high order DG solution is positive
a scaling limiter to compress the high order solution towards the cell solution average
The first argument assumes both a Forward Euler time-step structure and "standard" polynomial interpolation to evaluate interface values of the solution, neither of which is true for GaussSBP solvers, so I'm not sure we could prove positivity of Zhang-Shu limiting in this setting.
We could, however, implement the second part - a heuristic a-posteriori scaling to try to avoid negativity.
Hello,
The positivity preserving feature
PositivityPreservingLimiterZhangShu
does not seem to work forGaussSBP()
with the time stepperSSPRK43()
orCarpenterKennedy2N54
1D Compressible Euler example with
GaussSBP()
:The text was updated successfully, but these errors were encountered: