R-RK: Newton and sub-diagonal schemes only (#2480)

* start working on relaxation RK

* start doc

* docstrings + clean up

* spell out

* docstrings, remove secant

* examples

* tests

* typo

* increase allowed allocs

* typo

* test printing/show routines

* rename to `u_tmp`

* comment

* reduce elixirs and tests

* rename more explicitly

* fmt

* document

* add sciml base

* comments and format

* more tests

* docs

* fix docs

* refactor simple ssp

* fix

* test resize

* remove unused resize!

* Apply suggestions from code review

* Apply suggestions from code review

* Update src/time_integration/relaxation_methods/entropy_relaxation.jl

* Apply suggestions from code review

* Apply suggestions from code review

* shorten

* test

* test vals

* fix tet

* R-RK: Newton and Sub-Diagonal

* fmt

* docs

* fix

* abstract type

* remove

* restore

* fix options

* 3d example

* fd

* take errors into account

* fix test

* fmt

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update docs/src/time_integration.md

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* remove sciml

* rename to ralston

* rename

* remove 2d

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: DanielDoehring

README.md

@@ -48,6 +48,7 @@ installation and postprocessing procedures. Its features include:
   * CFL-based and error-based time step control
 * Custom explicit time integration schemes
   * Maximized linear stability via paired explicit Runge-Kutta methods
+  * Relaxation Runge-Kutta methods for entropy-conservative time integration
 * Native support for differentiable programming
   * Forward mode automatic differentiation via [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
 * Periodic and weakly-enforced boundary conditions

docs/src/index.md

@@ -44,6 +44,7 @@ installation and postprocessing procedures. Its features include:
   * CFL-based and error-based time step control
 * Custom explicit time integration schemes
   * Maximized linear stability via paired explicit Runge-Kutta methods
+  * Relaxation Runge-Kutta methods for entropy-conservative time integration
 * Native support for differentiable programming
   * Forward mode automatic differentiation via [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
 * Periodic and weakly-enforced boundary conditions

docs/src/time_integration.md

@@ -215,4 +215,91 @@ Then, the stable CFL number can be computed as described above.
 
 - [`Trixi.PairedExplicitRK2`](@ref): Second-order PERK method with at least two stages.
 - [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages.
-- [`Trixi.PairedExplicitRK4`](@ref): Fourth-order PERK method with at least five stages.
+- [`Trixi.PairedExplicitRK4`](@ref): Fourth-order PERK method with at least five stages.
+
+## Relaxation Runge-Kutta Methods for Entropy-Conservative Time Integration
+
+While standard Runge-Kutta methods (or in fact the whole broad class of general linear methods such as multistep, additive, and partitioned Runge-Kutta methods) preserve linear solution invariants such as mass, momentum and energy, (assuming evolution in conserved variables $\boldsymbol u = (\rho, \rho v_i, \rho e)$) they do in general not preserve nonlinear solution invariants such as entropy.
+
+### The Notion of Entropy
+
+For an ideal gas with isentropic exponent $\gamma$, the thermodynamic entropy is given by
+```math
+s_\text{therm}(\boldsymbol u) = \ln \left( \frac{p}{\rho^\gamma} \right)
+```
+where $p$ is the pressure, $\rho$ the density, and $\gamma$ the ratio of specific heats.
+The mathematical entropy is then given by
+```math
+s(\boldsymbol u) \coloneqq - \underbrace{\rho}_{\equiv u_1} \cdot s_\text{therm}(\boldsymbol u) = - \rho \cdot \log \left( \frac{p(\boldsymbol u)}{\rho^\gamma} \right) \: .
+```
+The total entropy $\eta$ is then obtained by integrating the mathematical entropy $s$ over the domain $\Omega$:
+```math
+\eta(t) \coloneqq \eta \big(\boldsymbol u(t, \boldsymbol x) \big) = \int_{\Omega} s \big (\boldsymbol u(t, \boldsymbol x) \big ) \, \text{d} \boldsymbol x \tag{1}
+```
+
+For a semidiscretized partial differential equation (PDE) of the form
+```math
+\begin{align*}
+\boldsymbol U(t_0) &= \boldsymbol U_0, \\
+\boldsymbol U'(t) &= \boldsymbol F\big(t, \boldsymbol U(t) \big) \tag{2}
+\end{align*}
+```
+one can construct a discrete equivalent $H$ to (1) which is obtained by computing the mathematical entropy $s$ at every node of the mesh and then integrating it over the domain $\Omega$ by applying a quadrature rule:
+```math
+H(t) \coloneqq H\big(\boldsymbol U(t)\big) = \int_{\Omega} s \big(\boldsymbol U(t) \big) \, \text{d} \Omega
+```
+
+For a suitable spatial discretization (2) entropy-conservative systems such as the Euler equations preserve the total entropy $H(t)$ over time, i.e., 
+```math
+\frac{\text{d}}{\text{d} t} H \big(\boldsymbol U(t) \big ) 
+= 
+\left \langle \frac{\partial H(\boldsymbol U)}{\partial \boldsymbol U}, \frac{\text{d}}{\text{d} t} \boldsymbol U(t) \right \rangle 
+\overset{(2)}{=}
+\left \langle \frac{\partial H(\boldsymbol U)}{\partial \boldsymbol U}, \boldsymbol F\big(t, \boldsymbol U(t) \big) \right \rangle = 0 \tag{3}
+```
+while entropy-stable discretiations of entropy-diffusive systems such as the Navier-Stokes equations ensure that the total entropy decays over time, i.e., 
+```math
+\left \langle \frac{\partial H(\boldsymbol U)}{\partial \boldsymbol U}, \boldsymbol F(t, \boldsymbol U) \right \rangle \leq 0 \tag{4}
+```
+
+### Ensuring Entropy-Conservation/Stability with Relaxation Runge-Kutta Methods
+
+Evolving the ordinary differential equation (ODE) for the entropy (2) with a Runge-Kutta scheme gives
+```math
+H_{n+1} = H_n + \Delta t \sum_{i=1}^S b_i \, \left\langle \frac{\partial 
+H(\boldsymbol U_{n, i})
+}{\partial \boldsymbol U}, \boldsymbol F(\boldsymbol U_{n, i}) \right\rangle \tag{5}
+```
+which preserves (3) and (4).
+In practice, however, we evolve the conserved variables $\boldsymbol U$ which results in 
+```math
+\boldsymbol U_{n+1} = \boldsymbol U_n + \Delta t \sum_{i=1}^S b_i \boldsymbol F(\boldsymbol U_{n, i})
+```
+and in particular for the entropy $H$
+```math
+H(\boldsymbol U_{n+1}) = H\left( \boldsymbol U_n + \Delta t \sum_{i=1}^S b_i \boldsymbol F(\boldsymbol U_{n, i}) \right) \neq H_{n+1} \text{ computed from (5)}
+```
+
+To resolve the difference $H(\boldsymbol U_{n+1}) - H_{n+1}$ Ketcheson, Ranocha and collaborators have introduced *relaxation* Runge-Kutta methods in a series of publications, see for instance
+- [Ketcheson (2019)](https://doi.org/10.1137/19M1263662): Relaxation Runge-Kutta Methods: Conservation and Stability for Inner-Product Norms
+- [Ranocha et al. (2020)](https://doi.org/10.1137/19M1263480): Relaxation Runge-Kutta methods: Fully discrete explicit entropy-stable schemes for the compressible Euler and Navier-Stokes equations
+- [Ranocha, Lóczi, and Ketcheson (2020)](https://doi.org/10.1007/s00211-020-01158-4): General relaxation methods for initial-value problems with application to multistep schemes
+
+Almost miraculously, it suffices to introduce a single parameter $\gamma$ in the final update step of the Runge-Kutta method to ensure that the properties of the spatial discretization are preserved, i.e., 
+```math
+H \big(\boldsymbol U_{n+1}( \gamma  ) \big) 
+\overset{!}{=} 
+H(\boldsymbol U_n) + 
+\gamma  \Delta t \sum_{i=1}^S b_i 
+\left \langle 
+\frac{\partial H(\boldsymbol U_{n, i})}{\partial \boldsymbol U_{n, i}}, \boldsymbol  F(\boldsymbol U_{n, i}) 
+\right  \rangle
+\tag{6}
+```
+This comes only at the price that one needs to solve the scalar nonlinear equation (6) for $\gamma$ at every time step.
+To do so, [`Trixi.RelaxationSolverNewton`](@ref) is implemented in Trixi.jl.
+These can then be supplied to the relaxation time algorithms such as [`Trixi.RelaxationRalston3`](@ref) and [`Trixi.RelaxationRK44`](@ref) via specifying the `relaxation_solver` keyword argument:
+```julia
+ode_algorithm = Trixi.RelaxationRK44(solver = Trixi.RelaxationSolverNewton())
+ode_algorithm = Trixi.RelaxationRalston3(solver = Trixi.RelaxationSolverNewton())
+```

examples/structured_1d_dgsem/elixir_euler_weak_blast_er.jl

@@ -0,0 +1,53 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+# Volume flux stabilizes the simulation - in contrast to standard DGSEM with 
+# `surface_flux = flux_ranocha` only which crashes.
+solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
+               volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
+
+coordinates_min = -2.0
+coordinates_max = 2.0
+cells_per_dimension = 32
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+initial_condition = initial_condition_weak_blast_wave
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_errors = Symbol[], # Switch off error computation
+                                     # Note: `entropy` defaults to mathematical entropy
+                                     analysis_integrals = (entropy,),
+                                     analysis_filename = "entropy_ER.dat",
+                                     save_analysis = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.25)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# Ensure exact entropy conservation by employing a relaxation Runge-Kutta method
+relaxation_solver = Trixi.RelaxationSolverNewton(max_iterations = 5,
+                                                 root_tol = eps(Float64),
+                                                 gamma_tol = eps(Float64))
+ode_alg = Trixi.RelaxationRK44(relaxation_solver = relaxation_solver)
+
+sol = Trixi.solve(ode, ode_alg,
+                  dt = 42.0, save_everystep = false, callback = callbacks);

examples/tree_3d_dgsem/elixir_advection_er.jl

@@ -0,0 +1,51 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (1.0, 1.0, 1.0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+solver = DGSEM(polydeg = 2, surface_flux = flux_central,
+               volume_integral = VolumeIntegralFluxDifferencing(flux_central)) # Entropy-conservative setup
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = (1.0, 1.0, 1.0)
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+ode = semidiscretize(semi, (0.0, 1.0))
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_errors = Symbol[], # Switch off error computation
+                                     # Note: `entropy` defaults to mathematical entropy
+                                     analysis_integrals = (entropy,),
+                                     analysis_filename = "entropy_ER.dat",
+                                     save_analysis = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+relaxation_solver = Trixi.RelaxationSolverNewton(max_iterations = 3,
+                                                 root_tol = eps(Float64),
+                                                 gamma_tol = eps(Float64))
+ode_alg = Trixi.RelaxationRalston3(relaxation_solver = relaxation_solver)
+
+sol = Trixi.solve(ode, ode_alg,
+                  dt = 42.0, save_everystep = false, callback = callbacks);