VolumeIntegralEntropyCorrection for TreeMesh

NEWS.md

@@ -20,6 +20,9 @@ This enables in particular adaptive mesh refinement for that solver-mesh combina
   `volume_integral_blend_high_order`, and `volume_integral_blend_low_order` besides the usual `indicator` argument.
   In particular, `volume_integral_default` may be e.g.  `VolumeIntegralWeakForm` or `VolumeIntegralAdaptive`, i.e.,
   the non-stabilized elements/cells are no longer restricted to `VolumeIntegralFluxDifferencing` only ([#2802]).
+- Added `IndicatorEntropyCorrection`. When combined with `VolumeIntegralAdaptive`, this blends together a stabilized and non-stabilized 
+  volume integral based on the violation of a volume entropy condition. `IndicatorEntropyCorrectionShockCapturingCombined` additionally 
+  guides the blending by taking the maximum of the entropy correction indicator and a shock capturing indicator ([#2764]). 
 - The second-order subcell volume integral is no longer limited to reconstruction in primitive variables.
   Instead, it is possible to reconstruct in custom variables, if functions `cons2recon` and `recon2cons` are provided to
   `VolumeIntegralPureLGLFiniteVolumeO2` and `VolumeIntegralShockCapturingRRG`([#2817]).

examples/tree_2d_dgsem/elixir_euler_peng_robinson_transcritical_mixing.jl → examples/structured_2d_dgsem/elixir_euler_peng_robinson_transcritical_mixing.jl

@@ -56,9 +56,15 @@ end
 initial_condition = initial_condition_transcritical_mixing
 
 volume_flux = flux_terashima_etal
-volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
-solver = DGSEM(polydeg = 3, volume_integral = volume_integral,
-               surface_flux = flux_lax_friedrichs)
+surface_flux = FluxPlusDissipation(volume_flux, DissipationLocalLaxFriedrichs())
+
+basis = LobattoLegendreBasis(3)
+volume_integral_default = VolumeIntegralFluxDifferencing(volume_flux)
+volume_integral_entropy_stable = VolumeIntegralPureLGLFiniteVolume(surface_flux)
+volume_integral = VolumeIntegralAdaptive(IndicatorEntropyCorrection(equations, basis),
+                                         volume_integral_default,
+                                         volume_integral_entropy_stable)
+solver = DGSEM(basis, surface_flux, volume_integral)
 
 cells_per_dimension = (32, 16)
 coordinates_min = (-0.5, -0.25)

examples/tree_1d_dgsem/elixir_euler_modified_sod_entropy_correction_amr.jl

@@ -0,0 +1,84 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations1D(1.4)
+
+function initial_condition_modified_sod(x, t, equations::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
+
+volume_flux = flux_central
+surface_flux = flux_lax_friedrichs
+basis = LobattoLegendreBasis(3)
+indicator_ec = IndicatorEntropyCorrection(equations, basis)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+indicator = IndicatorEntropyCorrectionShockCapturingCombined(indicator_entropy_correction = indicator_ec,
+                                                             indicator_shock_capturing = indicator_sc)
+
+volume_integral_default = VolumeIntegralFluxDifferencing(volume_flux)
+volume_integral_entropy_stable = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                                     volume_flux_fv = surface_flux)
+volume_integral = VolumeIntegralAdaptive(indicator,
+                                         volume_integral_default,
+                                         volume_integral_entropy_stable)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000,
+                periodicity = false)
+
+boundary_conditions = (x_neg = BoundaryConditionDirichlet(initial_condition),
+                       x_pos = BoundaryConditionDirichlet(initial_condition))
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 2000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorLöhner(semi, variable = first)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 3,
+                                      med_level = 4, med_threshold = 0.1,
+                                      max_level = 6, max_threshold = 0.15)
+amr_interval = 50
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = amr_interval)
+
+###############################################################################
+# run the simulation
+
+callbacks = CallbackSet(summary_callback, save_solution, amr_callback,
+                        analysis_callback, alive_callback)
+sol = solve(ode, SSPRK43();
+            abstol = 1e-6, reltol = 1e-4,
+            ode_default_options()..., callback = callbacks);

examples/tree_2d_dgsem/elixir_euler_nonideal_density_wave.jl

@@ -43,10 +43,12 @@ function Trixi.initial_condition_density_wave(x, t,
 end
 initial_condition = initial_condition_density_wave
 
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
 volume_flux = flux_terashima_etal
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
-solver = DGSEM(polydeg = 3, volume_integral = volume_integral,
-               surface_flux = flux_lax_friedrichs)
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)

src/Trixi.jl

@@ -273,6 +273,7 @@ export DG,
        VolumeIntegralAdaptive, IndicatorEntropyChange,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,
+       IndicatorEntropyCorrection, IndicatorEntropyCorrectionShockCapturingCombined,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
        MortarL2

src/auxiliary/math.jl

@@ -445,4 +445,8 @@ end
 @inline function maxmod(sl, sr)
     return 0.5f0 * (sign(sl) + sign(sr)) * max(abs(sl), abs(sr))
 end
+
+# Regularized approximation to the ratio a / b, which is numerically stable 
+# for b close to zero and approaches zero as b -> 0.
+@inline regularized_ratio(a, b) = a * b / (eps(typeof(b)) + b^2)
 end # @muladd

src/equations/nonideal_compressible_euler_1d.jl

@@ -328,4 +328,19 @@ end
     rho_e_internal = rho_e_total - 0.5f0 * rho_v1^2 / rho
     return rho_e_internal
 end
+
+"""
+    entropy_potential(u, orientation_or_normal_direction,
+                      equations::AbstractNonIdealCompressibleEulerEquations)
+
+Calculate the entropy potential, which for the compressible Euler equations with general
+EOS is ``p v_{\text{normal}} / T`` for the choice of [`entropy`](@ref) ``S(u) = -\rho s``. 
+"""
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::NonIdealCompressibleEulerEquations1D)
+    eos = equations.equation_of_state
+    V, v1, T = cons2thermo(u, equations)
+    p = pressure(V, T, eos)
+    return p * v1 / T
+end
 end # @muladd

src/equations/nonideal_compressible_euler_2d.jl

@@ -561,4 +561,25 @@ end
     rho_e_internal = rho_e_total - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho
     return rho_e_internal
 end
+
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::NonIdealCompressibleEulerEquations2D)
+    eos = equations.equation_of_state
+    V, v1, v2, T = cons2thermo(u, equations)
+    p = pressure(V, T, eos)
+    if orientation == 1
+        return p * v1 / T
+    else # if orientation == 2
+        return p * v2 / T
+    end
+end
+
+@inline function entropy_potential(u, normal_direction,
+                                   equations::NonIdealCompressibleEulerEquations2D)
+    eos = equations.equation_of_state
+    V, v1, v2, T = cons2thermo(u, equations)
+    v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
+    p = pressure(V, T, eos)
+    return p * v_normal / T
+end
 end # @muladd

src/solvers/dgsem/calc_volume_integral.jl

@@ -100,6 +100,83 @@ end
     return nothing
 end
 
+@inline function volume_integral_kernel!(du, u, element, mesh,
+                                         have_nonconservative_terms, equations,
+                                         volume_integral::VolumeIntegralEntropyCorrection,
+                                         dg::DGSEM, cache)
+    @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral
+    @unpack scaling = indicator
+    @unpack alpha = indicator.cache
+    du_element_threaded = indicator.cache.volume_integral_values_threaded
+
+    # run default volume integral 
+    volume_integral_kernel!(du, u, element, mesh,
+                            have_nonconservative_terms, equations,
+                            volume_integral_default, dg, cache)
+
+    # Check entropy production of "high order" volume integral. 
+    # 
+    # Note that, for `TreeMesh`, `dS_volume_integral` and `dS_true` are calculated
+    # on the reference element. For other mesh types, because ``dS_volume_integral`
+    # incorporates the scaled contravariant vectors, `dS_true` should 
+    # be calculated on the physical element instead.
+    #
+    # Minus sign because of the flipped sign of the volume term in the DG RHS.
+    # No scaling by inverse Jacobian here, as there is no Jacobian multiplication
+    # in `integrate_reference_element`.
+    dS_volume_integral = -entropy_change_reference_element(du, u, element,
+                                                           mesh, equations,
+                                                           dg, cache)
+
+    # Compute true entropy change given by surface integral of the entropy potential
+    dS_true = surface_integral_reference_element(entropy_potential, u, element,
+                                                 mesh, equations, dg, cache)
+
+    # This quantity should be ≤ 0 for an entropy stable volume integral, and 
+    # exactly zero for an entropy conservative volume integral. 
+    entropy_residual = dS_volume_integral - dS_true
+
+    if entropy_residual > 0
+        # Store "high order" result
+        du_FD_element = du_element_threaded[Threads.threadid()]
+        @views du_FD_element .= du[.., element]
+
+        # Reset pure flux-differencing volume integral 
+        # Note that this assumes that the volume terms are computed first,
+        # before any surface terms are added.
+        du[.., element] .= zero(eltype(du))
+
+        # Calculate entropy stable volume integral contribution
+        volume_integral_kernel!(du, u, element, mesh,
+                                have_nonconservative_terms, equations,
+                                volume_integral_stabilized, dg, cache)
+
+        dS_volume_integral_stabilized = -entropy_change_reference_element(du, u,
+                                                                          element,
+                                                                          mesh,
+                                                                          equations, dg,
+                                                                          cache)
+
+        # Calculate difference between high and low order FV entropy production;
+        # this should provide positive entropy dissipation if `entropy_residual > 0`, 
+        # assuming the stabilized volume integral is entropy stable.
+        entropy_dissipation = dS_volume_integral_stabilized - dS_volume_integral
+
+        # Calculate DG-FV blending factor 
+        ratio = regularized_ratio(-entropy_residual, entropy_dissipation)
+        alpha_element = min(1, scaling * ratio) # TODO: replacing this with a differentiable version of `min`
+
+        # Save blending coefficient for visualization
+        alpha[element] = alpha_element
+
+        # Blend the high order method back in 
+        @views du[.., element] .= alpha_element .* du[.., element] .+
+                                  (1 - alpha_element) .* du_FD_element
+    end
+
+    return nothing
+end
+
 function calc_volume_integral!(du, u, mesh,
                                have_nonconservative_terms, equations,
                                volume_integral, dg::DGSEM, cache)
@@ -153,4 +230,95 @@ function calc_volume_integral!(du, u, mesh,
 
     return nothing
 end
+
+function calc_volume_integral!(du, u, mesh,
+                               have_nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralEntropyCorrectionShockCapturingCombined,
+                               dg::DGSEM, cache)
+    (; volume_integral_default, volume_integral_stabilized, indicator) = volume_integral
+    (; indicator_entropy_correction, indicator_shock_capturing) = indicator
+    (; scaling) = indicator_entropy_correction
+    du_element_threaded = indicator_entropy_correction.cache.volume_integral_values_threaded
+
+    # Calculate DG-FV blending factors α a-priori for: u_{DG-FV} = u_DG * (1 - α) + u_FV * α
+    # Note that we also reuse the `alpha_shock_capturing` array to store the indicator values for visualization.
+    alpha_shock_capturing = @trixi_timeit timer() "blending factors" indicator_shock_capturing(u,
+                                                                                               mesh,
+                                                                                               equations,
+                                                                                               dg,
+                                                                                               cache)
+
+    @threaded for element in eachelement(dg, cache)
+        # run default volume integral 
+        volume_integral_kernel!(du, u, element, mesh,
+                                have_nonconservative_terms, equations,
+                                volume_integral_default, dg, cache)
+
+        # Check entropy production of "high order" volume integral. 
+        # 
+        # Note that, for `TreeMesh`, both volume and surface integrals are calculated
+        # on the reference element. For other mesh types, because the volume integral 
+        # incorporates the scaled contravariant vectors, the surface integral should 
+        # be calculated on the physical element instead.
+        #
+        # Minus sign because of the flipped sign of the volume term in the DG RHS.
+        # No scaling by inverse Jacobian here, as there is no Jacobian multiplication
+        # in `integrate_reference_element`.
+        dS_volume_integral = -entropy_change_reference_element(du, u, element,
+                                                               mesh, equations,
+                                                               dg, cache)
+
+        # Compute true entropy change given by surface integral of the entropy potential
+        dS_true = surface_integral_reference_element(entropy_potential, u, element,
+                                                     mesh, equations, dg, cache)
+
+        # This quantity should be ≤ 0 for an entropy stable volume integral, and 
+        # exactly zero for an entropy conservative volume integral. 
+        entropy_residual = dS_volume_integral - dS_true
+
+        if entropy_residual > 0
+            # Store "high order" result
+            du_FD_element = du_element_threaded[Threads.threadid()]
+            @views du_FD_element .= du[.., element]
+
+            # Reset pure flux-differencing volume integral 
+            # Note that this assumes that the volume terms are computed first,
+            # before any surface terms are added.
+            du[.., element] .= zero(eltype(du))
+
+            # Calculate entropy stable volume integral contribution
+            volume_integral_kernel!(du, u, element, mesh,
+                                    have_nonconservative_terms, equations,
+                                    volume_integral_stabilized, dg, cache)
+
+            dS_volume_integral_stabilized = -entropy_change_reference_element(du, u,
+                                                                              element,
+                                                                              mesh,
+                                                                              equations,
+                                                                              dg,
+                                                                              cache)
+
+            # Calculate difference between high and low order FV entropy production;
+            # this should provide positive entropy dissipation if `entropy_residual > 0`, 
+            # assuming the stabilized volume integral is entropy stable.
+            entropy_dissipation = dS_volume_integral_stabilized - dS_volume_integral
+
+            # Calculate DG-FV blending factor as the minimum between the entropy correction 
+            # indicator and shock capturing indicator
+            # TODO: replacing this with a differentiable version of `min`
+            ratio = regularized_ratio(-entropy_residual, entropy_dissipation)
+            alpha_element = min(1, max(alpha_shock_capturing[element], scaling * ratio))
+
+            # Save blending coefficient for visualization. Note that we overwrite the data 
+            # in `alpha_shock_capturing[element]`. 
+            alpha_shock_capturing[element] = alpha_element
+
+            # Blend the high order method back in 
+            @views du[.., element] .= alpha_element .* du[.., element] .+
+                                      (1 - alpha_element) .* du_FD_element
+        end
+    end
+
+    return nothing
+end
 end # @muladd

src/solvers/dgsem/dgsem.jl

@@ -73,7 +73,10 @@ end
 Base.summary(io::IO, dg::DGSEM) = print(io, "DGSEM(polydeg=$(polydeg(dg)))")
 
 include("compute_u_mean.jl")
+
 include("containers.jl")
+
 include("indicators.jl")
+include("special_volume_integrals.jl")
 include("calc_volume_integral.jl")
 end # @muladd