VolumeIntegralEntropyCorrection for TreeMesh

examples/tree_1d_dgsem/elixir_euler_modified_sod_entropy_correction.jl

@@ -0,0 +1,64 @@
+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
+        # this version of modified sod uses a 100x density and pressure jump
+        return prim2cons(SVector(0.0125, 0.0, 0.01), equations)
+    end
+end
+
+initial_condition = initial_condition_modified_sod
+
+volume_flux = flux_central
+surface_flux = flux_lax_friedrichs
+basis = LobattoLegendreBasis(3)
+volume_integral = VolumeIntegralEntropyCorrection(equations, basis;
+                                                  volume_flux_dg = volume_flux,
+                                                  volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, 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)
+
+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)
+
+###############################################################################
+# run the simulation
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback)
+sol = solve(ode, SSPRK43();
+            abstol = 1e-6, reltol = 1e-4,
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -252,7 +252,7 @@ export density, pressure, density_pressure, velocity, temperature,
        equilibrium_distribution,
        waterheight, waterheight_pressure
 export entropy, entropy_thermodynamic, entropy_math, entropy_guermond_etal,
-       energy_total, energy_kinetic, energy_internal,
+       energy_total, energy_kinetic, energy_internal, entropy_potential,
        energy_magnetic, cross_helicity, magnetic_field, divergence_cleaning_field,
        enstrophy, vorticity
 export lake_at_rest_error
@@ -271,6 +271,7 @@ export DG,
        VolumeIntegralShockCapturingHG, VolumeIntegralShockCapturingRRG,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,
+       VolumeIntegralEntropyCorrection,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
        MortarL2
@@ -317,7 +318,8 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
 
 export ControllerThreeLevel, ControllerThreeLevelCombined,
-       IndicatorLöhner, IndicatorLoehner, IndicatorMax
+       IndicatorLöhner, IndicatorLoehner, IndicatorMax,
+       IndicatorEntropyCorrection
 
 export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter
 

src/callbacks_step/analysis_dg1d.jl

@@ -119,6 +119,37 @@ function calc_error_norms(func, u, t, analyzer,
     return l2_error, linf_error
 end
 
+# integrate `du_local` against the entropy variables evaluated at `u`
+function integrate_against_entropy_variables(du_local, u, element,
+                                             mesh::AbstractMesh{1},
+                                             equations, dg::DGSEM, cache, args...)
+    (; weights) = dg.basis
+
+    # we integrate over the reference element since `du` is calculated by 
+    # a weak form-like volume integral, and should already be scaled by the 
+    # Jacobian factor. 
+
+    integral = zero(eltype(u))
+
+    for i in eachnode(dg)
+        du_node = get_node_vars(du_local, equations, dg, i)
+        u_node = get_node_vars(u, equations, dg, i, element)
+        integral = integral +
+                   weights[i] * dot(cons2entropy(u_node, equations), du_node)
+    end
+
+    return integral
+end
+
+# calculate surface integral of func(u, orientation, equations) * normal
+function surface_integral(func::Func, u, element, mesh::TreeMesh{1}, equations,
+                          dg::DGSEM, cache) where {Func}
+    u_left = get_node_vars(u, equations, dg, 1, element)
+    u_right = get_node_vars(u, equations, dg, nnodes(dg), element)
+    surface_integral = (func(u_right, 1, equations) - func(u_left, 1, equations))
+    return surface_integral
+end
+
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{1}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}

src/callbacks_step/analysis_dg2d.jl

@@ -185,6 +185,51 @@ function calc_error_norms(func, u, t, analyzer,
     return l2_error, linf_error
 end
 
+# integrate `du_local` against the entropy variables evaluated at `u`
+function integrate_against_entropy_variables(du_local, u, element,
+                                             mesh::AbstractMesh{2},
+                                             equations, dg::DGSEM, cache, args...)
+    (; weights) = dg.basis
+
+    jacobian_1d = inv(cache.elements.inverse_jacobian[element]) # O(h) 
+
+    integral = zero(eltype(u))
+
+    for j in eachnode(dg), i in eachnode(dg)
+        du_node = get_node_vars(du_local, equations, dg, i, j)
+        u_node = get_node_vars(u, equations, dg, i, j, element)
+        integral = integral +
+                   weights[i] * weights[j] *
+                   dot(cons2entropy(u_node, equations), du_node)
+    end
+
+    return integral * jacobian_1d
+end
+
+# calculate surface integral of func(u, equations) * normal
+function surface_integral(func::Func, u, element, mesh::TreeMesh{2}, equations,
+                          dg::DGSEM, cache) where {Func}
+    surface_integral = zero(real(dg))
+    for ii in eachnode(dg)
+        # x direction
+        u_left = get_node_vars(u, equations, dg, 1, ii, element)
+        u_right = get_node_vars(u, equations, dg, nnodes(dg), ii, element)
+        surface_integral = surface_integral +
+                           dg.basis.weights[ii] *
+                           (func(u_right, 1, equations) - func(u_left, 1, equations))
+
+        # y direction
+        u_left = get_node_vars(u, equations, dg, ii, 1, element)
+        u_right = get_node_vars(u, equations, dg, ii, nnodes(dg), element)
+        surface_integral = surface_integral +
+                           dg.basis.weights[ii] *
+                           (func(u_right, 2, equations) - func(u_left, 2, equations))
+    end
+
+    jacobian_1d = inv(cache.elements.inverse_jacobian[element]) # O(h) 
+    return surface_integral * jacobian_1d
+end
+
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}

src/equations/compressible_euler_1d.jl

@@ -1147,4 +1147,19 @@ of total energy and kinetic energy.
 @inline function energy_internal(cons, equations::CompressibleEulerEquations1D)
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
+
+@doc raw"""
+    entropy_potential(u, orientation::Int, equations::CompressibleEulerEquations1D)
+
+Calculate the entropy potential, which for the compressible Euler equations is simply 
+momentum for the choice of entropy ``S(u) = \frac{\rho s}{\gamma - 1}``. 
+
+"Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules 
+for hyperbolic conservation laws" by Chen and Shu (2017). 
+<http://dx.doi.org/10.1016/j.jcp.2017.05.025>    
+"""
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::CompressibleEulerEquations1D)
+    return u[2]
+end
 end # @muladd

src/equations/compressible_euler_2d.jl

@@ -2230,6 +2230,25 @@ end
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
 
+@doc raw"""
+    entropy_potential(u, orientation::Int, equations::CompressibleEulerEquations1D)
+
+Calculate the entropy potential, which for the compressible Euler equations is simply 
+the momentum for the choice of entropy ``S(u) = \frac{\rho s}{\gamma - 1}``. 
+
+"Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules 
+for hyperbolic conservation laws" by Chen and Shu (2017). 
+<http://dx.doi.org/10.1016/j.jcp.2017.05.025>    
+"""
+@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/equations/nonideal_compressible_euler_1d.jl

@@ -351,4 +351,18 @@ end
     rho_e = rho_e_total - 0.5f0 * rho_v1^2 / rho
     return rho_e
 end
+
+"""
+    entropy_potential(u, orientation::Int, equations::CompressibleEulerEquations1D)
+
+Calculate the entropy potential, which for the compressible Euler equations with general
+EOS is ``p v_i / T`` for the choice of entropy ``S(u) = \rho s``. 
+"""
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::NonIdealCompressibleEulerEquations1D)
+    eos = equations.equation_of_state
+    V, v1, T = cons2prim(u, equations)
+    p = pressure(V, T, eos)
+    return p * v1 / T
+end
 end # @muladd

src/solvers/dg.jl

@@ -280,6 +280,69 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralEntropyCorrection(equations, basis;
+                                    volume_flux_dg=flux_central,
+                                    volume_flux_fv=flux_lax_friedrichs, 
+                                    scaling = true)
+
+Entropy correction volume integral type for DG methods using a convex blending of
+the **first-order** finite volume method with numerical flux `volume_flux_fv` and the
+[`VolumeIntegralFluxDifferencing`](@ref) with volume flux `volume_flux_dg`.
+The amount of blending is determined by the violation of a cell entropy equality by
+the volume integral. 
+
+`scaling ≥ 1` arbitrarily scales the blending parameter by a constant, increasing 
+the amount of the subcell FV added in. This can be used to add shock capturing-like 
+behavior. 
+
+The use of `VolumeIntegralEntropyCorrection` requires either
+    `entropy_potential(u, orientation, equations)` for TreeMesh, or
+    `entropy_potential(u, normal_direction, equations)` for other mesh types
+to be defined. 
+"""
+struct VolumeIntegralEntropyCorrection{VolumeFluxDG, VolumeFluxFV, Indicator} <:
+       AbstractVolumeIntegralShockCapturing
+    volume_flux_dg::VolumeFluxDG # symmetric, e.g. split-form or entropy-conservative
+    volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
+    indicator::Indicator
+end
+
+function VolumeIntegralEntropyCorrection(equations, basis;
+                                         volume_flux_dg = flux_central,
+                                         volume_flux_fv = flux_lax_friedrichs,
+                                         scaling = true)
+    indicator = IndicatorEntropyCorrection(equations, basis; scaling)
+    return VolumeIntegralEntropyCorrection{typeof(volume_flux_dg),
+                                           typeof(volume_flux_fv),
+                                           typeof(indicator)}(volume_flux_dg,
+                                                              volume_flux_fv,
+                                                              indicator)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain",
+                   integral::VolumeIntegralEntropyCorrection)
+    @nospecialize integral # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        summary_header(io, "VolumeIntegralEntropyCorrection")
+        summary_line(io, "volume flux DG", integral.volume_flux_dg)
+        summary_line(io, "volume flux FV", integral.volume_flux_fv)
+        summary_line(io, "indicator", integral.indicator |> typeof |> nameof)
+        show(increment_indent(io), mime, integral.indicator)
+        summary_footer(io)
+    end
+end
+
+function get_element_variables!(element_variables, u, mesh, equations,
+                                volume_integral::VolumeIntegralEntropyCorrection,
+                                dg, cache)
+    element_variables[:indicator_shock_capturing] = volume_integral.indicator.cache.alpha
+    return nothing
+end
+
 # Abstract supertype for first-order `VolumeIntegralPureLGLFiniteVolume` and
 # second-order `VolumeIntegralPureLGLFiniteVolumeO2` subcell-based finite volume
 # volume integrals.

src/solvers/dgsem/calc_volume_integral.jl

@@ -155,4 +155,74 @@ function calc_volume_integral!(du, u, mesh,
 
     return nothing
 end
+
+@inline regularized_ratio(a, b) = a * b / (eps(b) + b^2)
+
+function calc_volume_integral!(du, u, mesh,
+                               have_nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralEntropyCorrection,
+                               dg::DGSEM, cache)
+    (; volume_flux_dg, volume_flux_fv, indicator) = volume_integral
+    du_element_threaded = indicator.cache.indicator_threaded
+    (; scaling) = indicator
+    (; alpha) = indicator.cache
+    resize!(alpha, nelements(dg, cache))
+
+    @threaded for element in eachelement(dg, cache)
+        flux_differencing_kernel!(du, u, element, mesh,
+                                  have_nonconservative_terms, equations,
+                                  volume_flux_dg, dg, cache)
+
+        # check entropy production of "high order" volume integral         
+        volume_integral_entropy_vars = integrate_against_entropy_variables(view(du, ..,
+                                                                                element),
+                                                                           u, element,
+                                                                           mesh,
+                                                                           equations,
+                                                                           dg, cache)
+        surface_integral_entropy_potential = surface_integral(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 = -(volume_integral_entropy_vars +
+                             surface_integral_entropy_potential)
+
+        if entropy_residual > 0
+            # Store "high order" result
+            du_element = du_element_threaded[Threads.threadid()]
+            @views du_element .= du[.., element]
+
+            # Reset pure flux-differencing volume integral 
+            du[.., element] .= zero(eltype(du))
+
+            # Calculate FV volume integral contribution
+            fv_kernel!(du, u, mesh,
+                       have_nonconservative_terms, equations,
+                       volume_flux_fv, dg, cache, element)
+
+            # calculate difference between high and low order FV integral;
+            # this should be entropy dissipative if entropy_residual > 0.
+            @views du_element .= (du_element .- du[.., element])
+
+            entropy_dissipation = integrate_against_entropy_variables(du_element, u,
+                                                                      element,
+                                                                      mesh, equations,
+                                                                      dg, cache)
+
+            # calculate blending factor 
+            ratio = regularized_ratio(-entropy_residual, entropy_dissipation)
+            theta = min(1, scaling * ratio) # TODO: replacing this with a differentiable version of `min`
+
+            # save blending coefficient for visualization
+            alpha[element] = theta
+
+            # blend the high order method back in 
+            @views du[.., element] .= du[.., element] .+ (1 - theta) .* du_element
+        end
+    end
+
+    return nothing
+end
 end # @muladd