VolumeIntegralAdaptive with IndicatorEntropyChange

examples/tree_2d_dgsem/elixir_euler_density_wave_adaptive_VI.jl

@@ -0,0 +1,73 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations2D(1.4)
+
+# We repeat test case for stability of EC fluxes from 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 
+# entropy-conserving flux-differencing volume integral.
+# If the entropy production of the weak form is lower than that of the
+# flux-differencing form, we use the flux-differencing form to stabilize the solution.
+indicator = IndicatorEntropyDiffusion(equations, basis)
+
+# Adaptive volume integral using the entropy production comparison indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff,
+                                         indicator = indicator)
+
+#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)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# 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)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9) # In paper, CFL = 0.05 is used
+
+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_VI.jl

@@ -0,0 +1,93 @@
+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 
+# entropy-conserving flux-differencing volume integral.
+# If the entropy production of the weak form is lower than that of the
+# flux-differencing form, we use the flux-differencing form to stabilize the solution.
+indicator = IndicatorEntropyDiffusion(equations, basis)
+
+# Adaptive volume integral using the entropy production comparison indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff,
+                                         indicator = indicator)
+
+#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)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# 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 = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -268,6 +268,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, VolumeIntegralShockCapturingRRG,
+       VolumeIntegralAdaptive, IndicatorEntropyDiffusion,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,

src/callbacks_step/analysis_dg2d.jl

@@ -185,6 +185,27 @@ function calc_error_norms(func, u, t, analyzer,
     return l2_error, linf_error
 end
 
+# used in `calc_entropy_change_element`
+function integrate_element_ref(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, 1, equations, dg, args...))
+
+    # Use quadrature to numerically integrate element.
+    # We do not multiply with 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 `calc_entropy_change_element`.
+    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
+
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}
@@ -214,8 +235,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/solvers/dg.jl

@@ -280,6 +280,58 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralAdaptive(;
+                           volume_integral_default = VolumeIntegralWeakForm(),
+                           volume_integral_stabilized = VolumeIntegralFluxDifferencing(flux_central),
+                           indicator = IndicatorEntropyDiffusion())
+
+!!! warning "Experimental code"
+    This code is experimental and may change in any future release.
+
+Possible combinations:
+- [`VolumeIntegralWeakForm`](@ref), [`VolumeIntegralFluxDifferencing`](@ref), and [`IndicatorEntropyDiffusion()`](@ref)
+"""
+struct VolumeIntegralAdaptive{VolumeIntegralDefault, VolumeIntegralStabilized,
+                              Indicator} <: AbstractVolumeIntegral
+    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
+    indicator::Indicator
+end
+
+function VolumeIntegralAdaptive(;
+                                volume_integral_default = VolumeIntegralWeakForm(),
+                                volume_integral_stabilized = VolumeIntegralFluxDifferencing(flux_central),
+                                indicator = IndicatorEntropyDiffusion())
+    return VolumeIntegralAdaptive{typeof(volume_integral_default),
+                                  typeof(volume_integral_stabilized),
+                                  typeof(indicator)}(volume_integral_default,
+                                                     volume_integral_stabilized,
+                                                     indicator)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain",
+                   integral::VolumeIntegralAdaptive)
+    @nospecialize integral # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        summary_header(io, "VolumeIntegralAdaptive")
+        summary_line(io, "volume integral default",
+                     integral.volume_integral_default)
+        summary_line(io, "volume integral stabilized",
+                     integral.volume_integral_stabilized)
+        if integral.indicator === nothing
+            summary_line(io, "indicator", integral.indicator)
+        else
+            summary_line(io, "indicator", integral.indicator |> typeof |> nameof)
+            show(increment_indent(io), mime, integral.indicator)
+        end
+        summary_footer(io)
+    end
+end
+
 # Abstract supertype for first-order `VolumeIntegralPureLGLFiniteVolume` and
 # second-order `VolumeIntegralPureLGLFiniteVolumeO2` subcell-based finite volume
 # volume integrals.

src/solvers/dgsem/calc_volume_integral.jl

@@ -12,6 +12,17 @@ function create_cache(mesh, equations,
     return NamedTuple()
 end
 
+function create_cache(mesh, equations,
+                      volume_integral::VolumeIntegralAdaptive,
+                      dg::DG, cache_containers, uEltype)
+    # This assumes that `volume_integral.volume_integral_default` needs no special cache!
+    @assert volume_integral.volume_integral_default isa VolumeIntegralWeakForm
+
+    return create_cache(mesh, equations,
+                        volume_integral.volume_integral_stabilized,
+                        dg, cache_containers, uEltype)
+end
+
 # The following `calc_volume_integral!` functions are
 # dimension and meshtype agnostic, i.e., valid for all 1D, 2D, and 3D meshes.
 
@@ -123,6 +134,68 @@ function calc_volume_integral!(du, u, mesh,
     return nothing
 end
 
+# Calculate ∫_e (∂S/∂u ⋅ ∂u/∂t) dΩ_e where "e" is an element
+function calc_entropy_change_element(du, u, element,
+                                     mesh::AbstractMesh{2}, equations, dg, cache)
+    return integrate_element_ref(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)
+        # Minus sign because of the flipped sign in the DG RHS.
+        # No scaling by inverse Jacobian here, as there is no Jacobian multiplication
+        # in `integrate_element_ref`.
+        -dot(cons2entropy(u_node, equations), du_node)
+    end
+end
+
+function calc_volume_integral!(du, u, mesh,
+                               have_nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralAdaptive{VolumeIntegralWeakForm,
+                                                                       VolumeIntegralFD,
+                                                                       Indicator},
+                               dg::DGSEM,
+                               cache) where {
+                                             VolumeIntegralFD <:
+                                             VolumeIntegralFluxDifferencing,
+                                             Indicator <: IndicatorEntropyDiffusion}
+    @unpack volume_integral_default, volume_integral_stabilized = volume_integral
+    @unpack du_element_threaded = volume_integral.indicator
+
+    @threaded for element in eachelement(dg, cache)
+        # Compute weak form volume integral
+        weak_form_kernel!(du, u, element, mesh,
+                          have_nonconservative_terms, equations,
+                          dg, cache)
+
+        # Compute entropy production of WF volume integral
+        entropy_delta_WF = calc_entropy_change_element(du, u, element,
+                                                       mesh, equations, dg, cache)
+        # Store weak form result
+        du_element_WF = du_element_threaded[Threads.threadid()]
+        @views du_element_WF .= du[.., element]
+
+        # Reset weak form volume integral 
+        du[.., element] .= zero(eltype(du))
+
+        # Recompute using entropy-conservative volume integral
+        flux_differencing_kernel!(du, u, element, mesh,
+                                  have_nonconservative_terms, equations,
+                                  volume_integral_stabilized.volume_flux,
+                                  dg, cache)
+
+        # Compute entropy production of FD volume integral
+        entropy_delta_FD = calc_entropy_change_element(du, u, element,
+                                                       mesh, equations, dg, cache)
+
+        entropy_delta = entropy_delta_WF - entropy_delta_FD
+        if entropy_delta < 0 # Use weak form if it is more stable
+            @views du[.., element] .= du_element_WF
+        end
+    end
+
+    return nothing
+end
+
 function calc_volume_integral!(du, u, mesh,
                                have_nonconservative_terms, equations,
                                volume_integral::VolumeIntegralPureLGLFiniteVolume,

src/solvers/dgsem/dgsem.jl

@@ -133,5 +133,6 @@ end
 end
 
 include("containers.jl")
+include("indicators.jl")
 include("calc_volume_integral.jl")
 end # @muladd