diff --git a/NEWS.md b/NEWS.md index cba0d3f86d..d561562b97 100644 --- a/NEWS.md +++ b/NEWS.md @@ -63,6 +63,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 diff --git a/Project.toml b/Project.toml index 285d29ce63..63febf9c65 100644 --- a/Project.toml +++ b/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" diff --git a/examples/structured_1d_dgsem/elixir_burgers_perk3.jl b/examples/structured_1d_dgsem/elixir_burgers_perk3.jl new file mode 100644 index 0000000000..4d031db7ae --- /dev/null +++ b/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() diff --git a/examples/structured_1d_dgsem/elixir_burgers_perk3_optimal_cfl.jl b/examples/structured_1d_dgsem/elixir_burgers_perk3_optimal_cfl.jl new file mode 100644 index 0000000000..b59543fbfd --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_burgers_perk3_optimal_cfl.jl @@ -0,0 +1,69 @@ +# 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) + +# Construct second 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 +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() diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl index 7db461a079..68befbcfe8 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl +++ b/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 diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl new file mode 100644 index 0000000000..4de8feb25a --- /dev/null +++ b/ext/TrixiNLsolveExt.jl @@ -0,0 +1,129 @@ +# 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 Chapter4.3, Proposition 3.2, Equation (3.3) from + # Hairer, Wanner: Solving Ordinary Differential Equations 2 + + # 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) + + # 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, 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 + 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 diff --git a/src/Trixi.jl b/src/Trixi.jl index d9b9245918..d51c25d9c9 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -320,6 +320,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 diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index ad385e6df2..9934208bf2 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -256,11 +256,11 @@ function calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode) 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. @@ -271,8 +271,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) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl new file mode 100644 index 0000000000..6f0a9a9629 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -0,0 +1,368 @@ +# 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. +using DelimitedFiles: readdlm + +@muladd begin +#! format: noindent + +# Initialize Butcher array abscissae c for PairedExplicitRK3 based on SSPRK33 base method +function compute_c_coeffs(num_stages, cS2) + c = zeros(num_stages) + + # Last timesteps as for SSPRK33, see motivation in Section 3.3 of + # https://doi.org/10.1016/j.jcp.2024.113223 + c[num_stages - 1] = 1.0f0 + c[num_stages] = 0.5f0 + + # Linear increasing timestep for remainder + for i in 2:(num_stages - 2) + c[i] = cS2 * (i - 1) / (num_stages - 3) + end + + return c +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 +# using a list of eigenvalues +function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, + eig_vals::Vector{ComplexF64}; + verbose = false, cS2) + # Initialize array of c + c = compute_c_coeffs(num_stages, cS2) + + # Initialize the array of our solution + a_unknown = zeros(num_stages - 2) + + # Special case of e = 3 + if num_stages == 3 + a_unknown = [0.25] + else + # Calculate coefficients of the stability polynomial in monomial form + consistency_order = 3 + dtmax = tspan[2] - tspan[1] + dteps = 1.0f-9 + + num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + + monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, + num_eig_vals, num_stages, + dtmax, dteps, + eig_vals; verbose) + monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, + num_stages) + + # Solve the nonlinear system of equations from monomial coefficient and + # Butcher array abscissae c to find Butcher matrix A + # This function is extended in TrixiNLsolveExt.jl + a_unknown = solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, + monomial_coeffs, cS2, c; + verbose) + end + # Fill A-matrix in P-ERK style + a_matrix = zeros(num_stages - 2, 2) + a_matrix[:, 1] = c[3:end] + a_matrix[:, 1] -= a_unknown + a_matrix[:, 2] = a_unknown + + return a_matrix, c, dt_opt +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 +# using provided values of coefficients a in A-matrix of Butcher tableau +function compute_PairedExplicitRK3_butcher_tableau(num_stages, + base_path_a_coeffs::AbstractString; + cS2) + + # Initialize array of c + c = compute_c_coeffs(num_stages, cS2) + + # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) + a_coeffs_max = num_stages - 2 + + a_matrix = zeros(a_coeffs_max, 2) + a_matrix[:, 1] = c[3:end] + + path_a_coeffs = joinpath(base_path_a_coeffs, + "a_" * string(num_stages) * ".txt") + + @assert isfile(path_a_coeffs) "Couldn't find file $path_a_coeffs" + a_coeffs = readdlm(path_a_coeffs, Float64) + num_a_coeffs = size(a_coeffs, 1) + + @assert num_a_coeffs == a_coeffs_max + # Fill A-matrix in P-ERK style + a_matrix[:, 1] -= a_coeffs + a_matrix[:, 2] = a_coeffs + + return a_matrix, c +end + +@doc raw""" + PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt; + cS2 = 1.0f0) + PairedExplicitRK3(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, cS2 = 1.0f0) + PairedExplicitRK3(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, cS2 = 1.0f0) + + Parameters: + - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method. + - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in + the Butcher tableau of the Runge Kutta method. + The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks. + - `dt_opt` (`Float64`): Optimal time step size for the simulation setup. + - `tspan`: Time span of the simulation. + - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. + - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the + equation has been semidiscretized. + - `verbose` (`Bool`, optional): Verbosity flag, default is false. + - `cS2` (`Float64`, optional): Value of c in the Butcher tableau at c_{s-2}, when + s is the number of stages, default is 1.0f0. + +The following structures and methods provide an implementation of +the third-order paired explicit Runge-Kutta (P-ERK) method +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). +The original paper is +- Nasab, Vermeire (2022) +Third-order Paired Explicit Runge-Kutta schemes for stiff systems of equations +[DOI: 10.1016/j.jcp.2022.111470](https://doi.org/10.1016/j.jcp.2022.111470) +While the changes to SSPRK33 base-scheme are described in +- Doehring, Schlottke-Lakemper, Gassner, Torrilhon (2024) +Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively Refined Meshes for Compressible Fluid Dynamics +[DOI: 10.1016/j.jcp.2024.113223](https://doi.org/10.1016/j.jcp.2024.113223) +""" +mutable struct PairedExplicitRK3 <: AbstractPairedExplicitRKSingle + const num_stages::Int # S + + a_matrix::Matrix{Float64} + c::Vector{Float64} + dt_opt::Float64 +end # struct PairedExplicitRK3 + +# Constructor for previously computed A Coeffs +function PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt; + cS2 = 1.0f0) + a_matrix, c = compute_PairedExplicitRK3_butcher_tableau(num_stages, + base_path_a_coeffs; + cS2) + + return PairedExplicitRK3(num_stages, a_matrix, c, dt_opt) +end + +# Constructor that computes Butcher matrix A coefficients from a semidiscretization +function PairedExplicitRK3(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, cS2 = 1.0f0) + eig_vals = eigvals(jacobian_ad_forward(semi)) + + return PairedExplicitRK3(num_stages, tspan, eig_vals; verbose, cS2) +end + +# Constructor that calculates the coefficients with polynomial optimizer from a list of eigenvalues +function PairedExplicitRK3(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, cS2 = 1.0f0) + a_matrix, c, dt_opt = compute_PairedExplicitRK3_butcher_tableau(num_stages, + tspan, + eig_vals; + verbose, cS2) + return PairedExplicitRK3(num_stages, a_matrix, c, dt_opt) +end + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77 +# This implements the interface components described at +# https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1 +# which are used in Trixi.jl. +mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F, Alg, + PairedExplicitRKOptions} <: + AbstractPairedExplicitRKSingleIntegrator + u::uType + du::uType + u_tmp::uType + t::RealT + tdir::RealT + dt::RealT # current time step + dtcache::RealT # manually set time step + iter::Int # current number of time steps (iteration) + p::Params # will be the semidiscretization from Trixi + sol::Sol # faked + f::F + alg::Alg # This is our own class written above; Abbreviation for ALGorithm + opts::PairedExplicitRKOptions + finalstep::Bool # added for convenience + dtchangeable::Bool + force_stepfail::Bool + # PairedExplicitRK stages: + k1::uType + k_higher::uType +end + +function init(ode::ODEProblem, alg::PairedExplicitRK3; + dt, callback = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(u0) + u_tmp = zero(u0) + + # PairedExplicitRK stages + k1 = zero(u0) + k_higher = zero(u0) + + t0 = first(ode.tspan) + tdir = sign(ode.tspan[end] - ode.tspan[1]) + iter = 0 + + integrator = PairedExplicitRK3Integrator(u0, du, u_tmp, t0, tdir, dt, dt, iter, + ode.p, + (prob = ode,), ode.f, alg, + PairedExplicitRKOptions(callback, + ode.tspan; + kwargs...), + false, true, false, + k1, k_higher) + + # initialize callbacks + if callback isa CallbackSet + for cb in callback.continuous_callbacks + error("Continuous callbacks are unsupported with paired explicit Runge-Kutta methods.") + end + for cb in callback.discrete_callbacks + cb.initialize(cb, integrator.u, integrator.t, integrator) + end + elseif !isnothing(callback) + error("unsupported") + end + + return integrator +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::PairedExplicitRK3; + dt, callback = nothing, kwargs...) + integrator = init(ode, alg, dt = dt, callback = callback; kwargs...) + + # Start actual solve + solve!(integrator) +end + +function solve!(integrator::PairedExplicitRK3Integrator) + @unpack prob = integrator.sol + + integrator.finalstep = false + + @trixi_timeit timer() "main loop" while !integrator.finalstep + step!(integrator) + end # "main loop" timer + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) +end + +function step!(integrator::PairedExplicitRK3Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end + + modify_dt_for_tstops!(integrator) + + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end + + @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin + # k1 + integrator.f(integrator.du, integrator.u, prob.p, integrator.t) + @threaded for i in eachindex(integrator.du) + integrator.k1[i] = integrator.du[i] * integrator.dt + end + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + # k2 + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[2] * integrator.dt) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + # Higher stages + for stage in 3:(alg.num_stages - 1) + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.a_matrix[stage - 2, 1] * + integrator.k1[i] + + alg.a_matrix[stage - 2, 2] * + integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[stage] * integrator.dt) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + end + + # Last stage + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.a_matrix[alg.num_stages - 2, 1] * + integrator.k1[i] + + alg.a_matrix[alg.num_stages - 2, 2] * + integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[alg.num_stages] * integrator.dt) + + @threaded for i in eachindex(integrator.u) + # "Own" PairedExplicitRK based on SSPRK33. + # Note that 'k_higher' carries the values of K_{S-1} + # and that we construct 'K_S' "in-place" from 'integrator.du' + integrator.u[i] += (integrator.k1[i] + integrator.k_higher[i] + + 4.0 * integrator.du[i] * integrator.dt) / 6.0 + end + end # PairedExplicitRK step timer + + integrator.iter += 1 + integrator.t += integrator.dt + + # handle callbacks + if callbacks isa CallbackSet + for cb in callbacks.discrete_callbacks + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) + end + end + end + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end +end + +# used for AMR (Adaptive Mesh Refinement) +function Base.resize!(integrator::PairedExplicitRK3Integrator, new_size) + resize!(integrator.u, new_size) + resize!(integrator.du, new_size) + resize!(integrator.u_tmp, new_size) + + resize!(integrator.k1, new_size) + resize!(integrator.k_higher, new_size) +end +end # @muladd diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl index b73ea75831..8f8cacadd7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -7,6 +7,12 @@ # Basic implementation of the second-order paired explicit Runge-Kutta (PERK) method include("methods_PERK2.jl") +include("methods_PERK3.jl") # Define all of the functions necessary for polynomial optimizations include("polynomial_optimizer.jl") + +# Add definitions of functions related to polynomial optimization by NLsolve here +# such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package +# extension or by the NLsolve-specific code loaded by Requires.jl +function solve_a_butcher_coeffs_unknown! end end # @muladd diff --git a/test/Project.toml b/test/Project.toml index c8ae33a40a..f8bcff947c 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,10 +10,12 @@ FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] @@ -28,10 +30,12 @@ FFMPEG = "0.4" ForwardDiff = "0.10.24" LinearAlgebra = "1" MPI = "0.20" +NLsolve = "4.5.1" OrdinaryDiffEq = "6.49.1" Plots = "1.19" Printf = "1" Random = "1" +StableRNGs = "1.0.2" Test = "1" [preferences.OrdinaryDiffEq] diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index f64b8c9c06..6bc4b9ccab 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -58,6 +58,57 @@ end end end +@trixi_testset "elixir_burgers_perk3.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3.jl"), + l2=[4.12066275835687e-6], + linf=[2.538190787615413e-5], + atol=1.0e-6) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end + +# Testing the third-order paired explicit Runge-Kutta (PERK) method with its optimal CFL number +@trixi_testset "elixir_burgers_perk3_optimal_cfl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3_optimal_cfl.jl"), + l2=[3.8156922097242205e-6], + linf=[2.1962957979626552e-5], + atol=1.0e-6) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end + +# Testing the third-order paired explicit Runge-Kutta (PERK) method without stepsize callback +@trixi_testset "elixir_burgers_perk3.jl(fixed time step)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3.jl"), + dt=2.0e-3, + tspan=(0.0, 2.0), + save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8), + callbacks=CallbackSet(summary_callback, save_solution, + analysis_callback, alive_callback), + l2=[5.726144786001842e-7], + linf=[3.430730019182704e-6]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end + @trixi_testset "elixir_euler_sedov.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[3.67478226e-01, 3.49491179e-01, 8.08910759e-01], diff --git a/test/test_unit.jl b/test/test_unit.jl index bccdcf8faa..e4d2ff24d9 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -10,6 +10,10 @@ using DelimitedFiles: readdlm using Convex: Convex using ECOS: Optimizer +# Use NLsolve to load the extension that extends functions for testing +# PERK Single p3 Constructors +using NLsolve: nlsolve + include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists @@ -1699,6 +1703,42 @@ end 0.13942836392866081 0.3605716360713392], atol = 1e-13) end +@testset "PERK Single p3 Constructors" begin + path_coeff_file = mktempdir() + Trixi.download("https://gist.githubusercontent.com/warisa-r/0796db36abcd5abe735ac7eebf41b973/raw/32889062fd5dcf7f450748f4f5f0797c8155a18d/a_8_8.txt", + joinpath(path_coeff_file, "a_8.txt")) + + ode_algorithm = Trixi.PairedExplicitRK3(8, path_coeff_file, 42) # dummy optimal time step (dt_opt plays no role in determining `a_matrix`) + + @test isapprox(ode_algorithm.a_matrix, + [0.33551678438002486 0.06448322158043965 + 0.49653494442225443 0.10346507941960345 + 0.6496890912144586 0.15031092070647037 + 0.789172498521197 0.21082750147880308 + 0.7522972036571336 0.2477027963428664 + 0.31192569908571666 0.18807430091428337], atol = 1e-13) + + Trixi.download("https://gist.githubusercontent.com/warisa-r/8d93f6a3ae0635e13b9f51ee32ab7fff/raw/54dc5b14be9288e186b745facb5bbcb04d1476f8/EigenvalueList_Refined2.txt", + joinpath(path_coeff_file, "spectrum.txt")) + + eig_vals = readdlm(joinpath(path_coeff_file, "spectrum.txt"), ComplexF64) + tspan = (0.0, 1.0) + ode_algorithm = Trixi.PairedExplicitRK3(13, tspan, vec(eig_vals)) + + @test isapprox(ode_algorithm.a_matrix, + [0.19258815508348048 0.007411847896751755 + 0.28723293302534425 0.012767078895584727 + 0.38017717306156973 0.01982283289889476 + 0.4706748926615584 0.02932510733844162 + 0.5575748313419124 0.04242519249994544 + 0.6390917539959684 0.06090823408310269 + 0.7124876783811593 0.08751233353976957 + 0.7736369902636945 0.12636298589444758 + 0.8161315438506253 0.18386845614937475 + 0.7532704353232954 0.24672956467670462 + 0.3116823911691762 0.18831760883082385], atol = 1e-13) +end + @testset "Sutherlands Law" begin function mu(u, equations) T_ref = 291.15