VolumeIntegralAdaptive with IndicatorEntropyChange

examples/tree_2d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl

@@ -0,0 +1,77 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations2D(1.4)
+
+# We repeat test case for linear stability of EC fluxes from the paper
+# - Gregor J. Gassner, Magnus Svärd, Florian J. Hindenlang (2020)
+#   Stability issues of entropy-stable and/or split-form high-order schemes
+#   [DOI: 10.1007/s10915-021-01720-8](https://doi.org/10.1007/s10915-021-01720-8)
+initial_condition = initial_condition_density_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_chandrashekar
+
+polydeg = 5
+basis = LobattoLegendreBasis(polydeg)
+
+volume_integral_weakform = VolumeIntegralWeakForm()
+volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(volume_flux)
+
+# This indicator compares the entropy production of the weak form to the 
+# true entropy evolution in that cell.
+# If the weak form dissipates more entropy than the true evolution
+# the indicator renders this admissible. Otherwise, the more stable
+# volume integral is to be used.
+indicator = IndicatorEntropyChange(maximum_entropy_increase = 0.0)
+
+# Adaptive volume integral using the entropy change indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(indicator = indicator,
+                                         volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff)
+
+#volume_integral = volume_integral_weakform # Stable, but unphysical entropy increase!
+#volume_integral = volume_integral_fluxdiff # Crashes!
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2, # 4 x 4 elements
+                n_cells_max = 30_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+# In the paper, CFL = 0.05 is used. Same observations recovered for CFL = 0.9, though.
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, ode_default_options()..., callback = callbacks);

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl

@@ -0,0 +1,96 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    RealT = eltype(x)
+    slope = 15
+    B = tanh(slope * x[2] + 7.5f0) - tanh(slope * x[2] - 7.5f0)
+    rho = 0.5f0 + 0.75f0 * B
+    v1 = 0.5f0 * (B - 1)
+    v2 = convert(RealT, 0.1) * sinpi(2 * x[1])
+    p = 1
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+volume_integral_weakform = VolumeIntegralWeakForm()
+volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(volume_flux)
+
+# This indicator compares the entropy production of the weak form to the 
+# true entropy evolution in that cell.
+# If the weak form dissipates more entropy than the true evolution
+# the indicator renders this admissible. Otherwise, the more stable
+# volume integral is to be used.
+indicator = IndicatorEntropyChange(maximum_entropy_increase = 0.0)
+
+# Adaptive volume integral using the entropy change indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(indicator = indicator,
+                                         volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff)
+
+#volume_integral = volume_integral_weakform # Crashes
+#volume_integral = volume_integral_fluxdiff # Crashes as well!
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 100_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.25)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_errors = Symbol[],
+                                     extra_analysis_integrals = (entropy,),
+                                     save_analysis = true)
+
+alive_callback = AliveCallback(alive_interval = 200)
+
+stepsize_callback = StepsizeCallback(cfl = 1.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 5e-3, ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -251,6 +251,7 @@ export density, pressure, density_pressure, velocity, temperature,
        equilibrium_distribution,
        waterheight, waterheight_pressure
 export entropy, entropy_thermodynamic, entropy_math, entropy_guermond_etal,
+       entropy_potential,
        energy_total, energy_kinetic, energy_internal, energy_internal_specific,
        energy_magnetic, cross_helicity, magnetic_field, divergence_cleaning_field,
        enstrophy, vorticity
@@ -268,6 +269,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, VolumeIntegralShockCapturingRRG,
+       VolumeIntegralAdaptive, IndicatorEntropyChange,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,

src/callbacks_step/analysis_dg2d.jl

@@ -185,6 +185,68 @@ function calc_error_norms(func, u, t, analyzer,
     return l2_error, linf_error
 end
 
+# Use quadrature to numerically integrate a single element.
+# We do not multiply by the Jacobian to stay in reference space.
+# This avoids the need to divide the RHS of the DG scheme by the Jacobian when computing
+# the time derivative of entropy, see `entropy_change_reference_element`.
+function integrate_reference_element(func::Func, u, element,
+                                     mesh::AbstractMesh{2}, equations, dg::DGSEM, cache,
+                                     args...) where {Func}
+    @unpack weights = dg.basis
+
+    # Initialize integral with zeros of the right shape
+    element_integral = zero(func(u, 1, 1, element, equations, dg, args...))
+
+    for j in eachnode(dg), i in eachnode(dg)
+        element_integral += weights[i] * weights[j] *
+                            func(u, i, j, element, equations, dg, args...)
+    end
+
+    return element_integral
+end
+
+# Calculate ∫_e (∂S/∂u ⋅ ∂u/∂t) dΩ_e where the result on element 'e' is kept in reference space
+# Note that ∂S/∂u = w(u) with entropy variables w
+function entropy_change_reference_element(du::AbstractArray{<:Any, 4}, u,
+                                          element,
+                                          mesh::AbstractMesh{2},
+                                          equations, dg::DGSEM, cache, args...)
+    return integrate_reference_element(u, element, mesh, equations, dg, cache,
+                                       du) do u, i, j, element, equations, dg, du
+        u_node = get_node_vars(u, equations, dg, i, j, element)
+        du_node = get_node_vars(du, equations, dg, i, j, element)
+
+        dot(cons2entropy(u_node, equations), du_node)
+    end
+end
+
+# calculate surface integral of func(u, equations) * normal on the reference element.
+function surface_integral(func::Func, u, element,
+                          mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
+                          args...) where {Func}
+    @unpack weights = dg.basis
+
+    u_tmp = get_node_vars(u, equations, dg, 1, 1, element)
+    surface_integral = zero(func(u_tmp, 1, equations))
+    for i in eachnode(dg)
+        # integrate along x direction, normal in y (2) direction
+        u_bottom = get_node_vars(u, equations, dg, i, 1, element)
+        u_top = get_node_vars(u, equations, dg, i, nnodes(dg), element)
+
+        surface_integral += weights[i] *
+                            (func(u_top, 2, equations) - func(u_bottom, 2, equations))
+
+        # integrate along y direction, normal in x (1) direction
+        u_left = get_node_vars(u, equations, dg, 1, i, element)
+        u_right = get_node_vars(u, equations, dg, nnodes(dg), i, element)
+
+        surface_integral += weights[i] *
+                            (func(u_right, 1, equations) - func(u_left, 1, equations))
+    end
+
+    return surface_integral
+end
+
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}
@@ -214,8 +276,8 @@ function integrate_via_indices(func::Func, u,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
                                            UnstructuredMesh2D, P4estMesh{2},
                                            T8codeMesh{2}},
-                               equations,
-                               dg::DGSEM, cache, args...; normalize = true) where {Func}
+                               equations, dg::DGSEM, cache,
+                               args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
     # Initialize integral with zeros of the right shape

src/equations/compressible_euler_2d.jl

@@ -2235,6 +2235,28 @@ end
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
 
+@doc raw"""
+    entropy_potential(u, orientation::Int, 
+                      equations::AbstractCompressibleEulerEquations)
+
+Calculate the entropy potential, which for the compressible Euler equations is simply 
+the momentum for the choice of mathematical [`entropy`](@ref) ``S(u) = \frac{\rho s}{\gamma - 1}``
+with thermodynamic entropy ``s = \ln(p) - \gamma \ln(\rho)``.
+
+## References
+- Eitan Tadmor (2003)
+  Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems
+  [DOI: 10.1017/S0962492902000156](https://doi.org/10.1017/S0962492902000156)
+"""
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::CompressibleEulerEquations2D)
+    if orientation == 1
+        return u[2]
+    else # if orientation == 2
+        return u[3]
+    end
+end
+
 # State validation for Newton-bisection method of subcell IDP limiting
 @inline function Base.isvalid(u, equations::CompressibleEulerEquations2D)
     if u[1] <= 0 || pressure(u, equations) <= 0

src/solvers/dg.jl

@@ -282,6 +282,80 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralAdaptive(;
+                           indicator = IndicatorEntropyChange(),
+                           volume_integral_default,
+                           volume_integral_stabilized)
+
+!!! warning "Experimental code"
+    This code is experimental and may change in any future release.
+
+The `indicator` is currently limited to [`IndicatorEntropyChange`](@ref).
+"""
+struct VolumeIntegralAdaptive{Indicator,
+                              VolumeIntegralDefault, VolumeIntegralStabilized} <:
+       AbstractVolumeIntegral
+    indicator::Indicator # A-posteriori indicator called after computation of `volume_integral_default`
+    volume_integral_default::VolumeIntegralDefault # Cheap(er) default volume integral to be used in non-critical regions
+    volume_integral_stabilized::VolumeIntegralStabilized # More expensive volume integral with stabilizing effect
+end
+
+function VolumeIntegralAdaptive(;
+                                indicator = IndicatorEntropyChange(),
+                                volume_integral_default,
+                                volume_integral_stabilized)
+    if !(indicator isa IndicatorEntropyChange)
+        throw(ArgumentError("`indicator` must be of type `IndicatorEntropyChange`."))
+    end
+
+    return VolumeIntegralAdaptive{typeof(indicator),
+                                  typeof(volume_integral_default),
+                                  typeof(volume_integral_stabilized)}(indicator,
+                                                                      volume_integral_default,
+                                                                      volume_integral_stabilized)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain",
+                   integral::VolumeIntegralAdaptive)
+    @nospecialize integral # reduce precompilation time
+    @unpack volume_integral_default, volume_integral_stabilized, indicator = integral
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        summary_header(io, "VolumeIntegralAdaptive")
+
+        summary_line(io, "volume integral default",
+                     volume_integral_default |> typeof |> nameof)
+        if !(volume_integral_default isa VolumeIntegralWeakForm)
+            show(increment_indent(io), mime, volume_integral_default)
+        end
+
+        summary_line(io, "volume integral stabilized",
+                     volume_integral_stabilized |> typeof |> nameof)
+        show(increment_indent(io), mime, volume_integral_stabilized)
+
+        summary_line(io, "indicator", indicator |> typeof |> nameof)
+        show(increment_indent(io), mime, indicator)
+
+        summary_footer(io)
+    end
+end
+
+function create_cache(mesh, equations,
+                      volume_integral::VolumeIntegralAdaptive,
+                      dg, cache_containers, uEltype)
+    cache_default = create_cache(mesh, equations,
+                                 volume_integral.volume_integral_default,
+                                 dg, cache_containers, uEltype)
+    cache_stabilized = create_cache(mesh, equations,
+                                    volume_integral.volume_integral_stabilized,
+                                    dg, cache_containers, uEltype)
+
+    return (; cache_default..., cache_stabilized...)
+end
+
 # Abstract supertype for first-order `VolumeIntegralPureLGLFiniteVolume` and
 # second-order `VolumeIntegralPureLGLFiniteVolumeO2` subcell-based finite volume
 # volume integrals.