Skip to content

flux_ranocha is not entropy conservative for CompressibleEulerMulticomponentEquations1D #3148

Description

@jlchan

The implementation of flux_ranocha for CompressibleEulerMulticomponentEquations1D is not entropy conservative when gammas or gas_constants are not identical. flux_chandrashekar seems OK though.

Additionally, the initial_condition_weak_blast_wave used in examples/tree_1d_dgsem/elixir_eulermulti_ec.jl results in rho1, rho2, rho3 which are identical up to a constant scaling - this also seems to preserve EC up to machine precision.

Here is a MWE where rho_i are not scalings of each other and gammas/gas_constants are not identical.

using OrdinaryDiffEqLowStorageRK
using Trixi

equations = CompressibleEulerMulticomponentEquations1D(gammas = (1.4, 5/3, 4/3),
                                                       gas_constants = (0.4, 0.6, 0.5))

function initial_condition_random(x, t, equations)
    rho1, rho2, rho3 = ntuple(i -> 1 + 0.1 * sin(i + x[1]), 3)
    v1 = 0.1 
    p = 1 + 0.1 * sin(x[1])
    q = SVector(v1, p, rho1, rho2, rho3)
    return prim2cons(q, equations)
end
initial_condition = initial_condition_random

volume_flux = flux_ranocha
solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-2.0,)
coordinates_max = (2.0,)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 10_000, periodicity = true)

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

tspan = (0.0, 1e-2)
ode = semidiscretize(semi, tspan)

analysis_interval = 10
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                     extra_analysis_integrals = (Trixi.density,))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 0.1)

callbacks = CallbackSet(analysis_callback, alive_callback,                        
                        stepsize_callback)

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

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
            dt = 1, # solve needs some value here but it will be overwritten by the stepsize_callback
            ode_default_options()..., callback = callbacks);

The analysis output yields ∑∂S/∂U ⋅ Uₜ : -6.80429623e-06

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions