elixir_euler_modified_sod.jl

using OrdinaryDiffEqSSPRK
using Trixi

equations = CompressibleEulerEquations1D(1.4)

"""
    initial_condition_modified_sod(x, t, equations::CompressibleEulerEquations1D)

Modified Sod shock tube problem, presented in Section 6.4 of Toro's book.
This problem consists of a left sonic rarefaction wave and is useful for testing whether numerical solutions
violate the entropy condition.
An entropy-satisfying solution should produce a smooth(!) rarefaction wave.

## References
- Toro (2009).
  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 3rd Edition.
  [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)

- Lin, Chan (2014)
  High order entropy stable discontinuous Galerkin spectral element methods through subcell limiting
  [DOI: 10.1016/j.jcp.2023.112677](https://doi.org/10.1016/j.jcp.2023.112677)
"""
function initial_condition_modified_sod(x, t, ::CompressibleEulerEquations1D)
    if x[1] < 0.3
        return prim2cons(SVector(1, 0.75, 1), equations)
    else
        return prim2cons(SVector(0.125, 0.0, 0.1), equations)
    end
end
initial_condition = initial_condition_modified_sod

# Using the weak form volume integral gives a wrong solution at the rarefaction wave!
#volume_integral = VolumeIntegralWeakForm()

# The entropy-conservative flux-differencing volume integral recovers the rarefaction wave correctly!
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)

solver = DGSEM(polydeg = 3, surface_flux = flux_hllc,
               volume_integral = volume_integral)

coordinates_min = 0.0
coordinates_max = 1.0

mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 6,
                n_cells_max = 30_000,
                periodicity = false)

# Dirichlet boundary condition is only valid for considered time interval.
# If the rarefaction wave reaches the boundary, this condition is no longer valid!
boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition),
                       x_pos = boundary_condition_do_nothing)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions = boundary_conditions)

tspan = (0.0, 0.2)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval = 100)

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

###############################################################################
# run the simulation

# For weak form run: Need to enforce positivity explicitly! Not required for flux differencing volume integral.
#=
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
                                                     variables = (Trixi.density, pressure))

ode_alg = SSPRK43(stage_limiter! = stage_limiter!)
=#
# Flux-differencing volume integral does not require positivity preservation for this test case.
ode_alg = SSPRK43()

sol = solve(ode, ode_alg;
            dt = 4e-4, adaptive = true,
            ode_default_options()..., callback = callbacks);