New Paired Explicit Runge-Kutta Integrator: Third Order

NEWS.md

@@ -4,6 +4,14 @@ 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 in the v0.9 lifecycle
+
+#### Added
+
+- New time integrator `PairedExplicitRK3`, implementing the third-order paired explicit Runge-Kutta
+  method with [Convex.jl](https://github.com/jump-dev/Convex.jl), [ECOS.jl](https://github.com/jump-dev/ECOS.jl),
+  and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008])
+
 ## Changes when updating to v0.9 from v0.8.x
 
 #### Added

Project.toml

@@ -36,6 +36,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
@@ -55,10 +56,12 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 
 [extensions]
 TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
+TrixiNLsolveExt = "NLsolve"
 
 [compat]
 Accessors = "0.1.12"
@@ -82,6 +85,7 @@ LoopVectorization = "0.12.151"
 MPI = "0.20"
 Makie = "0.19, 0.20"
 MuladdMacro = "0.2.2"
+NLsolve = "4.5.1"
 Octavian = "0.3.21"
 OffsetArrays = "1.12"
 P4est = "0.4.9"
@@ -96,6 +100,7 @@ Requires = "1.1"
 SciMLBase = "1.90, 2"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
+StableRNGs = "1.0.2"
 StartUpDG = "0.17.7, 1.1.5"
 Static = "0.8.7"
 StaticArrayInterface = "1.4"
@@ -116,3 +121,4 @@ julia = "1.8"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"

examples/structured_1d_dgsem/elixir_burgers_perk3.jl

@@ -0,0 +1,70 @@
+# Convex and ECOS are imported because they are used for finding the optimal time step and optimal 
+# monomial coefficients in the stability polynomial of P-ERK time integrators.
+using Convex, ECOS
+
+# NLsolve is imported to solve the system of nonlinear equations to find the coefficients
+# in the Butcher tableau in the third order P-ERK time integrator.
+using NLsolve
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the (inviscid) Burgers equation
+
+equations = InviscidBurgersEquation1D()
+
+initial_condition = initial_condition_convergence_test
+
+# Create DG solver with polynomial degree = 4 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (0.0,) # minimum coordinate
+coordinates_max = (1.0,) # maximum coordinate
+cells_per_dimension = (64,)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 200
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.1,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+# Construct third order paired explicit Runge-Kutta method with 8 stages for given simulation setup.
+# Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used 
+# in calculating the polynomial coefficients in the ODE algorithm.                                     
+ode_algorithm = Trixi.PairedExplicitRK3(8, tspan, semi)
+
+cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
+# 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)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/tree_1d_dgsem/elixir_advection_perk2.jl

@@ -1,5 +1,8 @@
 
+# Convex and ECOS are imported because they are used for finding the optimal time step and optimal 
+# monomial coefficients in the stability polynomial of P-ERK time integrators.
 using Convex, ECOS
+
 using OrdinaryDiffEq
 using Trixi
 

ext/TrixiNLsolveExt.jl

@@ -0,0 +1,128 @@
+# Package extension for adding NLsolve-based features to Trixi.jl
+module TrixiNLsolveExt
+
+# Required for finding coefficients in Butcher tableau in the third order of 
+# P-ERK scheme integrators
+if isdefined(Base, :get_extension)
+    using NLsolve: nlsolve
+else
+    # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
+    using ..NLsolve: nlsolve
+end
+
+# We use a random initialization of the nonlinear solver.
+# To make the tests reproducible, we need to seed the RNG
+using StableRNGs: StableRNG, rand
+
+# Use functions that are to be extended and additional symbols that are not exported
+using Trixi: Trixi, compute_c_coeffs, @muladd
+
+# 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
+
+# Compute residuals for nonlinear equations to match a stability polynomial with given coefficients,
+# in order to find A-matrix in the Butcher-Tableau
+function PairedExplicitRK3_butcher_tableau_objective_function!(c_eq, a_unknown,
+                                                               num_stages,
+                                                               num_stage_evals,
+                                                               monomial_coeffs,
+                                                               cS2)
+    c_ts = compute_c_coeffs(num_stages, cS2) # ts = timestep
+    # For explicit methods, a_{1,1} = 0 and a_{2,1} = c_2 (Butcher's condition)
+    a_coeff = [0, c_ts[2], a_unknown...]
+    # Equality constraint array that ensures that the stability polynomial computed from 
+    # the to-be-constructed Butcher-Tableau matches the monomial coefficients of the 
+    # optimized stability polynomial.
+    # For details, see Chapter 4.3, Proposition 3.2, Equation (3.3) from 
+    # Hairer, Wanner: Solving Ordinary Differential Equations 2
+    # DOI: 10.1007/978-3-662-09947-6
+
+    # Lower-order terms: Two summands present
+    for i in 1:(num_stage_evals - 4)
+        term1 = a_coeff[num_stage_evals - 1]
+        term2 = a_coeff[num_stage_evals]
+        for j in 1:i
+            term1 *= a_coeff[num_stage_evals - 1 - j]
+            term2 *= a_coeff[num_stage_evals - j]
+        end
+        term1 *= c_ts[num_stages - 2 - i] * 1 / 6 # 1 / 6 = b_{S-1}
+        term2 *= c_ts[num_stages - 1 - i] * 2 / 3 # 2 / 3 = b_S
+
+        c_eq[i] = monomial_coeffs[i] - (term1 + term2)
+    end
+
+    # Highest coefficient: Only one term present
+    i = num_stage_evals - 3
+    term2 = a_coeff[num_stage_evals]
+    for j in 1:i
+        term2 *= a_coeff[num_stage_evals - j]
+    end
+    term2 *= c_ts[num_stages - 1 - i] * 2 / 3 # 2 / 3 = b_S
+
+    c_eq[i] = monomial_coeffs[i] - term2
+    # Third-order consistency condition (Cf. eq. (27) from https://doi.org/10.1016/j.jcp.2022.111470
+    c_eq[num_stage_evals - 2] = 1 - 4 * a_coeff[num_stage_evals] -
+                                a_coeff[num_stage_evals - 1]
+end
+
+# Find the values of the a_{i, i-1} in the Butcher tableau matrix A by solving a system of
+# non-linear equations that arise from the relation of the stability polynomial to the Butcher tableau.
+# For details, see Proposition 3.2, Equation (3.3) from 
+# Hairer, Wanner: Solving Ordinary Differential Equations 2
+function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, monomial_coeffs,
+                                               c_s2, c;
+                                               verbose, max_iter = 100000)
+
+    # Define the objective_function
+    function objective_function!(c_eq, x)
+        return PairedExplicitRK3_butcher_tableau_objective_function!(c_eq, x,
+                                                                     num_stages,
+                                                                     num_stages,
+                                                                     monomial_coeffs,
+                                                                     c_s2)
+    end
+
+    # RealT is determined as the type of the first element in monomial_coeffs to ensure type consistency
+    RealT = typeof(monomial_coeffs[1])
+
+    # To ensure consistency and reproducibility of results across runs, we use 
+    # a seeded random initial guess.
+    rng = StableRNG(555)
+
+    for _ in 1:max_iter
+        # Due to the nature of the nonlinear solver, different initial guesses can lead to 
+        # small numerical differences in the solution.
+
+        x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stages - 2)
+
+        sol = nlsolve(objective_function!, x0, method = :trust_region,
+                      ftol = 4.0e-16, # Enforce objective up to machine precision
+                      iterations = 10^4, xtol = 1.0e-13, autodiff = :forward)
+
+        a_unknown = sol.zero # Retrieve solution (root = zero)
+
+        # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) 
+        # and also c[i] - a[i, i-1] >= 0.0 since all coefficients should be non-negative
+        # to avoid downwinding of numerical fluxes.
+        is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown) &&
+                       all(!isnan(c[i] - a_unknown[i - 2]) &&
+                           c[i] - a_unknown[i - 2] >= 0 for i in eachindex(c) if i > 2)
+
+        if is_sol_valid
+            return a_unknown
+        else
+            if verbose
+                println("Solution invalid. Restart the process of solving non-linear system of equations again.")
+            end
+        end
+    end
+
+    error("Maximum number of iterations ($max_iter) reached. Cannot find valid sets of coefficients.")
+end
+end # @muladd
+
+end # module TrixiNLsolveExt

src/Trixi.jl

@@ -19,7 +19,8 @@ module Trixi
 # (standard library packages first, other packages next, all of them sorted alphabetically)
 
 using Accessors: @reset
-using LinearAlgebra: LinearAlgebra, Diagonal, diag, dot, mul!, norm, cross, normalize, I,
+using LinearAlgebra: LinearAlgebra, Diagonal, diag, dot, eigvals, mul!, norm, cross,
+                     normalize, I,
                      UniformScaling, det
 using Printf: @printf, @sprintf, println
 using SparseArrays: AbstractSparseMatrix, AbstractSparseMatrixCSC, sparse, droptol!,
@@ -40,6 +41,7 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   get_proposed_dt, set_proposed_dt!,
                   terminate!, remake, add_tstop!, has_tstop, first_tstop
 
+using DelimitedFiles: readdlm
 using Downloads: Downloads
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
@@ -320,6 +322,12 @@ function __init__()
         end
     end
 
+    @static if !isdefined(Base, :get_extension)
+        @require NLsolve="2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
+            include("../ext/TrixiNLsolveExt.jl")
+        end
+    end
+
     # FIXME upstream. This is a hacky workaround for
     #       https://github.com/trixi-framework/Trixi.jl/issues/628
     #       https://github.com/trixi-framework/Trixi.jl/issues/1185