Merge branch 'main' into t8codemesh-fv

trixi-framework · Oct 15, 2024 · a8f35ef · a8f35ef
2 parents 9c10471 + cdd552f
commit a8f35ef
Show file tree

Hide file tree

Showing 89 changed files with 1,754 additions and 550 deletions.
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.24.3
+        uses: crate-ci/typos@v1.25.0
diff --git a/NEWS.md b/NEWS.md
@@ -4,13 +4,39 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.9 from v0.8.x
+
+#### Added
+
+- Boundary conditions are now supported on nonconservative terms ([#2062]).
+
+#### Changed
+
+- We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
+  `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`.
+  ([#2069])
+- In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
+  Users should not be affected by this.
+- Nonconservative terms depend only on `normal_direction_average` instead of both 
+  `normal_direction_average` and `normal_direction_ll`, such that the function signature is now 
+  identical with conservative fluxes. This required a change of the `normal_direction` in
+  `flux_nonconservative_powell` ([#2062]).
+
+#### Deprecated
+
+#### Removed
+
+
 ## Changes in the v0.8 lifecycle
 
 #### Changed
 
 - The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product
   of the Jacobian and the conserved variables instead of the conserved variables only
   to make AMR fully conservative ([#2028]). This may change AMR results slightly.
+- Subcell (IDP) limiting is now officially supported and not marked as experimental
+  anymore (see `VolumeIntegralSubcellLimiting`).
 
 ## Changes when updating to v0.8 from v0.7.x
 

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.8.11-DEV"
+version = "0.9.2-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

diff --git a/README.md b/README.md
@@ -19,28 +19,6 @@
   <img width="300px" src="https://trixi-framework.github.io/assets/logo.png">
 </p>
 
-***
-**Trixi.jl at JuliaCon 2024**<br/>
-At this year's JuliaCon in Eindhoven, Netherlands, we will be present with several contributions
-from the Trixi Framework ecosystem:
-
-* [**Julia for Particle-Based Multiphysics with TrixiParticles.jl**](https://pretalx.com/juliacon2024/talk/TPFF8L/),<br/>
-  [*Erik Faulhaber*](https://github.com/efaulhaber/), [*Niklas Neher*](https://github.com/lasnikas/),
-  10th July 2024, 11:00–11:30, Function (4.1)
-* [**Towards Aerodynamic Simulations in Julia with Trixi.jl**](https://pretalx.com/juliacon2024/talk/XH8KBG/),<br/>
-  [*Daniel Doehring*](https://github.com/danieldoehring/),
-  10th July 2024, 15:30–16:00, While Loop (4.2)
-* [**libtrixi: serving legacy codes in earth system modeling with fresh Julia CFD**](https://pretalx.com/juliacon2024/talk/SXC7LA/),<br/>
-  [*Benedict Geihe*](https://github.com/benegee/),
-  12th July 2024, 14:00–17:00, Function (4.1)
-
-The last talk is part of the 
-[Julia for High-Performance Computing](https://juliacon.org/2024/minisymposia/hpc/)
-minisymposium, which this year is hosted by our own [*Hendrik Ranocha*](https://github.com/ranocha/).
-
-We are looking forward to seeing you there ♥️
-***
-
 **Trixi.jl** is a numerical simulation framework for conservation
 laws written in [Julia](https://julialang.org). A key objective for the
 framework is to be useful to both scientists and students. Therefore, next to
@@ -58,8 +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
-  * Positivity-preserving 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/benchmark/Project.toml b/benchmark/Project.toml
@@ -8,4 +8,4 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 BenchmarkTools = "0.5, 0.7, 1.0"
 OrdinaryDiffEq = "5.65, 6"
 PkgBenchmark = "0.2.10"
-Trixi = "0.4, 0.5, 0.6, 0.7, 0.8"
+Trixi = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
@@ -35,7 +35,7 @@ plot!(Shape([(-1.7,-11), (1.7,-11), (1.7,-17), (-1.7,-17)]), linecolor=:blue, fi
 annotate!(1.5, -14, ("Trixi.jl", 12, :blue, :center))
 
 plot!(Shape([(-1.4,-14.5), (1.4,-14.5), (1.4,-16.5), (-1.4,-16.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
-annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, initial_condition, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
+annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
 
 
 

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/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl
@@ -51,7 +51,11 @@
 # The Newton-bisection algorithm is an iterative method and requires some parameters.
 # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and
 # relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are
-# sufficient in most cases and therefore used as default. Additionally, there is the parameter
+# sufficient in most cases and therefore used as default. If the implemented bounds checking
+# functionality indicates problems with the limiting (see [below](@ref subcell_bounds_check))
+# the Newton method with the chosen parameters might not manage to converge. If so, adapting
+# the mentioned parameters helps fix that.
+# Additionally, there is the parameter
 # `gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation
 # of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`,
 # as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.]
@@ -244,7 +248,7 @@ plot(sol)
 # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116)
 
 
-# ## Bounds checking
+# ## [Bounds checking](@id subcell_bounds_check)
 # Subcell limiting is based on the fulfillment of target bounds - either global or local.
 # Although the implementation works and has been thoroughly tested, there are some cases where
 # these bounds are not met.

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -30,7 +30,12 @@ installation and postprocessing procedures. Its features include:
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
   * 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_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
@@ -85,11 +85,11 @@ analysis_interval = 2000
 l_inf = 1.0 # Length of airfoil
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
@@ -89,11 +89,11 @@ rho_inf = 1.4
 u_inf = 0.38
 l_inf = 1.0 # Diameter of circle
 
-drag_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+drag_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            DragCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+lift_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            LiftCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -120,23 +120,23 @@ summary_callback = SummaryCallback()
 analysis_interval = 2000
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        DragCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        LiftCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),