From cdd552ff5ee342e8e9343cdc10a75aebb10df73f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 14 Oct 2024 11:55:47 +0200 Subject: [PATCH] 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 <74359358+bennibolm@users.noreply.github.com> * 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 <74359358+bennibolm@users.noreply.github.com> --- README.md | 4 +- docs/literate/src/files/shock_capturing.jl | 141 +++++++++++++++++- docs/src/index.md | 4 + .../elixir_mhd_amr_entropy_bounded.jl | 128 ++++++++++++++++ ...elixir_euler_blast_wave_entropy_bounded.jl | 80 ++++++++++ ...uler_colliding_flow_amr_entropy_bounded.jl | 109 ++++++++++++++ src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + .../entropy_bounded_limiter.jl | 70 +++++++++ .../entropy_bounded_limiter_1d.jl | 72 +++++++++ .../entropy_bounded_limiter_2d.jl | 76 ++++++++++ .../entropy_bounded_limiter_3d.jl | 76 ++++++++++ test/test_p4est_3d.jl | 37 +++++ test/test_tree_1d_euler.jl | 15 ++ test/test_tree_2d_euler.jl | 28 ++++ 15 files changed, 836 insertions(+), 7 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl create mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_1d.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_2d.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_3d.jl diff --git a/README.md b/README.md index 30e9da243d..70c28799f3 100644 --- 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 index dd6698c2a8..62388416d8 100644 --- 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 index 0e4749dde3..88752cc7f1 100644 --- 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 new file mode 100644 index 0000000000..00ab84bf4c --- /dev/null +++ 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 new file mode 100644 index 0000000000..fa81defb6c --- /dev/null +++ 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 diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl new file mode 100644 index 0000000000..6e7cac9a17 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.001 # almost isothermal when gamma reaches 1 +equations = CompressibleEulerEquations2D(gamma) + +# This is a hand made colliding flow setup without reference. Features Mach=70 inflow from both +# sides, with relative low temperature, such that pressure keeps relatively small +# Computed with gamma close to 1, to simulate isothermal gas +function initial_condition_colliding_flow_astro(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF) + # domain size is [-64,+64]^2 + @unpack gamma = equations + # the quantities are chosen such, that they are as close as possible to the astro examples + # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...) + rho = 0.0247 + c = 0.2 + p = c^2 / gamma * rho + vel = 13.907432274789372 + slope = 1.0 + v1 = -vel * tanh(slope * x[1]) + # add small initial disturbance to the field, but only close to the interface + if abs(x[1]) < 10 + v1 = v1 * (1 + 0.01 * sin(pi * x[2])) + end + v2 = 0.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_colliding_flow_astro + +boundary_conditions = (x_neg = BoundaryConditionDirichlet(initial_condition_colliding_flow_astro), + x_pos = BoundaryConditionDirichlet(initial_condition_colliding_flow_astro), + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +surface_flux = flux_lax_friedrichs # HLLC needs more shock capturing (alpha_max) +volume_flux = flux_ranocha # works with Chandrashekar flux as well +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# shock capturing necessary for this tough example, however alpha_max = 0.5 is fine +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.0001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-64.0, -64.0) +coordinates_max = (64.0, 64.0) + +# only refinement in a patch. Needs x=-17/+17 to trigger refinement due to coarse base mesh +refinement_patches = ((type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)) + #(type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)), # very high resolution, takes about 1000s on 2 cores + ) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + refinement_patches = refinement_patches, + periodicity = (false, true), + n_cells_max = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 25.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + +stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -1.3e-4) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + dt = 1e-2, save_everystep = false, adaptive = true, + callback = callbacks); + +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index d9b9245918..0cedec782d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -274,7 +274,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, export ControllerThreeLevel, ControllerThreeLevelCombined, IndicatorLöhner, IndicatorLoehner, IndicatorMax -export PositivityPreservingLimiterZhangShu +export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured, ode_default_options diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index d5abc1d227..f4f6e72e84 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -8,4 +8,5 @@ include("positivity_zhang_shu.jl") include("subcell_limiter_idp_correction.jl") include("subcell_bounds_check.jl") +include("entropy_bounded_limiter.jl") end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter.jl b/src/callbacks_stage/entropy_bounded_limiter.jl new file mode 100644 index 0000000000..f8dc9a612d --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter.jl @@ -0,0 +1,70 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +mutable struct EntropyBoundedLimiter{RealT <: Real} + exp_entropy_decrease_max::RealT # < 0 +end + +""" + EntropyBoundedLimiter{RealT <: Real}(; exp_entropy_decrease_max = 1f-13) + +Entropy-bounded limiter by +- Lv, Ihme (2015) + Entropy-bounded discontinuous Galerkin scheme for Euler equations + [doi: 10.1016/j.jcp.2015.04.026](https://doi.org/10.1016/j.jcp.2015.04.026) + +This is an ideal-gas specific limiter that bounds the (unphysical) decrease of the thermodynamic entropy +per element from one time step (or Runge-Kutta stage) to the next. +The parameter `exp_entropy_decrease_max` is the maximum allowed exponentiated +entropy decrease per element at each element's node. + +In the original version of the paper, this value is set to zero to +ensure that the entropy does not decrease, i.e., guarantee entropy stability in the sense of +- Tadmor (1986) + A minimum entropy principle in the gas dynamics equations + [doi: 10.1016/0168-9274(86)90029-2](https://doi.org/10.1016/0168-9274(86)90029-2) + +This, however, leads in general to very diffusive solutions for timesteps violating +a CFL condition (Lemma 3 in Lv, Ihme (2015)) which is required for entropy stability in the mean values. +Since most practical simulations will employ a significantly larger timestep, one can relax the +strict entropy increase requirement by setting `exp_entropy_decrease_max` to a negative value. +The limiter acts if the exponentiated entropy decrease on an element is larger than +`exp_entropy_decrease_max`. +This means that if the change in exponentiated entropy lies *below* `exp_entropy_decrease_max` (i.e., larger in absolute value) the limiter takes action. +The choice of the tolerated exponentiated entropy decrease is a problem-specific parameter +which balances the trade-off between accuracy and stability. +""" +function EntropyBoundedLimiter(; + exp_entropy_decrease_max::RealT = -1.0f-13) where {RealT <: + Real} + @assert exp_entropy_decrease_max<0 "Supplied `exp_entropy_decrease_max` expected to be negative" + EntropyBoundedLimiter{RealT}(exp_entropy_decrease_max) +end + +function (limiter!::EntropyBoundedLimiter)(u_ode, integrator, + semi::AbstractSemidiscretization, + t) + @trixi_timeit timer() "entropy-bounded limiter" begin + @assert :uprev in fieldnames(typeof(integrator)) "EntropyBoundedLimiter requires `uprev` for computation of previous entropy" + + u = wrap_array(u_ode, semi) + u_prev = wrap_array(integrator.uprev, semi) + limiter_entropy_bounded!(u, u_prev, limiter!.exp_entropy_decrease_max, + mesh_equations_solver_cache(semi)...) + end +end + +# Exponentiated entropy change for the thermodynamic entropy (see `entropy_thermodynamic`) +# of an ideal gas with constant gamma. +@inline function exp_entropy_change(p, rho, gamma, exp_entropy_prev) + return p - rho^gamma * exp_entropy_prev +end + +include("entropy_bounded_limiter_1d.jl") +include("entropy_bounded_limiter_2d.jl") +include("entropy_bounded_limiter_3d.jl") +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_1d.jl b/src/callbacks_stage/entropy_bounded_limiter_1d.jl new file mode 100644 index 0000000000..217261d39f --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_1d.jl @@ -0,0 +1,72 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{1}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. + # Limiting only if entropy DECREASE below a user defined threshold is detected. + d_exp_s_min < exp_entropy_decrease_max || continue + + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, element)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # Note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_2d.jl b/src/callbacks_stage/entropy_bounded_limiter_2d.jl new file mode 100644 index 0000000000..36d70e4f2e --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_2d.jl @@ -0,0 +1,76 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{2}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, j, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. + # Limiting only if entropy DECREASE below a user defined threshold is detected. + d_exp_s_min < exp_entropy_decrease_max || continue + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element)) + total_volume = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, element))) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, j, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_3d.jl b/src/callbacks_stage/entropy_bounded_limiter_3d.jl new file mode 100644 index 0000000000..6f194f0814 --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_3d.jl @@ -0,0 +1,76 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{3}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, j, k, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. Avoid division by ("near") zero + d_exp_s_min < exp_entropy_decrease_max || continue + + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element)) + total_volume = zero(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, k, element))) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, j, k, element) + end + end + + return nothing +end +end # @muladd diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 55f7ab7602..5cec16300f 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -594,6 +594,43 @@ end end end +@trixi_testset "elixir_mhd_amr_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_amr_entropy_bounded.jl"), + l2=[ + 0.005430176785094096, + 0.006185803468926062, + 0.012158513265762224, + 0.006185144232789619, + 0.03509140423905665, + 0.004968215426326584, + 0.006553519141867704, + 0.005008885124643863, + 5.165777182726578e-6 + ], + linf=[ + 0.1864317840224794, + 0.2041246899193812, + 0.36992946717578445, + 0.2327158690965257, + 1.0368624176126007, + 0.1846308291826353, + 0.2062255411778191, + 0.18955666546331185, + 0.0005208969502913304 + ], + tspan=(0.0, 0.04), + coverage_override=(maxiters = 6, initial_refinement_level = 1, + base_level = 1, max_level = 2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_linearizedeuler_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_convergence.jl"), diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 7490855252..1386c5f5c4 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -409,6 +409,21 @@ end end end +@trixi_testset "elixir_euler_blast_wave_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_entropy_bounded.jl"), + l2=[0.9689207881108007, 0.1617708899929322, 1.3847895715669456], + linf=[2.95591859210077, 0.3135723412205586, 2.3871554358655365]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "test_quasi_1D_entropy" begin a = 0.9 u_1D = SVector(1.1, 0.2, 2.1) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 6685c289a9..479a07b97a 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -732,6 +732,34 @@ end end end +@trixi_testset "elixir_euler_colliding_flow_amr_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_colliding_flow_amr_entropy_bounded.jl"), + l2=[ + 0.04120588220419942, + 0.09868046588789257, + 7.446553779796626e-7, + 5.5746513268066105 + ], + linf=[ + 0.3478655090378702, + 0.864011305195807, + 5.419432288048388e-5, + 47.284459667934684 + ], + tspan=(0.0, 1.0), + dt=2.5e-2, adaptive=false, + coverage_override=(maxiters = 10^3,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_astro_jet_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_astro_jet_amr.jl"), l2=[