Update Documentation: Optimized Schemes, PERK Section, and PairedExplicitRK2 Tutorial

docs/Project.toml

@@ -1,6 +1,8 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
+Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
+ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HOHQMesh = "e4f4c7b8-17cb-445a-93c5-f69190ed6c8c"
@@ -16,7 +18,9 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 [compat]
 CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11"
 Changelog = "1.1"
+Convex = "0.16"
 Documenter = "1"
+ECOS = "1.1.2"
 ForwardDiff = "0.10"
 HOHQMesh = "0.1, 0.2"
 LaTeXStrings = "1.2"

docs/src/time_integration.md

@@ -1,11 +1,13 @@
 # [Time integration methods](@id time-integration)
 
+## Methods from OrdinaryDiffEq.jl
+
 Trixi.jl is compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/).
-In particular, explicit Runge-Kutta methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+In particular, [explicit Runge-Kutta methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Runge-Kutta-Methods) from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 are tested extensively.
 Interesting classes of time integration schemes are
-- [Explicit low-storage Runge-Kutta methods](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))
+- [Explicit low-storage Runge-Kutta methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Low-Storage-Methods)
+- [Strong stability preserving (SSP) methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))
 
 Some common options for `solve` from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 are the following. Further documentation can be found in the
@@ -21,7 +23,7 @@ are the following. Further documentation can be found in the
   out of memory or start to swap).
 - You can set the maximal number of time steps via `maxiters=...`.
 - SSP methods and many low-storage methods from OrdinaryDiffEq.jl support
-  `stage_limiter!`s and `step_limiter!`s, e.g., [`PositivityPreservingLimiterZhangShu`](@ref)
+  `stage_limiter!`s and `step_limiter!`s, e.g., [`PositivityPreservingLimiterZhangShu`](@ref) and [`EntropyBoundedLimiter`](@ref)
   from Trixi.jl.
 - If you start Julia with multiple threads and want to use them also in the time 
   integration method from OrdinaryDiffEq.jl, you need to pass the keyword argument
@@ -31,7 +33,7 @@ are the following. Further documentation can be found in the
   For more information on using thread-based parallelism in Trixi.jl, please refer to 
   [Shared-memory parallelization with threads](@ref).
 - If you use error-based step size control (see also the section on
-  [error-based adaptive step sizes](@ref adaptive_step_sizes) together with MPI, you need to
+  [error-based adaptive step sizes](@ref adaptive_step_sizes)) together with MPI, you need to
   pass `internalnorm=ode_norm` and you should pass `unstable_check=ode_unstable_check` to
   OrdinaryDiffEq's [`solve`](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/), which are both
   included in [`ode_default_options`](@ref).
@@ -42,3 +44,172 @@ are the following. Further documentation can be found in the
     of stages, e.g. to allow for interpolation (dense output), root-finding for continuous callbacks,
     and error-based time step control. In general, you often should not need to worry about this if you
     use Trixi.jl.
+
+## Custom Optimized Schemes
+
+### Stabilized Explicit Runge-Kutta Methods
+
+Optimized explicit schemes aim to maximize the timestep $\Delta t$ for a given simulation setup.
+Formally, this boils down to an optimization problem of the form 
+```math
+\underset{P_{p;S} \, \in \, \mathcal{P}_{p;S}}{\max} \Delta t \text{ such that } \big \vert P_{p;S}(\Delta t \lambda_m) \big \vert \leq 1, \quad  m = 1 , \dots , M \tag{1}
+```
+where $p$ denotes the order of consistency of the scheme, $S$ is the number of stage evaluations and $M$ denotes the number of eigenvalues $\lambda_m$ of the Jacobian matrix $J \coloneqq \frac{\partial \boldsymbol F}{\partial \boldsymbol U}$ of the right-hand side of the [semidiscretized PDE](https://trixi-framework.github.io/Trixi.jl/stable/overview/#overview-semidiscretizations):
+```math
+\dot{\boldsymbol U} = \boldsymbol F(\boldsymbol U) \tag{2} \: .
+```
+In particular, for $S > p$ the Runge-Kutta method includes some free coefficients which may be used to adapt the domain of absolute stability to the problem at hand.
+Since Trixi.jl [supports exact computation of the Jacobian $J$ by means of automatic differentiation](https://trixi-framework.github.io/Trixi.jl/stable/tutorials/differentiable_programming/), we have access to the Jacobian of a given simulation setup.
+For small (say, up to roughly $10^4$ DoF) systems, the spectrum $\boldsymbol \sigma = \left \{ \lambda_m \right \}_{m=1, \dots, M}$ can be computed directly using [`LinearAlgebra.eigvals(J)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigvals).
+For larger systems, we recommend the procedure outlined in section 4.1 of [Doehring et al. (2024)](https://doi.org/10.1016/j.jcp.2024.113223). This approach computes a reduced set of (estimated) eigenvalues $\widetilde{\boldsymbol \sigma}$ around the convex hull of the spectrum by means of the [Arnoldi method](https://github.com/JuliaLinearAlgebra/Arpack.jl).
+
+The optimization problem (1) can be solved using the algorithms described in [Ketcheson, Ahmadia (2012)](http://dx.doi.org/10.2140/camcos.2012.7.247) for a moderate number of stages $S$ or [Doehring, Gassner, Torrilhon (2024)](https://doi.org/10.1007/s10915-024-02478-5) for a large number of stages $S$.
+In Trixi.jl, the former approach is implemented by means of convex optimization using the [Convex.jl](https://github.com/jump-dev/Convex.jl) package.
+
+The resulting stability polynomial $P_{p;S}$ is then used to construct an optimized Runge-Kutta method.
+Trixi.jl implements the [Paired-Explicit Runge-Kutta (PERK)](https://doi.org/10.1016/j.jcp.2019.05.014) method, a low-storage, multirate-ready method with optimized domain of absolute stability.
+
+### Paired-Explicit Runge-Kutta (PERK) Schemes
+
+Paired-Explicit Runge-Kutta (PERK) or `PairedExplicitRK` schemes are an advanced class of numerical methods designed to efficiently solve ODEs.
+In the [original publication](https://doi.org/10.1016/j.jcp.2019.05.014), second-order schemes were introduced, which have been extended to [third](https://doi.org/10.1016/j.jcp.2022.111470) and [fourth](https://doi.org/10.48550/arXiv.2408.05470) order in subsequent works.
+
+By construction, PERK schemes are suited for integrating multirate systems, i.e., systems with varying characteristics speeds throughout the domain.
+Nevertheless, due to their optimized stability properties and low-storage nature, the PERK schemes are also highly efficient when applied as standalone methods. In Trixi.jl, the _standalone_ PERK integrators are implemented such that all stages of the method are active.
+
+#### Tutorial: Using `PairedExplicitRK2`
+
+In this tutorial, we will demonstrate how you can use the second-order PERK time integrator. You need the packages `Convex.jl` and `ECOS.jl`, so be sure they are added to your environment.
+
+First, you need to load the necessary packages:
+
+```@example PERK-example-1
+using Convex, ECOS
+using OrdinaryDiffEq
+using Trixi
+```
+
+Then, define the ODE problem and the semidiscretization setup. For this example, we will use a simple advection problem.
+
+```@example PERK-example-1
+# Define the mesh
+cells_per_dimension = 100
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = StructuredMesh(cells_per_dimension, 
+                      coordinates_min, coordinates_max)
+
+# Define the equation and initial condition
+advection_velocity = 1.0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+# Define the solver
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# Define the semidiscretization
+semi = SemidiscretizationHyperbolic(mesh, 
+                                    equations, initial_condition, 
+                                    solver)
+```
+
+After that, we will define the necessary [callbacks](@ref callbacks-id) for the simulation. Callbacks are used to perform actions at specific points during the integration process.
+
+```@example PERK-example-1
+# Define some standard callbacks
+summary_callback  = SummaryCallback()
+alive_callback    = AliveCallback()
+analysis_callback = AnalysisCallback(semi, interval = 200)
+# For this optimized method we can use a relatively large CFL number
+stepsize_callback = StepsizeCallback(cfl = 2.5)
+
+# Create a CallbackSet to collect all callbacks
+callbacks = CallbackSet(summary_callback, 
+                        alive_callback, 
+                        analysis_callback, 
+                        stepsize_callback)
+```
+
+Now, we define the ODE problem by specifying the time span over which the ODE will be solved. 
+The `tspan` parameter is a tuple `(t_start, t_end)` that defines the start and end times for the simulation. 
+The `semidiscretize` function is used to create the ODE problem from the simulation setup.
+
+```@example PERK-example-1
+# Define the time span
+tspan = (0.0, 1.0)
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, tspan)
+```
+
+Next, we will construct the time integrator. In order to do this, you need the following components:
+
+  - Number of stages: The number of stages $S$ in the Runge-Kutta method. 
+  In this example, we use `6` stages.
+  - Time span (`tspan`): A tuple `(t_start, t_end)` that defines the time span over which the ODE will be solved. 
+  This defines the bounds for the bisection routine for the optimal timestep $\Delta t$ used in calculating the polynomial coefficients at optimization stage. 
+  This variable is already defined in step 5.
+  - Semidiscretization (`semi`): The semidiscretization setup that includes the mesh, equations, initial condition, and solver. In this example, this variable is already defined in step 3.
+  In the background, we compute from `semi` the Jacobian $J$ evaluated at the initial condition using [`jacobian_ad_forward`](https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.jacobian_ad_forward-Tuple{Trixi.AbstractSemidiscretization}).
+  This is then followed by the computation of the spectrum $\boldsymbol \sigma(J)$ using `LinearAlgebra.eigvals`.
+  Equipped with the spectrum, the optimal stability polynomial is computed, from which the corresponding Runge-Kutta method is derived. Other constructors (if the coefficients $\boldsymbol{\alpha}$ of the stability polynomial are already available, or if a reduced spectrum $\widetilde{\boldsymbol{\sigma}}$ should be used) are discussed below.
+
+```@example PERK-example-1
+# Construct second order-explicit Runge-Kutta method with 6 stages for given simulation setup (`semi`)
+# `tspan` provides the bounds for the bisection routine that is used to calculate the maximum timestep.
+ode_algorithm = Trixi.PairedExplicitRK2(6, tspan, semi)
+```
+
+With everything set up, you can now use `Trixi.solve` to solve the ODE problem. The `solve` function takes the ODE problem, the time integrator, and some options such as the time step (`dt`), whether to save every step (`save_everystep`), and the callbacks.
+
+```@example PERK-example-1
+# Solve the ODE problem using PERK2
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 1.0, # overwritten by `stepsize_callback`
+                  save_everystep = false, callback = callbacks)
+```
+
+##### Advanced constructors:
+There are two additional constructors for the `PairedExplicitRK2` method besides the one taking in a semidiscretization `semi`:
+  - `PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString)` constructs a `num_stages`-stage method from the given optimal monomial_coeffs $\boldsymbol \alpha$.
+  These are expected to be present in the provided directory in the form of a `gamma_<S>.txt` file, where `<S>` is the number of stages `num_stages`.
+  This constructor is useful when the optimal coefficients cannot be obtained using the optimization routine by Ketcheson and Ahmadia, possibly due to a large number of stages $S$.
+  - `PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64})` constructs a `num_stages`-stage using the optimization approach by Ketcheson and Ahmadia for the (reduced) spectrum `eig_vals`.
+  The use-case for this constructor would be a large system, for which the computation of all eigenvalues is infeasible.
+
+#### Automatic computation of a stable CFL Number
+
+In the previous tutorial the CFL number was set manually to $2.5$.
+To avoid the manual trial-and-error process behind this, instantiations of `AbstractPairedExplicitRK` methods can automatically compute the stable CFL number for a given simulation setup using the [`Trixi.calculate_cfl`](@ref) function.
+When constructing the time integrator from a semidiscretization `semi`, 
+```@example PERK-example-1
+# Construct third-order paired-explicit Runge-Kutta method with 8 stages for given simulation setup.                               
+ode_algorithm = Trixi.PairedExplicitRK3(8, tspan, semi)
+```
+the maximum timestep `dt` is stored by the `ode_algorithm`.
+This can then be used to compute the stable CFL number for the given simulation setup:
+```@example PERK-example-1
+cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
+```
+For nonlinear problems, the spectrum will in general change over the course of the simulation.
+Thus, it is often necessary to reduce the optimal `cfl_number` by a safety factor:
+```@example PERK-example-1
+# For non-linear problems, the CFL number should be reduced by a safety factor
+# as the spectrum changes (in general) over the course of a simulation
+stepsize_callback = StepsizeCallback(cfl = 0.85 * cfl_number)
+```
+If the optimal monomial coefficients are precomputed, the user needs to provide the obtained maximum timestep `dt_opt` from the optimization at construction stage.
+The corresponding constructor has signature
+```julia
+PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString,
+                  dt_opt = nothing; cS2 = 1.0f0)
+```
+Then, the stable CFL number can be computed as described above.
+
+#### Currently implemented PERK methods
+
+##### Single/Standalone methods
+
+- [`Trixi.PairedExplicitRK2`](@ref): Second-order PERK method with at least two stages.
+- [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages.