Add tutorial for subcell IDP limiting

docs/literate/src/files/index.jl

@@ -51,71 +51,75 @@
 # explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler
 # equations.
 
-# ### [6 Non-periodic boundary conditions](@ref non_periodic_boundaries)
+# ### [6 Subcell limiting with the IDP Limiter](@ref subcell_shock_capturing)
+#-
+# TODO
+
+# ### [7 Non-periodic boundary conditions](@ref non_periodic_boundaries)
 #-
 # Thus far, all examples used periodic boundaries. In Trixi.jl, you can also set up a simulation with
 # non-periodic boundaries. This tutorial presents the implementation of the classical Dirichlet
 # boundary condition with a following example. Then, other non-periodic boundaries are mentioned.
 
-# ### [7 DG schemes via `DGMulti` solver](@ref DGMulti_1)
+# ### [8 DG schemes via `DGMulti` solver](@ref DGMulti_1)
 #-
 # This tutorial is about the more general DG solver [`DGMulti`](@ref), introduced [here](@ref DGMulti).
 # We are showing some examples for this solver, for instance with discretization nodes by Gauss or
 # triangular elements. Moreover, we present a simple way to include pre-defined triangulate meshes for
 # non-Cartesian domains using the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl).
 
-# ### [8 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
+# ### [9 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
 #-
 # Supplementary to the previous tutorial about DG schemes via the `DGMulti` solver we now present
 # the possibility for `DGMulti` to use other SBP schemes via the package
 # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
 # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin
 # (CGSEM) method.
 
-# ### [9 Upwind FD SBP schemes](@ref upwind_fdsbp)
+# ### [10 Upwind FD SBP schemes](@ref upwind_fdsbp)
 #-
 # General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but
 # also with a general `DG` solver. In particular, upwind finite difference SBP
 # methods can be used together with the `TreeMesh`. Similar to general SBP
 # schemes in the `DGMulti` framework, the interface is based on the package
 # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
 
-# ### [10 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
+# ### [11 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
 #-
 # This tutorial explains how to add a new physics model using the example of the cubic conservation
 # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
 # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
 # and the ODE problem is solved by OrdinaryDiffEq's `solve` method.
 
-# ### [11 Adding a non-conservative equation](@ref adding_nonconservative_equation)
+# ### [12 Adding a non-conservative equation](@ref adding_nonconservative_equation)
 #-
 # In this part, another physics model is implemented, the nonconservative linear advection equation.
 # We run two different simulations with different levels of refinement and compare the resulting errors.
 
-# ### [12 Parabolic terms](@ref parabolic_terms)
+# ### [13 Parabolic terms](@ref parabolic_terms)
 #-
 # This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g.,
 # to solve the advection-diffusion equation.
 
-# ### [13 Adding new parabolic terms](@ref adding_new_parabolic_terms)
+# ### [14 Adding new parabolic terms](@ref adding_new_parabolic_terms)
 #-
 # This tutorial describes how new parabolic terms can be implemented using Trixi.jl.
 
-# ### [14 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
+# ### [15 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
 #-
 # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while
 # not wasting resources for less interesting parts of the domain. This leads to much more efficient
 # simulations. This tutorial presents the implementation strategy of AMR in Trixi.jl, including the use of
 # different indicators and controllers.
 
-# ### [15 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
+# ### [16 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
 #-
 # In this tutorial, the use of Trixi.jl's structured curved mesh type [`StructuredMesh`](@ref) is explained.
 # We present the two basic option to initialize such a mesh. First, the curved domain boundaries
 # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is
 # defined by passing the transformation mapping.
 
-# ### [16 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
+# ### [17 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
 #-
 # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
 # functionality of Trixi.jl. This begins by running and visualizing an available unstructured
@@ -124,26 +128,26 @@
 # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.
 # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`.
 
-# ### [17 P4est mesh from gmsh](@ref p4est_from_gmsh)
+# ### [18 P4est mesh from gmsh](@ref p4est_from_gmsh)
 #-
 # This tutorial describes how to obtain a [`P4estMesh`](@ref) from an existing mesh generated
 # by [`gmsh`](https://gmsh.info/) or any other meshing software that can export to the Abaqus
 # input `.inp` format. The tutorial demonstrates how edges/faces can be associated with boundary conditions based on the physical nodesets.
 
-# ### [18 Explicit time stepping](@ref time_stepping)
+# ### [19 Explicit time stepping](@ref time_stepping)
 #-
 # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
 # It explains how to use their algorithms and presents two types of time step choices - with error-based
 # and CFL-based adaptive step size control.
 
-# ### [19 Differentiable programming](@ref differentiable_programming)
+# ### [20 Differentiable programming](@ref differentiable_programming)
 #-
 # This part deals with some basic differentiable programming topics. For example, a Jacobian, its
 # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for
 # a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl
 # at the end.
 
-# ### [20 Custom semidiscretization](@ref custom_semidiscretization)
+# ### [21 Custom semidiscretization](@ref custom_semidiscretization)
 #-
 # This tutorial describes the [semidiscretiations](@ref overview-semidiscretizations) of Trixi.jl
 # and explains how to extend them for custom tasks.

docs/literate/src/files/subcell_shock_capturing.jl

@@ -0,0 +1,104 @@
+#src # Subcell limiting with the IDP Limiter
+
+# In the previous tutorial, the element-wise limiting with [`IndicatorHennemannGassner`](@ref)
+# and [`VolumeIntegralShockCapturingHG`](@ref) was explained. This tutorial contains a short
+# introduction to the idea and implementation of subcell shock capturing approaches in Trixi.jl,
+# which is also based on the DGSEM scheme in flux differencing formulation.
+# Trixi.jl contains the a-posteriori invariant domain-preserving (IDP) limiter which was
+# introduced by [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.1016/j.compfluid.2022.105627)
+# and [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876). It is a flux-corrected
+# transport-type (FCT) limiter and is implemented using [`SubcellLimiterIDP`](@ref) and
+# [`VolumeIntegralSubcellLimiting`](@ref).
+# Since it is an a-posteriori limiter you have to apply a correction stage after each Runge-Kutta
+# stage. This is done by passing the stage callback [`SubcellLimiterIDPCorrection`](@ref) to the
+# time integration method.
+
+# ## Time integration method
+# As mentioned before, the IDP limiting is an a-posteriori limiter. Its limiting process
+# guarantees the target bounds for a simple Euler evolution. To still achieve a high-order
+# approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta method
+# which can be written as convex combination of these forward Euler steps.
+#-
+# Due to this functionality of the limiting procedure the correcting stage and therefore the stage
+# callbacks has to be applied to the solution after the forward Euler step and before further
+# computation. Unfortunately, the `solve(...)` routines of
+# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), which is normally used in
+# Trixi.jl for the time integration, does not support calculations via callback at this point
+# in the simulation.
+#-
+# Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern
+# time integration SSPRK method called with
+# ````julia
+# Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)`.
+# ````
+#-
+# Right now, only the third-order SSPRK method [`Trixi.SimpleSSPRK33`](@ref) is implemented.
+
+# TODO: Some comments about
+# - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)))
+# - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds)
+
+using Trixi
+
+# # `SubcellLimiterIDP`
+# The IDP limiter supports several options of limiting which are passed very flexible as parameters to
+# the limiter individually.
+
+# ## Global bounds
+# First, there is the use of global bounds. If enabled, they enforce physical admissibility
+# conditions, such as non-negativity of variables.
+# This can be done for conservative variables, where the limiter is of a one-sided Zalesak-type, and
+# general non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds.
+
+# ### Conservative variables
+# The procedure to enforce global bounds for a conservative variables is as follows:
+# If you want to guarantee non-negativity for the density of compressible Euler equations,
+# you pass the specific quantity name of the conservative variable.
+equations = CompressibleEulerEquations2D(1.4)
+
+# The quantity name of the density is `rho` shich is how we enable its limiting.
+# ````julia
+# positivity_variables_cons = ["rho"]
+# ````
+
+# The quantity names are passed as a vector to allow several quantities.
+# This is for instance used if you want to limit the density of two different components using
+# the multicomponent compressible Euler equations.
+equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
+                                                       gas_constants = (0.287, 1.578))
+
+# Then, we just pass both quantity names.
+# ````julia
+# positivity_variables_cons = ["rho1", "rho2"]
+# ````
+
+# Alternatively, it is possible to all limit all density variables with a general command using
+# ````julia
+# positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)]
+# ````
+
+# ### Non-linear variables
+# To allow limitation for all possible non-linear variables including on-the-fly defined ones,
+# you directly pass function here.
+# For instance, if you want to enforce non-negativity for the pressure, do as follows.
+# ````julia
+# positivity_variables_nonlinear = [pressure]
+# ````
+
+# ## Local bounds (Shock capturing)
+# Second, Trixi.jl supports the limiting with local bounds for conservative variables. They
+# allow to avoid spurious  oscillations within the global bounds and to improve the
+# shock-capturing capabilities of the method. The corresponding numerical admissibility
+# conditions are frequently formulated as local maximum or minimum principles.
+# There are different option to choose local bounds. We calculate them using the low-order FV
+# solution.
+
+# As for the limiting with global bounds you are passing the quantity names of the conservative
+# variables you want to limit. So, to limit the density with lower and upper local bounds pass
+# the following.
+# ````julia
+# local_minmax_variables_cons = ["rho"]
+# ````
+
+# ## Exemplary simulation
+# TODO

docs/make.jl

@@ -60,6 +60,7 @@ files = [
     "Introduction to DG methods" => "scalar_linear_advection_1d.jl",
     "DGSEM with flux differencing" => "DGSEM_FluxDiff.jl",
     "Shock capturing with flux differencing and stage limiter" => "shock_capturing.jl",
+    "Subcell limiting with the IDP Limiter" => "subcell_shock_capturing.jl",
     "Non-periodic boundaries" => "non_periodic_boundaries.jl",
     "DG schemes via `DGMulti` solver" => "DGMulti_1.jl",
     "Other SBP schemes (FD, CGSEM) via `DGMulti` solver" => "DGMulti_2.jl",