New Paired Explicit Runge-Kutta Integrator: Third Order

NEWS.md

@@ -29,6 +29,9 @@ for human readability.
 - New time integrator `PairedExplicitRK2`, implementing the second-order paired explicit Runge-Kutta
   method with [Convex.jl](https://github.com/jump-dev/Convex.jl) and [ECOS.jl](https://github.com/jump-dev/ECOS.jl) ([#1908])
 - Add subcell limiting support for `StructuredMesh` ([#1946]).
+- 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.7 from v0.6.x
 

Project.toml

@@ -34,6 +34,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"
@@ -50,13 +51,15 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [weakdeps]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 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]
-TrixiMakieExt = "Makie"
 TrixiConvexECOSExt = ["Convex", "ECOS"]
+TrixiMakieExt = "Makie"
+TrixiNLsolveExt = "NLsolve"
 
 [compat]
 CodeTracking = "1.0.5"
@@ -79,6 +82,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"
@@ -92,6 +96,7 @@ Requires = "1.1"
 SciMLBase = "1.90, 2"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
+StableRNGs = "1.0.2"
 StartUpDG = "0.17.7"
 Static = "0.8.7"
 StaticArrayInterface = "1.4"
@@ -109,6 +114,7 @@ UUIDs = "1.6"
 julia = "1.8"
 
 [extras]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 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,66 @@
+# 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 a 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)
+
+stepsize_callback = StepsizeCallback(cfl = 3.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# Optimize 8-stage, third order P-ERK scheme for this semidiscretization
+ode_algorithm = Trixi.PairedExplicitRK3(8, tspan, semi)
+
+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(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 Chapter4.3, Proposition 3.2, Equation (3.3) from 
+    # Hairer, Wanner: Solving Ordinary Differential Equations 2
+    c_eq = zeros(num_stage_evals - 2) # Add equality constraint that cS2 is equal to 1
+    # 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]
+
+    return c_eq
+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(x)
+        return PairedExplicitRK3_butcher_tableau_objective_function(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)
+
+    # Flag for criteria going beyond satisfaction of non-linear equations
+    is_sol_valid = false
+
+    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)
+
+        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
+        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

@@ -319,6 +319,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

src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl

@@ -233,11 +233,11 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F,
 end
 
 """
-    add_tstop!(integrator::PairedExplicitRK2Integrator, t)
+    add_tstop!(integrator::AbstractPairedExplicitRKSingleIntegrator, t)
 Add a time stop during the time integration process.
 This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution.
 """
-function add_tstop!(integrator::PairedExplicitRK2Integrator, t)
+function add_tstop!(integrator::AbstractPairedExplicitRKSingleIntegrator, t)
     integrator.tdir * (t - integrator.t) < zero(integrator.t) &&
         error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
     # We need to remove the first entry of tstops when a new entry is added.
@@ -248,8 +248,8 @@ function add_tstop!(integrator::PairedExplicitRK2Integrator, t)
     push!(integrator.opts.tstops, integrator.tdir * t)
 end
 
-has_tstop(integrator::PairedExplicitRK2Integrator) = !isempty(integrator.opts.tstops)
-first_tstop(integrator::PairedExplicitRK2Integrator) = first(integrator.opts.tstops)
+has_tstop(integrator::AbstractPairedExplicitRKSingleIntegrator) = !isempty(integrator.opts.tstops)
+first_tstop(integrator::AbstractPairedExplicitRKSingleIntegrator) = first(integrator.opts.tstops)
 
 # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771)
 function Base.getproperty(integrator::PairedExplicitRK, field::Symbol)