Entropy bounded limiter (#2078)

* start working * get EB limiter to work * Refine EB limiter * fmt * 2D/3D Limiters * introduce parameter of tolerated entropy decrease * change to exp_entropy_diff * fmt * export * docstring * 3d example * exam,ples * examples * examples * fmt * adaptive * examples * fmt * tests + docu * fmt + spelling * comment * Apply suggestions from code review Co-authored-by: Benjamin Bolm <[email protected]> * Apply suggestions from code review * increase => change * fmt * spelling * revisit * revert * back to decrease * fmt * typo * comments * theta * longer sim * removed unused code * for coverage * Apply suggestions from code review * longer runtime for coverage * vals * reproducible test values * shorten * changes * testvals * test vals * Update src/callbacks_stage/entropy_bounded_limiter_2d.jl * Apply suggestions from code review * Update docs/literate/src/files/shock_capturing.jl * more iters * Apply suggestions from code review * comments * mention step_limiter! capability * Update docs/literate/src/files/shock_capturing.jl * Update src/callbacks_stage/entropy_bounded_limiter.jl * Update src/callbacks_stage/entropy_bounded_limiter.jl * Apply suggestions from code review --------- Co-authored-by: Benjamin Bolm <[email protected]>
trixi-framework · Oct 14, 2024 · cdd552f · cdd552f
1 parent a05a68d
commit cdd552f
Show file tree

Hide file tree

Showing 15 changed files with 836 additions and 7 deletions.
diff --git a/README.md b/README.md
@@ -36,9 +36,11 @@ installation and postprocessing procedures. Its features include:
 * Discontinuous Galerkin methods
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
+  * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl)
+* Advanced limiting strategies
   * Positivity-preserving limiting
   * Subcell invariant domain-preserving (IDP) limiting
-  * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl)
+  * Entropy-bounded limiting
 * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/)
   * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)
   * [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))

diff --git a/docs/literate/src/files/shock_capturing.jl b/docs/literate/src/files/shock_capturing.jl
@@ -96,19 +96,19 @@
 # [Zhang, Shu (2011)](https://doi.org/10.1098/rspa.2011.0153).
 
 # It works the following way. For every passed (scalar) variable and for every DG element we calculate
-# the minimal value $value_{min}$. If this value falls below the given threshold $\varepsilon$,
+# the minimal value $value_\text{min}$. If this value falls below the given threshold $\varepsilon$,
 # the approximation is slightly adapted such that the minimal value of the relevant variable lies
 # now above the threshold.
 # ```math
-# \underline{u}^{new} = \theta * \underline{u} + (1-\theta) * u_{mean}
+# \underline{u}^\text{new} = \theta * \underline{u} + (1-\theta) * u_\text{mean}
 # ```
 # where $\underline{u}$ are the collected pointwise evaluation coefficients in element $e$ and
-# $u_{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination
+# $u_\text{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination
 # of these two values with factor
 # ```math
-# \theta = \frac{value_{mean} - \varepsilon}{value_{mean} - value_{min}},
+# \theta = \frac{value_\text{mean} - \varepsilon}{value_\text{mean} - value_\text{min}},
 # ```
-# where $value_{mean}$ is the relevant variable evaluated for the mean value $u_{mean}$.
+# where $value_\text{mean}$ is the relevant variable evaluated for the mean value $u_\text{mean}$.
 
 # The adapted approximation keeps the exact same mean value, but the relevant variable is now greater
 # or equal the threshold $\varepsilon$ at every node in every element.
@@ -225,6 +225,137 @@ sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false
 using Plots
 plot(sol)
 
+# # Entropy bounded limiter
+
+# As argued in the description of the positivity preserving limiter above it might sometimes be
+# necessary to apply advanced techniques to ensure a physically meaningful solution.
+# Apart from the positivity of pressure and density, the physical entropy of the system should increase 
+# over the course of a simulation, see e.g. [this](https://doi.org/10.1016/0168-9274(86)90029-2) paper by Tadmor where this property is 
+# shown for the compressible Euler equations.
+# As this is not necessarily the case for the numerical approximation (especially for the high-order, non-diffusive DG discretizations), 
+# Lv and Ihme devised an a-posteriori limiter in [this paper](https://doi.org/10.1016/j.jcp.2015.04.026) which can be applied after each Runge-Kutta stage.
+# This limiter enforces a non-decrease in the physical, thermodynamic entropy $S$ 
+# by bounding the entropy decrease (entropy increase is always tolerated) $\Delta S$ in each grid cell.
+# 
+# This translates into a requirement that the entropy of the limited approximation $S\Big(\mathcal{L}\big[\boldsymbol u(\Delta t) \big] \Big)$ should be 
+# greater or equal than the previous iterates' entropy $S\big(\boldsymbol u(0) \big)$, enforced at each quadrature point:
+# ```math
+# S\Big(\mathcal{L}\big[\boldsymbol u(\Delta t, \boldsymbol{x}_i) \big] \Big) \overset{!}{\geq} S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big), \quad i = 1, \dots, (k+1)^d
+# ```
+# where $k$ denotes the polynomial degree of the element-wise solution and $d$ is the spatial dimension.
+# For an ideal gas, the thermodynamic entropy $S(\boldsymbol u) = S(p, \rho)$ is given by
+# ```math
+# S(p, \rho) = \ln \left( \frac{p}{\rho^\gamma} \right) \: .
+# ```
+# Thus, the non-decrease in entropy can be reformulated as
+# ```math
+# p(\boldsymbol{x}_i) - e^{ S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big)} \cdot \rho(\boldsymbol{x}_i)^\gamma \overset{!}{\geq} 0, \quad i = 1, \dots, (k+1)^d \: .
+# ```
+# In a practical simulation, we might tolerate a maximum (exponentiated) entropy decrease per element, i.e., 
+# ```math
+# \Delta e^S \coloneqq \min_{i} \left\{ p(\boldsymbol{x}_i) - e^{ S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big)} \cdot \rho(\boldsymbol{x}_i)^\gamma \right\} < c
+# ```
+# with hyper-parameter $c$ which is to be specified by the user.
+# The default value for the corresponding parameter $c=$ `exp_entropy_decrease_max` is set to $-10^{-13}$, i.e., slightly less than zero to 
+# avoid spurious limiter actions for cells in which the entropy remains effectively constant.
+# Other values can be specified by setting the `exp_entropy_decrease_max` keyword in the constructor of the limiter:
+# ```julia
+# stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max=-1e-9)
+# ```
+# Smaller values (larger in absolute value) for `exp_entropy_decrease_max` relax the entropy increase requirement and are thus less diffusive.
+# On the other hand, for larger values (smaller in absolute value) of `exp_entropy_decrease_max` the limiter acts more often and the solution becomes more diffusive.
+#
+# In particular, we compute again a limiting parameter $\vartheta \in [0, 1]$ which is then used to blend the 
+# unlimited nodal values $\boldsymbol u$ with the mean value $\boldsymbol u_{\text{mean}}$ of the element:
+# ```math
+# \mathcal{L} [\boldsymbol u](\vartheta) \coloneqq (1 - \vartheta) \boldsymbol u + \vartheta \cdot \boldsymbol u_{\text{mean}}
+# ```
+# For the exact definition of $\vartheta$ the interested user is referred to section 4.4 of the paper by Lv and Ihme.
+# Note that therein the limiting parameter is denoted by $\epsilon$, which is not to be confused with the threshold $\varepsilon$ of the Zhang-Shu limiter.
+
+# As for the positivity preserving limiter, the entropy bounded limiter may be applied after every Runge-Kutta stage.
+# Both fixed timestep methods such as [`CarpenterKennedy2N54`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) and 
+# adaptive timestep methods such as [`SSPRK43`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) are supported.
+# We would like to remark that of course every `stage_limiter!` can also be used as a `step_limiter!`, i.e., 
+# acting only after the full time step has been taken.
+
+# As an example, we consider a variant of the [1D medium blast wave example](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl)
+# wherein the shock capturing method discussed above is employed to handle the shock.
+# Here, we use a standard DG solver with HLLC surface flux:
+using Trixi
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hllc)
+
+# The remaining setup is the same as in the standard example:
+
+using OrdinaryDiffEq
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    ## Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
+    ## Set up polar coordinates
+    inicenter = SVector(0.0)
+    x_norm = x[1] - inicenter[1]
+    r = abs(x_norm)
+    ## The following code is equivalent to
+    ## phi = atan(0.0, x_norm)
+    ## cos_phi = cos(phi)
+    ## in 1D but faster
+    cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)
+
+    ## Calculate primitive variables
+    rho = r > 0.5 ? 1.0 : 1.1691
+    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
+    p = r > 0.5 ? 1.0E-3 : 1.245
+
+    return prim2cons(SVector(rho, v1, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+coordinates_min = (-2.0,)
+coordinates_max = (2.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+# We specify the `stage_limiter!` supplied to the classic SSPRK33 integrator
+# with strict entropy increase enforcement
+# (recall that the tolerated exponentiated entropy decrease is set to -1e-13).
+stage_limiter! = EntropyBoundedLimiter()
+
+# We run the simulation with the SSPRK33 method and the entropy bounded limiter:
+sol = solve(ode, SSPRK33(stage_limiter!);
+            dt = 1.0,
+            callback = callbacks);
+
+using Plots
+plot(sol)
 
 # ## Package versions
 

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -32,6 +32,10 @@ installation and postprocessing procedures. Its features include:
   * Positivity-preserving limiting
   * Subcell invariant domain-preserving (IDP) limiting
   * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl)
+* Advanced limiting strategies
+  * Positivity-preserving limiting
+  * Subcell invariant domain-preserving (IDP) limiting
+  * Entropy-bounded limiting
 * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/)
   * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)
   * [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))

diff --git a/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl b/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl
@@ -0,0 +1,128 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations3D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+
+Weak magnetic blast wave setup taken from Section 6.1 of the paper:
+- A. M. Rueda-Ramírez, S. Hennemann, F. J. Hindenlang, A. R. Winters, G. J. Gassner (2021)
+  An entropy stable nodal discontinuous Galerkin method for the resistive MHD
+  equations. Part II: Subcell finite volume shock capturing
+  [doi: 10.1016/j.jcp.2021.110580](https://doi.org/10.1016/j.jcp.2021.110580)
+"""
+function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+    # Center of the blast wave is selected for the domain [0, 3]^3
+    inicenter = (1.5, 1.5, 1.5)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    delta_0 = 0.1
+    r_0 = 0.3
+    lambda = exp(5.0 / delta_0 * (r - r_0))
+
+    prim_inner = SVector(1.2, 0.1, 0.0, 0.1, 0.9, 1.0, 1.0, 1.0, 0.0)
+    prim_outer = SVector(1.2, 0.2, -0.4, 0.2, 0.3, 1.0, 1.0, 1.0, 0.0)
+    prim_vars = (prim_inner + lambda * prim_outer) / (1.0 + lambda)
+
+    return prim2cons(prim_vars, equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+polydeg = 3
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but with slightly less warping.
+# The mapping will be interpolated at tree level, and then refined without changing
+# the geometry interpolant.
+function mapping(xi_, eta_, zeta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+    zeta = 1.5 * zeta_ + 1.5
+
+    y = eta +
+        3 / 11 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+         cos(0.5 * pi * (2 * eta - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    x = xi +
+        3 / 11 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+         cos(2 * pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    z = zeta +
+        3 / 11 * (cos(0.5 * pi * (2 * x - 3) / 3) *
+         cos(pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (2, 2, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3,
+                 mapping = mapping,
+                 initial_refinement_level = 2,
+                 periodicity = true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi,
+                                variable = density_pressure)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 2,
+                                      max_level = 4, max_threshold = 0.15)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+cfl = 3.0
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+# Stage limiter required for this high CFL                        
+stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -5e-3)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+
+A medium blast wave taken from
+- Sebastian Hennemann, Gregor J. Gassner (2020)
+  A provably entropy stable subcell shock capturing approach for high order split form DG
+  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+"""
+function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
+    # Set up polar coordinates
+    inicenter = SVector(0.0)
+    x_norm = x[1] - inicenter[1]
+    r = abs(x_norm)
+    # The following code is equivalent to
+    # phi = atan(0.0, x_norm)
+    # cos_phi = cos(phi)
+    # in 1D but faster
+    cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)
+
+    # Calculate primitive variables
+    rho = r > 0.5 ? 1.0 : 1.1691
+    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
+    p = r > 0.5 ? 1.0E-3 : 1.245
+
+    return prim2cons(SVector(rho, v1, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+# Note: We do not need to use the shock-capturing methodology here,
+# in contrast to the standard `euler_blast_wave.jl` example.
+solver = DGSEM(polydeg = 3, surface_flux = flux_hllc)
+
+coordinates_min = (-2.0,)
+coordinates_max = (2.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+stage_limiter! = EntropyBoundedLimiter()
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK33(stage_limiter!);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            callback = callbacks);
+
+summary_callback() # print the timer summary