diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index 234939a8017..46ce5d71e36 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -72,7 +72,7 @@ jobs: - uses: julia-actions/cache@v1 - uses: julia-actions/julia-downgrade-compat@v1 with: - skip: LinearAlgebra,Printf,SparseArrays,UUIDs,DiffEqBase + skip: LinearAlgebra,Printf,SparseArrays,UUIDs,DiffEqBase,DelimitedFiles projects: ., test - uses: julia-actions/julia-buildpkg@v1 env: diff --git a/NEWS.md b/NEWS.md index e2902229f71..ebc8d9cda39 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,12 +7,13 @@ for human readability. ## Changes in the v0.7 lifecycle #### Added -- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension - to 1D and 3D on `TreeMesh` ([#1855], [#1873]). +- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh` ([#1855], [#1873]). - Implementation of 1D Linearized Euler Equations ([#1867]). - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients ([#1812]). - Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed ([#1835]). - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` ([#1792]). +- 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]) ## Changes when updating to v0.7 from v0.6.x diff --git a/Project.toml b/Project.toml index 9ebb9c307a1..b605cb2c6b1 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.7.15-pre" CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" @@ -50,17 +51,23 @@ 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" [extensions] TrixiMakieExt = "Makie" +TrixiConvexECOSExt = ["Convex", "ECOS"] [compat] CodeTracking = "1.0.5" ConstructionBase = "1.3" +Convex = "0.15.4" DataStructures = "0.18.15" +DelimitedFiles = "1" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" +ECOS = "1.1.2" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.24" @@ -103,3 +110,5 @@ julia = "1.8" [extras] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +Convex = "f65535da-76fb-5f13-bab9-19810c17039a" +ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl new file mode 100644 index 00000000000..87a61515f82 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl @@ -0,0 +1,66 @@ + +using Convex, ECOS +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = 1.0 +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = -1.0 # minimum coordinate +coordinates_max = 1.0 # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 20.0 +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 2.5) + +alive_callback = AliveCallback(alive_interval = analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, + alive_callback, + analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# Construct second order paired explicit Runge-Kutta method with 6 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.PairedExplicitRK2(6, 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/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl new file mode 100644 index 00000000000..fac127699ce --- /dev/null +++ b/ext/TrixiConvexECOSExt.jl @@ -0,0 +1,158 @@ +# Package extension for adding Convex-based features to Trixi.jl +module TrixiConvexECOSExt + +# Required for coefficient optimization in P-ERK scheme integrators +if isdefined(Base, :get_extension) + using Convex: MOI, solve!, Variable, minimize, evaluate + using ECOS: Optimizer +else + # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl + using ..Convex: MOI, solve!, Variable, minimize, evaluate + using ..ECOS: Optimizer +end + +# Use other necessary libraries +using LinearAlgebra: eigvals + +# Use functions that are to be extended and additional symbols that are not exported +using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @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 + +# Undo normalization of stability polynomial coefficients by index factorial +# relative to consistency order. +function Trixi.undo_normalization!(gamma_opt, consistency_order, num_stage_evals) + for k in (consistency_order + 1):num_stage_evals + gamma_opt[k - consistency_order] = gamma_opt[k - consistency_order] / + factorial(k) + end + return gamma_opt +end + +# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency +# order, including contributions from free coefficients for higher orders, and +# return the maximum absolute value +function stability_polynomials!(pnoms, consistency_order, num_stage_evals, + normalized_powered_eigvals_scaled, + gamma) + num_eig_vals = length(pnoms) + + # Initialize with zero'th order (z^0) coefficient + for i in 1:num_eig_vals + pnoms[i] = 1.0 + end + + # First `consistency_order` terms of the exponential + for k in 1:consistency_order + for i in 1:num_eig_vals + pnoms[i] += normalized_powered_eigvals_scaled[i, k] + end + end + + # Contribution from free coefficients + for k in (consistency_order + 1):num_stage_evals + pnoms += gamma[k - consistency_order] * normalized_powered_eigvals_scaled[:, k] + end + + # For optimization only the maximum is relevant + return maximum(abs(pnoms)) +end + +#= +The following structures and methods provide a simplified implementation to +discover optimal stability polynomial for a given set of `eig_vals` +These are designed for the one-step (i.e., Runge-Kutta methods) integration of initial value ordinary +and partial differential equations. + +- Ketcheson and Ahmadia (2012). +Optimal stability polynomials for numerical integration of initial value problems +[DOI: 10.2140/camcos.2012.7.247](https://doi.org/10.2140/camcos.2012.7.247) +=# + +# Perform bisection to optimize timestep for stability of the polynomial +function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, + num_stage_evals, + dtmax, dteps, eig_vals; + verbose = false) + dtmin = 0.0 + dt = -1.0 + abs_p = -1.0 + + # Construct stability polynomial for each eigenvalue + pnoms = ones(Complex{Float64}, num_eig_vals, 1) + + # Init datastructure for monomial coefficients + gamma = Variable(num_stage_evals - consistency_order) + + normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) + + for j in 1:num_stage_evals + fac_j = factorial(j) + for i in 1:num_eig_vals + normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j + end + end + + normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) + + if verbose + println("Start optimization of stability polynomial \n") + end + + # Bisection on timestep + while dtmax - dtmin > dteps + dt = 0.5 * (dtmax + dtmin) + + # Compute stability polynomial for current timestep + for k in 1:num_stage_evals + dt_k = dt^k + for i in 1:num_eig_vals + normalized_powered_eigvals_scaled[i, k] = dt_k * + normalized_powered_eigvals[i, + k] + end + end + + # Use last optimal values for gamma in (potentially) next iteration + problem = minimize(stability_polynomials!(pnoms, consistency_order, + num_stage_evals, + normalized_powered_eigvals_scaled, + gamma)) + + solve!(problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + "feastol" => 1e-9, + "abstol" => 1e-9, + "reltol" => 1e-9, + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent_solver = true) + + abs_p = problem.optval + + if abs_p < 1 + dtmin = dt + else + dtmax = dt + end + end + + if verbose + println("Concluded stability polynomial optimization \n") + end + + return evaluate(gamma), dt +end +end # @muladd + +end # module TrixiConvexECOSExt diff --git a/src/Trixi.jl b/src/Trixi.jl index f3977f1f058..3a882d0962c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -310,6 +310,14 @@ function __init__() end end + @static if !isdefined(Base, :get_extension) + @require Convex="f65535da-76fb-5f13-bab9-19810c17039a" begin + @require ECOS="e2685f51-7e38-5353-a97d-a921fd2c8199" begin + include("../ext/TrixiConvexECOSExt.jl") + end + 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 new file mode 100644 index 00000000000..b3b917dc18d --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -0,0 +1,410 @@ +# 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 +using LinearAlgebra: eigvals + +@muladd begin +#! format: noindent + +# Abstract base type for both single/standalone and multi-level +# PERK (Paired-Explicit Runge-Kutta) time integration schemes +abstract type AbstractPairedExplicitRK end +# Abstract base type for single/standalone PERK time integration schemes +abstract type AbstractPairedExplicitRKSingle <: AbstractPairedExplicitRK end + +# Compute the coefficients of the A matrix in the Butcher tableau using +# stage scaling factors and monomial coefficients +function compute_a_coeffs(num_stage_evals, stage_scaling_factors, monomial_coeffs) + a_coeffs = copy(monomial_coeffs) + + for stage in 1:(num_stage_evals - 2) + a_coeffs[stage] /= stage_scaling_factors[stage] + for prev_stage in 1:(stage - 1) + a_coeffs[stage] /= a_coeffs[prev_stage] + end + end + + return reverse(a_coeffs) +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2 +# using a list of eigenvalues +function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, + bS, cS; verbose = false) + # c Vector from Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages + c[k] = cS * (k - 1) / (num_stages - 1) + end + stage_scaling_factors = bS * reverse(c[2:(end - 1)]) + + # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) + coeffs_max = num_stages - 2 + + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] + + consistency_order = 2 + + dtmax = tspan[2] - tspan[1] + dteps = 1e-9 # Hyperparameter of the optimization, might be too large for systems requiring very small timesteps + + 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) + + num_monomial_coeffs = length(monomial_coeffs) + @assert num_monomial_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + return a_matrix, c +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2 +# using provided monomial coefficients file +function compute_PairedExplicitRK2_butcher_tableau(num_stages, + base_path_monomial_coeffs::AbstractString, + bS, cS) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages + c[k] = cS * (k - 1) / (num_stages - 1) + end + stage_scaling_factors = bS * reverse(c[2:(end - 1)]) + + # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) + coeffs_max = num_stages - 2 + + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] + + path_monomial_coeffs = joinpath(base_path_monomial_coeffs, + "gamma_" * string(num_stages) * ".txt") + + @assert isfile(path_monomial_coeffs) "Couldn't find file" + monomial_coeffs = readdlm(path_monomial_coeffs, Float64) + num_monomial_coeffs = size(monomial_coeffs, 1) + + @assert num_monomial_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + return a_matrix, c +end + +@doc raw""" + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + bS = 1.0, cS = 0.5) + PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, bS = 1.0, cS = 0.5) + PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, bS = 1.0, cS = 0.5) + Parameters: + - `num_stages` (`Int`): Number of stages in the PERK method. + - `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing + monomial coefficients of the stability polynomial of PERK method. + The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks. + - `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. + - `bS` (`Float64`, optional): Value of b in the Butcher tableau at b_s, when + s is the number of stages, default is 1.0. + - `cS` (`Float64`, optional): Value of c in the Butcher tableau at c_s, when + s is the number of stages, default is 0.5. + +The following structures and methods provide a minimal implementation of +the second-order paired explicit Runge-Kutta (PERK) method +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +- Brian Vermeire (2019). + Paired explicit Runge-Kutta schemes for stiff systems of equations + [DOI: 10.1016/j.jcp.2019.05.014](https://doi.org/10.1016/j.jcp.2019.05.014) +""" +mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle + const num_stages::Int + + a_matrix::Matrix{Float64} + c::Vector{Float64} + b1::Float64 + bS::Float64 + cS::Float64 +end # struct PairedExplicitRK2 + +# Constructor that reads the coefficients from a file +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + bS = 1.0, cS = 0.5) + a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, + base_path_monomial_coeffs, + bS, cS) + + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) +end + +# Constructor that calculates the coefficients with polynomial optimizer from a +# semidiscretization +function PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, + bS = 1.0, cS = 0.5) + eig_vals = eigvals(jacobian_ad_forward(semi)) + + return PairedExplicitRK2(num_stages, tspan, eig_vals; verbose, bS, cS) +end + +# Constructor that calculates the coefficients with polynomial optimizer from a +# list of eigenvalues +function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, + bS = 1.0, cS = 0.5) + a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, + eig_vals, tspan, + bS, cS; + verbose) + + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) +end + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 +mutable struct PairedExplicitRKOptions{Callback} + callback::Callback # callbacks; used in Trixi + adaptive::Bool # whether the algorithm is adaptive; ignored + dtmax::Float64 # ignored + maxiters::Int # maximal number of time steps + tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored +end + +function PairedExplicitRKOptions(callback, tspan; maxiters = typemax(Int), kwargs...) + PairedExplicitRKOptions{typeof(callback)}(callback, false, Inf, maxiters, + [last(tspan)]) +end + +abstract type PairedExplicitRK end +abstract type AbstractPairedExplicitRKSingleIntegrator <: PairedExplicitRK 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. +mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, Alg, + PairedExplicitRKOptions} <: + AbstractPairedExplicitRKSingleIntegrator + u::uType + du::uType + u_tmp::uType + t::RealT + dt::RealT # current time step + dtcache::RealT # ignored + 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 + # PairedExplicitRK2 stages: + k1::uType + k_higher::uType +end + +# Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) +function Base.getproperty(integrator::PairedExplicitRK, field::Symbol) + if field === :stats + return (naccept = getfield(integrator, :iter),) + end + # general fallback + return getfield(integrator, field) +end + +function init(ode::ODEProblem, alg::PairedExplicitRK2; + dt, callback = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(u0) + u_tmp = zero(u0) + + # PairedExplicitRK2 stages + k1 = zero(u0) + k_higher = zero(u0) + + t0 = first(ode.tspan) + iter = 0 + + integrator = PairedExplicitRK2Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, + ode.p, + (prob = ode,), ode.f, alg, + PairedExplicitRKOptions(callback, + ode.tspan; + kwargs...), + false, + k1, k_higher) + + # initialize callbacks + if callback isa CallbackSet + for cb in callback.continuous_callbacks + error("unsupported") + 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::PairedExplicitRK2; + dt, callback = nothing, kwargs...) + integrator = init(ode, alg, dt = dt, callback = callback; kwargs...) + + # Start actual solve + solve_steps!(integrator) +end + +function solve_steps!(integrator::PairedExplicitRK2Integrator) + @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::PairedExplicitRK2Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + integrator.finalstep = false + + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end + + # 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.u) + 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) + # Construct current state + @threaded for i in eachindex(integrator.u) + 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 + + @threaded for i in eachindex(integrator.u) + integrator.u[i] += alg.b1 * integrator.k1[i] + + alg.bS * integrator.k_higher[i] + end + end # PairedExplicitRK2 step + + integrator.iter += 1 + integrator.t += integrator.dt + + @trixi_timeit timer() "Step-Callbacks" begin + # handle callbacks + if callbacks isa CallbackSet + foreach(callbacks.discrete_callbacks) do cb + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) + end + return nothing + 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 + +# get a cache where the RHS can be stored +get_du(integrator::PairedExplicitRK) = integrator.du +get_tmp_cache(integrator::PairedExplicitRK) = (integrator.u_tmp,) + +# some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u +u_modified!(integrator::PairedExplicitRK, ::Bool) = false + +# used by adaptive timestepping algorithms in DiffEq +function set_proposed_dt!(integrator::PairedExplicitRK, dt) + integrator.dt = dt +end + +function get_proposed_dt(integrator::PairedExplicitRK) + return integrator.dt +end + +# stop the time integration +function terminate!(integrator::PairedExplicitRK) + integrator.finalstep = true + empty!(integrator.opts.tstops) +end + +# used for AMR (Adaptive Mesh Refinement) +function Base.resize!(integrator::PairedExplicitRK2Integrator, 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 new file mode 100644 index 00000000000..b73ea758312 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -0,0 +1,12 @@ +# 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 + +# Basic implementation of the second-order paired explicit Runge-Kutta (PERK) method +include("methods_PERK2.jl") +# Define all of the functions necessary for polynomial optimizations +include("polynomial_optimizer.jl") +end # @muladd diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl new file mode 100644 index 00000000000..bfd53ba2eaf --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -0,0 +1,28 @@ +# Filter out eigenvalues with positive real parts, those with negative imaginary +# parts due to eigenvalues' symmetry around the real axis, or the eigenvalues +# that are smaller than a specified threshold. +function filter_eig_vals(eig_vals, threshold = 1e-12; verbose = false) + filtered_eig_vals = Complex{Float64}[] + + for eig_val in eig_vals + if real(eig_val) < 0 && imag(eig_val) > 0 && abs(eig_val) >= threshold + push!(filtered_eig_vals, eig_val) + end + end + + filtered_eig_vals_count = length(eig_vals) - length(filtered_eig_vals) + + if verbose + println("$filtered_eig_vals_count eigenvalue(s) are not passed on because " * + "they either are in magnitude smaller than $threshold, have positive " * + "real parts, or have negative imaginary parts.\n") + end + + return length(filtered_eig_vals), filtered_eig_vals +end + +# Add definitions of functions related to polynomial optimization by Convex and ECOS here +# such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package +# extension or by the Convex and ECOS-specific code loaded by Requires.jl +function undo_normalization! end +function bisect_stability_polynomial end diff --git a/src/time_integration/time_integration.jl b/src/time_integration/time_integration.jl index c1e53527121..d19a1fcc37c 100644 --- a/src/time_integration/time_integration.jl +++ b/src/time_integration/time_integration.jl @@ -16,4 +16,5 @@ end include("methods_2N.jl") include("methods_3Sstar.jl") include("methods_SSP.jl") +include("paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl") end # @muladd diff --git a/test/Project.toml b/test/Project.toml index 1491d7a5c5f..5fc2bb18bdf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,10 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Convex = "f65535da-76fb-5f13-bab9-19810c17039a" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -16,7 +19,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" +Convex = "0.15.4" +DelimitedFiles = "1" Downloads = "1" +ECOS = "1.1.2" ExplicitImports = "1.0.1" FFMPEG = "0.4" ForwardDiff = "0.10.24" diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index a580a3b5600..afa92efeddb 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -81,6 +81,20 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_advection_perk2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), + l2=[0.014139244532882265], + linf=[0.019997568971592217]) + # 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 end end # module diff --git a/test/test_trixi.jl b/test/test_trixi.jl index 78195825886..9bfd73ea28d 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -155,7 +155,10 @@ macro test_nowarn_mod(expr, additional_ignore_content = String[]) # TODO: Silence warning introduced by Flux v0.13.13. Should be properly fixed. r"┌ Warning: Layer with Float32 parameters got Float64 input.+\n│.+\n│.+\n│.+\n└ @ Flux.+\n", # NOTE: These warnings arose from Julia 1.10 onwards - r"WARNING: Method definition .* in module .* at .* overwritten .*.\n"] + r"WARNING: Method definition .* in module .* at .* overwritten .*.\n", + # Warnings from third party packages + r"┌ Warning: Problem status ALMOST_INFEASIBLE; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n", + r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n"] append!(ignore_content, $additional_ignore_content) for pattern in ignore_content stderr_content = replace(stderr_content, pattern => "") diff --git a/test/test_unit.jl b/test/test_unit.jl index 90ee21030d3..de13d41e931 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -3,6 +3,13 @@ module TestUnit using Test using Trixi +using DelimitedFiles: readdlm + +# Use Convex and ECOS to load the extension that extends functions for testing +# PERK Single p2 Constructors +using Convex: Convex +using ECOS: Optimizer + include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists @@ -1634,6 +1641,39 @@ end @test mesh.boundary_faces[:entire_boundary] == [1, 2] end +@testset "PERK Single p2 Constructors" begin + path_coeff_file = mktempdir() + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", + joinpath(path_coeff_file, "gamma_6.txt")) + + ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file) + + @test isapprox(ode_algorithm.a_matrix, + [0.12405417889682908 0.07594582110317093 + 0.16178873711001726 0.13821126288998273 + 0.16692313960864164 0.2330768603913584 + 0.12281292901258256 0.37718707098741744], atol = 1e-13) + + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/c7a89eaaa857e87dde055f78eae9b94a/raw/2937f8872ffdc08e0dcf444ee35f9ebfe18735b0/Spectrum_2D_IsentropicVortex_CEE.txt", + joinpath(path_coeff_file, "spectrum_2d.txt")) + + eig_vals = readdlm(joinpath(path_coeff_file, "spectrum_2d.txt"), ComplexF64) + tspan = (0.0, 1.0) + ode_algorithm = Trixi.PairedExplicitRK2(12, tspan, vec(eig_vals)) + + @test isapprox(ode_algorithm.a_matrix, + [0.06453812656705388 0.02637096434203703 + 0.09470601372266194 0.04165762264097442 + 0.12332877820057538 0.05848940361760645 + 0.1498701503275483 0.07740257694517898 + 0.173421149536068 0.09930612319120471 + 0.19261978147927503 0.12556203670254315 + 0.2052334022622969 0.15840296137406676 + 0.2073489042901963 0.2017420048007128 + 0.19135142349998963 0.2631940310454649 + 0.13942836392940833 0.3605716360705917], atol = 1e-13) +end + @testset "Sutherlands Law" begin function mu(u, equations) T_ref = 291.15