From befd590bfa93576601ccab569bbb62f534dd30af Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 29 Jan 2024 10:06:40 +0100 Subject: [PATCH 01/64] Templates for PERK p2, p3 --- src/Trixi.jl | 2 + src/auxiliary/auxiliary.jl | 14 + .../methods_PERK2.jl | 308 ++++++++++++++++++ .../methods_PERK3.jl | 248 ++++++++++++++ .../paired_explicit_runge_kutta.jl | 11 + src/time_integration/time_integration.jl | 1 + 6 files changed, 584 insertions(+) create mode 100644 src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl create mode 100644 src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl create mode 100644 src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index e18b2f6415c..33f014fb9f8 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -274,6 +274,8 @@ export trixi_include, examples_dir, get_examples, default_example, export ode_norm, ode_unstable_check +export PERK, PERK3 + export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 1f7d30d6aa8..1ee36238276 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -345,4 +345,18 @@ function register_error_hints() return nothing end + +function read_file(FilePath::AbstractString, DataType::Type=Float64) + @assert isfile(FilePath) "Couldn't find file" + Data = zeros(DataType, 0) + open(FilePath, "r") do File + while !eof(File) + LineContent = readline(File) + append!(Data, parse(DataType, LineContent)) + end + end + NumLines = length(Data) + + return NumLines, Data + end end # @muladd 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..e57807e9238 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -0,0 +1,308 @@ +# 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 + +abstract type PERK end +# Abstract base type for single/standalone P-ERK time integration schemes +abstract type PERKSingle <: PERK end + +function ComputeACoeffs(NumStageEvals::Int, + SE_Factors::Vector{Float64}, MonCoeffs::Vector{Float64}) + ACoeffs = MonCoeffs + + for stage in 1:NumStageEvals - 2 + ACoeffs[stage] /= SE_Factors[stage] + for prev_stage in 1:stage-1 + ACoeffs[stage] /= ACoeffs[prev_stage] + end + end + + return reverse(ACoeffs) +end + +function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, bS::Float64, cEnd::Float64) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(NumStages) + for k in 2:NumStages + c[k] = cEnd * (k - 1)/(NumStages - 1) + end + #= + for k in 2:NumStages + c[k] = (k - 1)/(2.0*(NumStages - 1)) + end + =# + println("Timestep-split: "); display(c); println("\n") + SE_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) + CoeffsMax = NumStages - 2 + + AMatrix = zeros(CoeffsMax, 2) + AMatrix[:, 1] = c[3:end] + + + PathMonCoeffs = BasePathMonCoeffs * "gamma_" * string(NumStages) * ".txt" + NumMonCoeffs, MonCoeffs = read_file(PathMonCoeffs, Float64) + @assert NumMonCoeffs == CoeffsMax + A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) + + + #= + # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) + PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * ".txt" + NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) + @assert NumMonCoeffs == CoeffsMax + =# + + AMatrix[:, 1] -= A + AMatrix[:, 2] = A + + println("A matrix: "); display(AMatrix); println() + + return AMatrix, c +end + +""" + PERK2() + +The following structures and methods provide a minimal implementation of +the paired explicit Runge-Kutta method (https://doi.org/10.1016/j.jcp.2019.05.014) +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +This is using the same interface as OrdinaryDiffEq.jl, copied from file "methods_2N.jl" for the +CarpenterKennedy2N{54, 43} methods. +""" + +mutable struct PERK2 <: PERKSingle + const NumStages::Int + + AMatrix::Matrix{Float64} + c::Vector{Float64} + bS::Float64 + b1::Float64 + cEnd::Float64 + + # Constructor for previously computed A Coeffs + function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64, cEnd_::Float64) + + newPERK2 = new(NumStages_) + + newPERK2.AMatrix, newPERK2.c = + ComputePERK2_ButcherTableau(NumStages_, BasePathMonCoeffs_, bS_, cEnd_) + + newPERK2.b1 = one(bS_) - bS_ + newPERK2.bS = bS_ + newPERK2.cEnd = cEnd_ + return newPERK2 + end +end # struct PERK2 + + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 +mutable struct PERK_IntegratorOptions{Callback} + callback::Callback # callbacks; used in Trixi + adaptive::Bool # whether the algorithm is adaptive; ignored + dtmax::Float64 # ignored + maxiters::Int # maximal numer of time steps + tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored +end + +function PERK_IntegratorOptions(callback, tspan; maxiters=typemax(Int), kwargs...) + PERK_IntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, [last(tspan)]) +end + +abstract type PERK_Integrator end +abstract type PERKSingle_Integrator <: PERK_Integrator 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 PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK2_IntegratorOptions} <: PERKSingle_Integrator + 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::PERK_IntegratorOptions + finalstep::Bool # added for convenience + # PERK2 stages: + k1::uType + k_higher::uType + t_stage::RealT + du_ode_hyp::uType # TODO: Not best solution since this is not needed for hyperbolic problems +end + +# Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) +function Base.getproperty(integrator::PERK_Integrator, field::Symbol) + if field === :stats + return (naccept = getfield(integrator, :iter),) + end + # general fallback + return getfield(integrator, field) +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::PERK2; + dt, callback=nothing, kwargs...) + + u0 = copy(ode.u0) + du = zero(u0) #previously: similar(u0) + u_tmp = zero(u0) + + # PERK2 stages + k1 = zero(u0) + k_higher = zero(u0) + + t0 = first(ode.tspan) + iter = 0 + + integrator = PERK2_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, + (prob=ode,), ode.f, alg, + PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false, + k1, k_higher, t0, du_ode_hyp) + + # 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 + + solve!(integrator) +end + +function solve!(integrator::PERK2_Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + integrator.finalstep = false + + @trixi_timeit timer() "main loop" while !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 + + # k2 + integrator.t_stage = integrator.t + alg.c[2] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + # Higher stages + for stage = 3:alg.NumStages + integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.AMatrix[stage - 2, 1] * integrator.k1[i] + + alg.AMatrix[stage - 2, 2] * integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @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] += integrator.k_higher[i] + integrator.u[i] += alg.b1 * integrator.k1[i] + alg.bS * integrator.k_higher[i] + end + end # PERK2 step + + 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 # "main loop" timer + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) +end + +# get a cache where the RHS can be stored +get_du(integrator::PERK_Integrator) = integrator.du +get_tmp_cache(integrator::PERK_Integrator) = (integrator.u_tmp,) + +# some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u +u_modified!(integrator::PERK_Integrator, ::Bool) = false + +# used by adaptive timestepping algorithms in DiffEq +function set_proposed_dt!(integrator::PERK_Integrator, dt) + integrator.dt = dt +end + +function get_proposed_dt(integrator::PERK_Integrator) + return integrator.dt +end + +# stop the time integration +function terminate!(integrator::PERK_Integrator) + integrator.finalstep = true + empty!(integrator.opts.tstops) +end + +# used for AMR (Adaptive Mesh Refinement) +function Base.resize!(integrator::PERK2_Integrator, 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 \ No newline at end of file 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 00000000000..61cefbe8ac0 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -0,0 +1,248 @@ +# 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 + +function ComputePERK3_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, cS2::Float64) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(NumStages) + for k in 2:NumStages-2 + c[k] = cS2 * (k - 1)/(NumStages - 3) # Equidistant timestep distribution (similar to PERK2) + end + + # Original third order proposed PERK + #= + c[NumStages - 1] = 1.0/3.0 + c[NumStages] = 1.0 + =# + # Own third order PERK based on SSPRK33 + c[NumStages - 1] = 1.0 + c[NumStages] = 0.5 + + println("Timestep-split: "); display(c); println("\n") + + # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) + CoeffsMax = NumStages - 2 + + AMatrix = zeros(CoeffsMax, 2) + AMatrix[:, 1] = c[3:end] + + PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * "_" * string(NumStages) * ".txt" + NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) + @assert NumMonCoeffs == CoeffsMax + + AMatrix[:, 1] -= A + AMatrix[:, 2] = A + + println("A matrix: "); display(AMatrix); println() + + return AMatrix, c +end + +""" + PERK3() + +The following structures and methods provide a minimal implementation of +the third order paired explicit Runge-Kutta method (https://www.sciencedirect.com/science/article/pii/S0021999122005320) +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +This is using the same interface as OrdinaryDiffEq.jl, copied from file "methods_2N.jl" for the +CarpenterKennedy2N{54, 43} methods. +""" + +mutable struct PERK3 <: PERKSingle + const NumStages::Int + + AMatrix::Matrix{Float64} + c::Vector{Float64} + + # Constructor for previously computed A Coeffs + function PERK3(NumStages_::Int, BasePathMonCoeffs_::AbstractString, cS2_::Float64=1.0) + + newPERK3 = new(NumStages_) + + newPERK3.AMatrix, newPERK3.c = + ComputePERK3_ButcherTableau(NumStages_, BasePathMonCoeffs_, cS2_) + return newPERK3 + end +end # struct PERK3 + + +# 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 PERK3_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK_IntegratorOptions} + 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::PERK_IntegratorOptions + finalstep::Bool # added for convenience + # PERK stages: + k1::uType + k_higher::uType + k_S1::uType # Required for third order + t_stage::RealT +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::PERK3; + dt, callback=nothing, kwargs...) + + u0 = copy(ode.u0) + du = zero(u0) #previously: similar(u0) + u_tmp = zero(u0) + + # PERK stages + k1 = zero(u0) + k_higher = zero(u0) + k_S1 = zero(u0) + + t0 = first(ode.tspan) + iter = 0 + + integrator = PERK3_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, + (prob=ode,), ode.f, alg, + PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false, + k1, k_higher, k_S1, t0) + + # 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 + + solve!(integrator) +end + +function solve!(integrator::PERK3_Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + integrator.finalstep = false + + #@trixi_timeit timer() "main loop" while !integrator.finalstep + while !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 + + # k2 + integrator.t_stage = integrator.t + alg.c[2] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + if alg.NumStages == 3 + @threaded for i in eachindex(integrator.du) + integrator.k_S1[i] = integrator.k_higher[i] + end + end + + # Higher stages + for stage = 3:alg.NumStages + integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.AMatrix[stage - 2, 1] * integrator.k1[i] + + alg.AMatrix[stage - 2, 2] * integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + # IDEA: Stop for loop at NumStages -1 to avoid if (maybe more performant?) + if stage == alg.NumStages - 1 + @threaded for i in eachindex(integrator.du) + integrator.k_S1[i] = integrator.k_higher[i] + end + end + end + + @threaded for i in eachindex(integrator.u) + # Original proposed PERK3 + #integrator.u[i] += 0.75 * integrator.k_S1[i] + 0.25 * integrator.k_higher[i] + # Own PERK based on SSPRK33 + integrator.u[i] += (integrator.k1[i] + integrator.k_S1[i] + 4.0 * integrator.k_higher[i])/6.0 + end + end # PERK 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 # "main loop" timer + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) +end + +# used for AMR (Adaptive Mesh Refinement) +function Base.resize!(integrator::PERK3_Integrator, 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) + resize!(integrator.k_S1, new_size) +end + +end # @muladd \ No newline at end of file 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..03689b75321 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -0,0 +1,11 @@ +# 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 + +include("methods_PERK2.jl") +include("methods_PERK3.jl") + +end # @muladd \ No newline at end of file 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 From 0b378f79a789b85fa51cc101e907972cd118b06d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 29 Jan 2024 10:45:50 +0100 Subject: [PATCH 02/64] example --- .../tree_1d_dgsem/elixir_advection_PERK2.jl | 66 +++++++++++++++++++ src/Trixi.jl | 2 +- .../methods_PERK2.jl | 12 +--- 3 files changed, 70 insertions(+), 10 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_advection_PERK2.jl 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..07aa047d751 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -0,0 +1,66 @@ + +using OrdinaryDiffEq +using Trixi +using Downloads: download + +############################################################################### +# 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 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# 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_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 2.5) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +MonCoeffsFile = "gamma_6.txt" + +download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", +MonCoeffsFile) + +ode_algorithm = PERK2(6, "./") + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +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() + +using Plots +plot(sol) \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index 33f014fb9f8..0e96bb8dfd8 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -274,7 +274,7 @@ export trixi_include, examples_dir, get_examples, default_example, export ode_norm, ode_unstable_check -export PERK, PERK3 +export PERK2, PERK3 export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure 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 e57807e9238..3d247fd65e4 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -29,11 +29,6 @@ function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::Abstract for k in 2:NumStages c[k] = cEnd * (k - 1)/(NumStages - 1) end - #= - for k in 2:NumStages - c[k] = (k - 1)/(2.0*(NumStages - 1)) - end - =# println("Timestep-split: "); display(c); println("\n") SE_Factors = bS * reverse(c[2:end-1]) @@ -86,7 +81,7 @@ mutable struct PERK2 <: PERKSingle cEnd::Float64 # Constructor for previously computed A Coeffs - function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64, cEnd_::Float64) + function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5) newPERK2 = new(NumStages_) @@ -121,7 +116,7 @@ abstract type PERKSingle_Integrator <: PERK_Integrator end # 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 PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK2_IntegratorOptions} <: PERKSingle_Integrator +mutable struct PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK_IntegratorOptions} <: PERKSingle_Integrator u::uType du::uType u_tmp::uType @@ -139,7 +134,6 @@ mutable struct PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK2_I k1::uType k_higher::uType t_stage::RealT - du_ode_hyp::uType # TODO: Not best solution since this is not needed for hyperbolic problems end # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) @@ -169,7 +163,7 @@ function solve(ode::ODEProblem, alg::PERK2; integrator = PERK2_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, (prob=ode,), ode.f, alg, PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher, t0, du_ode_hyp) + k1, k_higher, t0) # initialize callbacks if callback isa CallbackSet From 42c74d1954eb1c2ad7e8c1beaebc6fafa4e01051 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 2 Feb 2024 10:00:27 +0100 Subject: [PATCH 03/64] constructor changed --- .../tree_1d_dgsem/elixir_advection_PERK2.jl | 7 +- .../methods_PERK2.jl | 23 ++- .../polynomial_optimizer.jl | 147 ++++++++++++++++++ 3 files changed, 166 insertions(+), 11 deletions(-) create mode 100644 src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 07aa047d751..7c22fed20c0 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -47,12 +47,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation -MonCoeffsFile = "gamma_6.txt" - -download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", -MonCoeffsFile) - -ode_algorithm = PERK2(6, "./") +ode_algorithm = PERK2(6, semi) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = Trixi.solve(ode, ode_algorithm, 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 3d247fd65e4..b13b8204f6b 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,6 +2,7 @@ # 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. +include("polynomial_optimizer.jl") @muladd begin abstract type PERK end @@ -22,7 +23,7 @@ function ComputeACoeffs(NumStageEvals::Int, return reverse(ACoeffs) end -function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, bS::Float64, cEnd::Float64) +function ComputePERK2_ButcherTableau(NumStages::Int, semi::AbstractSemidiscretization, bS::Float64, cEnd::Float64) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(NumStages) @@ -39,8 +40,20 @@ function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::Abstract AMatrix[:, 1] = c[3:end] - PathMonCoeffs = BasePathMonCoeffs * "gamma_" * string(NumStages) * ".txt" - NumMonCoeffs, MonCoeffs = read_file(PathMonCoeffs, Float64) + ConsOrder = 2 + filter_thres = 1e-12 + dtMax = 1.0 + dtEps = 1e-9 + + J = jacobian_ad_forward(semi) + EigVals = eigvals(J) + + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + + MonCoeffs, PWorstCase, dtOpt = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + + NumMonCoeffs = length(MonCoeffs) @assert NumMonCoeffs == CoeffsMax A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) @@ -81,12 +94,12 @@ mutable struct PERK2 <: PERKSingle cEnd::Float64 # Constructor for previously computed A Coeffs - function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5) + function PERK2(NumStages_::Int, semi_::AbstractSemidiscretization, bS_::Float64=1.0, cEnd_::Float64=0.5) newPERK2 = new(NumStages_) newPERK2.AMatrix, newPERK2.c = - ComputePERK2_ButcherTableau(NumStages_, BasePathMonCoeffs_, bS_, cEnd_) + ComputePERK2_ButcherTableau(NumStages_, semi_, bS_, cEnd_) newPERK2.b1 = one(bS_) - bS_ newPERK2.bS = bS_ 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..9bf795d2b3e --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -0,0 +1,147 @@ +using LinearAlgebra +using Convex +using ECOS +const MOI = Convex.MOI + +function filter_Eigvals(EigVals::Array{Complex{Float64}}, threshold:: Float64) + filtered_eigvals_counter = 0 + #FilteredEigVals \is not an Array but the code seems to work fine regardless + FilteredEigVals = Complex{Float64}[] + for EigVal in EigVals + if abs(EigVal) < threshold + filtered_eigvals_counter += 1 + else + push!(FilteredEigVals, EigVal) + end + end + + println("Total number of $filtered_eigvals_counter eigenvalues has not been recorded because they were smaller than $threshold") + + return length(FilteredEigVals), FilteredEigVals +end + +function ReadInEigVals(Path_toEvalFile::AbstractString) + + NumEigVals = -1 # Declare and set to some value + open(Path_toEvalFile, "r") do EvalFile + NumEigVals = countlines(EvalFile) + end + + EigVals = Array{Complex{Float64}}(undef, NumEigVals) + + LineIndex = 0 + open(Path_toEvalFile, "r") do EvalFile + # read till end of file + while !eof(EvalFile) + + # read a new / next line for every iteration + LineContent = readline(EvalFile) + + EigVals[LineIndex + 1] = parse(Complex{Float64}, LineContent) + + LineIndex += 1 + end + end + + return NumEigVals, EigVals +end + +function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) + + for i in 1:NumEigVals + pnoms[i] = 1.0 + end + + for k in 1:ConsOrder + for i in 1:NumEigVals + pnoms[i] += normalized_powered_EigvalsScaled[i,k] + end + #pnoms += 1/factorial(k) * EigValsScaled.^k + end + + for k in ConsOrder + 1:NumStages + #Attemp to use nested loop instead of [:,k] has been made. Convex.jl doesn't accept it. + pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] + end + + return maximum(abs(pnoms)) +end + + +function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) + #dtMin = Base.big(0) + dtMin = 0.0 + + #dt = Base.big(-1.0) + dt = -1.0 + + #AbsP = Base.big(-1.0) + AbsP = -1.0 + + pnoms = ones(Complex{Float64}, NumEigVals, 1) + + gamma = Variable(NumStages - ConsOrder) # Init datastructure for results + + normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStages) + + for j in 1:NumStages + fac_j = factorial(j) + for i in 1:NumEigVals + normalized_powered_Eigvals[i, j] = EigVals[i]^j / fac_j + end + end + + normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStages) + + while dtMax - dtMin > dtEps + #dt = Base.big(0.5) * (dtMax + dtMin) + dt = 0.5 * (dtMax + dtMin) + + for k in 1:NumStages + dt_k = dt^k + for i in 1:NumEigVals + normalized_powered_EigvalsScaled[i,k] = dt_k * normalized_powered_Eigvals[i,k] + end + end + + # Use last optimal values for gamm0 in (potentially) next iteration + problem = minimize(Polynoms(ConsOrder, NumStages, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) + + Convex.solve!( + problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + #"eps" => 1e9, # 1e-13 + "feastol" => 1e-9, # 1e-9 + "abstol" => 1e-9, # 1e-9 + "reltol" => 1e-9, # 1e-9 + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent_solver = false + ) + + AbsP = problem.optval + + println("Current MaxAbsP: ", AbsP, "\nCurrent dt: ", dt, "\n") + + #if AbsP < Base.big(1.0) + if AbsP < 1.0 + dtMin = dt + else + dtMax = dt + end + end + + return evaluate(gamma), AbsP, dt +end + +function undo_normalization(ConsOrder::Int, NumStages::Int, gammaOpt) + for k in ConsOrder + 1:NumStages + gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) + end + return gammaOpt +end \ No newline at end of file From 04784df2bc3aace52f346e6be5fc47f06ac38f63 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 2 Feb 2024 11:26:06 +0100 Subject: [PATCH 04/64] Modified so that both of the constructor stay --- .../tree_1d_dgsem/elixir_advection_PERK2.jl | 7 +++ .../methods_PERK2.jl | 54 ++++++++++++++++++- .../polynomial_optimizer.jl | 10 +--- 3 files changed, 62 insertions(+), 9 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 7c22fed20c0..84ded3528b8 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -47,6 +47,13 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation +MonCoeffsFile = "gamma_6.txt" + +download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", +MonCoeffsFile) + +#ode_algorithm = PERK2(6, "./") + ode_algorithm = PERK2(6, semi) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks 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 b13b8204f6b..11cbeba7154 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -73,6 +73,44 @@ function ComputePERK2_ButcherTableau(NumStages::Int, semi::AbstractSemidiscretiz return AMatrix, c end +function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, bS::Float64, cEnd::Float64) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(NumStages) + for k in 2:NumStages + c[k] = cEnd * (k - 1)/(NumStages - 1) + end + println("Timestep-split: "); display(c); println("\n") + SE_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) + CoeffsMax = NumStages - 2 + + AMatrix = zeros(CoeffsMax, 2) + AMatrix[:, 1] = c[3:end] + + + PathMonCoeffs = BasePathMonCoeffs * "gamma_" * string(NumStages) * ".txt" + NumMonCoeffs, MonCoeffs = read_file(PathMonCoeffs, Float64) + @assert NumMonCoeffs == CoeffsMax + A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) + + + #= + # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) + PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * ".txt" + NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) + @assert NumMonCoeffs == CoeffsMax + =# + + AMatrix[:, 1] -= A + AMatrix[:, 2] = A + + println("A matrix: "); display(AMatrix); println() + + return AMatrix, c +end + """ PERK2() @@ -93,7 +131,21 @@ mutable struct PERK2 <: PERKSingle b1::Float64 cEnd::Float64 - # Constructor for previously computed A Coeffs + #Constructor that read the coefficients from the file + function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5) + + newPERK2 = new(NumStages_) + + newPERK2.AMatrix, newPERK2.c = + ComputePERK2_ButcherTableau(NumStages_, BasePathMonCoeffs_, bS_, cEnd_) + + newPERK2.b1 = one(bS_) - bS_ + newPERK2.bS = bS_ + newPERK2.cEnd = cEnd_ + return newPERK2 + end + + #Constructor that calculate the coefficients with polynomial optimizer function PERK2(NumStages_::Int, semi_::AbstractSemidiscretization, bS_::Float64=1.0, cEnd_::Float64=0.5) newPERK2 = new(NumStages_) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 9bf795d2b3e..7c22114dea7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -5,7 +5,7 @@ const MOI = Convex.MOI function filter_Eigvals(EigVals::Array{Complex{Float64}}, threshold:: Float64) filtered_eigvals_counter = 0 - #FilteredEigVals \is not an Array but the code seems to work fine regardless + #FilteredEigVals is not an Array but the code seems to work fine regardless FilteredEigVals = Complex{Float64}[] for EigVal in EigVals if abs(EigVal) < threshold @@ -56,11 +56,9 @@ function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_po for i in 1:NumEigVals pnoms[i] += normalized_powered_EigvalsScaled[i,k] end - #pnoms += 1/factorial(k) * EigValsScaled.^k end for k in ConsOrder + 1:NumStages - #Attemp to use nested loop instead of [:,k] has been made. Convex.jl doesn't accept it. pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] end @@ -69,13 +67,11 @@ end function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) - #dtMin = Base.big(0) + dtMin = 0.0 - #dt = Base.big(-1.0) dt = -1.0 - #AbsP = Base.big(-1.0) AbsP = -1.0 pnoms = ones(Complex{Float64}, NumEigVals, 1) @@ -94,7 +90,6 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStages) while dtMax - dtMin > dtEps - #dt = Base.big(0.5) * (dtMax + dtMin) dt = 0.5 * (dtMax + dtMin) for k in 1:NumStages @@ -128,7 +123,6 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float println("Current MaxAbsP: ", AbsP, "\nCurrent dt: ", dt, "\n") - #if AbsP < Base.big(1.0) if AbsP < 1.0 dtMin = dt else From 55f41653f2706a4db5ff1291eec28e88b5ed4cf9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 09:06:18 +0100 Subject: [PATCH 05/64] bring back eps --- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 7c22114dea7..f0abfcc9ceb 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -107,7 +107,7 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float # Parameters taken from default values for EiCOS MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, "delta" => 2e-7, - #"eps" => 1e9, # 1e-13 + "eps" => 1e-13, # 1e-13 "feastol" => 1e-9, # 1e-9 "abstol" => 1e-9, # 1e-9 "reltol" => 1e-9, # 1e-9 From 6a31a0b3adc80abe4dc421707777883c57a53fe1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 16:06:18 +0100 Subject: [PATCH 06/64] correct val --- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index f0abfcc9ceb..b42049248d5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -107,7 +107,7 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float # Parameters taken from default values for EiCOS MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, "delta" => 2e-7, - "eps" => 1e-13, # 1e-13 + "eps" => 1e13, # 1e13 "feastol" => 1e-9, # 1e-9 "abstol" => 1e-9, # 1e-9 "reltol" => 1e-9, # 1e-9 @@ -116,7 +116,7 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float "reltol_inacc" => 5e-5, "nitref" => 9, "maxit" => 100, - "verbose" => 3); silent_solver = false + "verbose" => 3); silent_solver = true ) AbsP = problem.optval From 87399dc00081ce39ceae78b7fddd1d2176e54581 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 14 Feb 2024 09:57:40 +0100 Subject: [PATCH 07/64] add constructor of PERK3 --- .../methods_PERK2.jl | 1 - .../methods_PERK3.jl | 121 ++++++++++++++++++ .../paired_explicit_runge_kutta.jl | 1 + .../polynomial_optimizer.jl | 22 ++-- 4 files changed, 133 insertions(+), 12 deletions(-) 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 11cbeba7154..f8e061c0872 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,7 +2,6 @@ # 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. -include("polynomial_optimizer.jl") @muladd begin abstract type PERK end diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 61cefbe8ac0..fd3d1705cd6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -2,8 +2,117 @@ # 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 NLsolve @muladd begin +function F_own(aUnknown, NumStages, NumStageEvals, MonCoeffs) + cS2 = 1 # cS2 = c_ts(S-2) + c_ts = c_Own(cS2, NumStages) # ts = timestep + + c_eq = zeros(NumStageEvals - 2) # Add equality constraint that cS2 is equal to 1 + # Both terms should be present + for i in 1:(NumStageEvals - 4) + term1 = aUnknown[NumStageEvals - 1] + term2 = aUnknown[NumStageEvals] + for j in 1:i + term1 *= aUnknown[NumStageEvals - 1 - j] + term2 *= aUnknown[NumStageEvals - j] + end + term1 *= c_ts[NumStages - 2 - i] * 1/6 + term2 *= c_ts[NumStages - 1 - i] * 4/6 + + c_eq[i] = MonCoeffs[i] - (term1 + term2) + end + + # Highest coefficient: Only one term present + i = NumStageEvals - 3 + term2 = aUnknown[NumStageEvals] + for j in 1:i + term2 *= aUnknown[NumStageEvals - j] + end + term2 *= c_ts[NumStages - 1 - i] * 4/6 + + c_eq[i] = MonCoeffs[i] - term2 + + c_eq[NumStageEvals - 2] = 1.0 - 4 * aUnknown[NumStageEvals] - aUnknown[NumStageEvals - 1] + + return c_eq +end + +function c_Own(cS2, NumStages) + c_ts = zeros(NumStages); + + # Last timesteps as for SSPRK33 + c_ts[NumStages] = 0.5; + c_ts[NumStages - 1] = 1; + + # Linear increasing timestep for remainder + for i = 2:NumStages-2 + c_ts[i] = cS2 * (i-1)/(NumStages - 3); + end + return c_ts +end + +function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSemidiscretization) + + #Initialize array of c + c_ts = c_Own(1.0, NumStages) + + #Initialize the array of our solution + aUnknown = zeros(NumStageEvals) + + #Case of e = 3 + if NumStageEvals == 3 + aUnknown = [0, c_ts[2], 0.25] + + else + #Calculate MonCoeffs from polynomial Optimizer + ConsOrder = 3 + dtMax = 1.0 + dtEps = 1e-9 + filter_thres = 1e-12 + + #Get EigVals + J = jacobian_ad_forward(semi) + EigVals = eigvals(J) + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + + MonCoeffs, _, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + + #Define the objective_function + function objective_function(x) + return F_own(x, NumStages, NumStageEvals, MonCoeffs) + end + + #call nlsolver to solve until the answer is not NaN or negative values + is_sol_valid = false + while is_sol_valid == false + #Initialize initial guess + x0 = 0.1 .* rand(NumStageEvals) + x0[1] = Base.big(0) + x0[2] = c_ts[2] + + sol = nlsolve(objective_function, method = :trust_region, x0, + ftol = 1e-15, iterations = 10^4, xtol = 1e-13) + + aUnknown = sol.zero + + is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) + end + end + + println("aUnknown") + println(aUnknown[3:end]) #To debug + + AMatrix = zeros(NumStages - 2, 2) + AMatrix[:, 1] = c_ts[3:end] + AMatrix[:, 1] -= aUnknown[3:end] + AMatrix[:, 2] = aUnknown[3:end] + + return AMatrix, c_ts +end + function ComputePERK3_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, cS2::Float64) # c Vector form Butcher Tableau (defines timestep per stage) @@ -67,6 +176,18 @@ mutable struct PERK3 <: PERKSingle ComputePERK3_ButcherTableau(NumStages_, BasePathMonCoeffs_, cS2_) return newPERK3 end + + #Constructor that compute A coefficients from semidiscretization + function PERK3(NumStages_::Int, NumStageEvals_ ::Int, semi_::AbstractSemidiscretization) + + newPERK3 = new(NumStages_) + + newPERK3.AMatrix, newPERK3.c = + ComputePERK3_ButcherTableau(NumStages_, NumStageEvals_, semi_) + + return newPERK3 + end + end # struct PERK3 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 03689b75321..9bb3d0b817c 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,5 +7,6 @@ include("methods_PERK2.jl") include("methods_PERK3.jl") +include("polynomial_optimizer.jl") end # @muladd \ No newline at end of file diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 7c22114dea7..1ad29cd2f0c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -46,7 +46,7 @@ function ReadInEigVals(Path_toEvalFile::AbstractString) return NumEigVals, EigVals end -function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) +function Polynoms(ConsOrder::Int, NumStageEvals::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) for i in 1:NumEigVals pnoms[i] = 1.0 @@ -58,7 +58,7 @@ function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_po end end - for k in ConsOrder + 1:NumStages + for k in ConsOrder + 1:NumStageEvals pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] end @@ -66,7 +66,7 @@ function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_po end -function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) +function Bisection(ConsOrder::Int, NumEigVals::Int, NumStageEvals::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) dtMin = 0.0 @@ -76,23 +76,23 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float pnoms = ones(Complex{Float64}, NumEigVals, 1) - gamma = Variable(NumStages - ConsOrder) # Init datastructure for results + gamma = Variable(NumStageEvals - ConsOrder) # Init datastructure for results - normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStages) + normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStageEvals) - for j in 1:NumStages + for j in 1:NumStageEvals fac_j = factorial(j) for i in 1:NumEigVals normalized_powered_Eigvals[i, j] = EigVals[i]^j / fac_j end end - normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStages) + normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStageEvals) while dtMax - dtMin > dtEps dt = 0.5 * (dtMax + dtMin) - for k in 1:NumStages + for k in 1:NumStageEvals dt_k = dt^k for i in 1:NumEigVals normalized_powered_EigvalsScaled[i,k] = dt_k * normalized_powered_Eigvals[i,k] @@ -100,7 +100,7 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float end # Use last optimal values for gamm0 in (potentially) next iteration - problem = minimize(Polynoms(ConsOrder, NumStages, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) + problem = minimize(Polynoms(ConsOrder, NumStageEvals, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) Convex.solve!( problem, @@ -133,8 +133,8 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float return evaluate(gamma), AbsP, dt end -function undo_normalization(ConsOrder::Int, NumStages::Int, gammaOpt) - for k in ConsOrder + 1:NumStages +function undo_normalization(ConsOrder::Int, NumStageEvals::Int, gammaOpt) + for k in ConsOrder + 1:NumStageEvals gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) end return gammaOpt From e518386946efe48d9c7d306e19eca4adcec47566 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 14 Feb 2024 16:42:20 +0100 Subject: [PATCH 08/64] minor fixes --- .../paired_explicit_runge_kutta/methods_PERK3.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index fd3d1705cd6..79f090711a8 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -70,14 +70,14 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem ConsOrder = 3 dtMax = 1.0 dtEps = 1e-9 - filter_thres = 1e-12 + filter_threshold = 1e-12 #Get EigVals J = jacobian_ad_forward(semi) EigVals = eigvals(J) - NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_threshold) - MonCoeffs, _, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs, dt, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) #Define the objective_function @@ -90,7 +90,7 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem while is_sol_valid == false #Initialize initial guess x0 = 0.1 .* rand(NumStageEvals) - x0[1] = Base.big(0) + x0[1] = 0.0 x0[2] = c_ts[2] sol = nlsolve(objective_function, method = :trust_region, x0, @@ -98,6 +98,7 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem aUnknown = sol.zero + #First check the values of a[i, i-1] which are acquired by nonlinear solver and c[i] - a[i, i-1]. is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) end end From a1222d7d37bb65e4114eee323e999902807b3df8 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 15 Feb 2024 20:43:22 +0100 Subject: [PATCH 09/64] function and variable name adjustments --- .../methods_PERK3.jl | 205 +++++++++--------- .../polynomial_optimizer.jl | 188 ++++++++-------- 2 files changed, 198 insertions(+), 195 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 79f090711a8..5901d3b5e65 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -5,150 +5,149 @@ using NLsolve @muladd begin -function F_own(aUnknown, NumStages, NumStageEvals, MonCoeffs) - cS2 = 1 # cS2 = c_ts(S-2) - c_ts = c_Own(cS2, NumStages) # ts = timestep +function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) + c_s2 = 1 # c_s2 = c_ts(S-2) + c_ts = c_own(c_s2, num_stages) # ts = timestep - c_eq = zeros(NumStageEvals - 2) # Add equality constraint that cS2 is equal to 1 + c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_s2 is equal to 1 # Both terms should be present - for i in 1:(NumStageEvals - 4) - term1 = aUnknown[NumStageEvals - 1] - term2 = aUnknown[NumStageEvals] + for i in 1:(num_stage_evals - 4) + term1 = a_unknown[num_stage_evals - 1] + term2 = a_unknown[num_stage_evals] for j in 1:i - term1 *= aUnknown[NumStageEvals - 1 - j] - term2 *= aUnknown[NumStageEvals - j] + term1 *= a_unknown[num_stage_evals - 1 - j] + term2 *= a_unknown[num_stage_evals - j] end - term1 *= c_ts[NumStages - 2 - i] * 1/6 - term2 *= c_ts[NumStages - 1 - i] * 4/6 + term1 *= c_ts[num_stages - 2 - i] * 1/6 + term2 *= c_ts[num_stages - 1 - i] * 4/6 - c_eq[i] = MonCoeffs[i] - (term1 + term2) + c_eq[i] = mon_coeffs[i] - (term1 + term2) end # Highest coefficient: Only one term present - i = NumStageEvals - 3 - term2 = aUnknown[NumStageEvals] + i = num_stage_evals - 3 + term2 = a_unknown[num_stage_evals] for j in 1:i - term2 *= aUnknown[NumStageEvals - j] + term2 *= a_unknown[num_stage_evals - j] end - term2 *= c_ts[NumStages - 1 - i] * 4/6 + term2 *= c_ts[num_stages - 1 - i] * 4/6 - c_eq[i] = MonCoeffs[i] - term2 + c_eq[i] = mon_coeffs[i] - term2 - c_eq[NumStageEvals - 2] = 1.0 - 4 * aUnknown[NumStageEvals] - aUnknown[NumStageEvals - 1] + c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - a_unknown[num_stage_evals - 1] return c_eq end -function c_Own(cS2, NumStages) - c_ts = zeros(NumStages); +function c_own(c_s2, num_stages) + c_ts = zeros(num_stages) # Last timesteps as for SSPRK33 - c_ts[NumStages] = 0.5; - c_ts[NumStages - 1] = 1; + c_ts[num_stages] = 0.5 + c_ts[num_stages - 1] = 1 # Linear increasing timestep for remainder - for i = 2:NumStages-2 - c_ts[i] = cS2 * (i-1)/(NumStages - 3); + for i in 2:num_stages-2 + c_ts[i] = c_s2 * (i-1)/(num_stages - 3) end return c_ts end -function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSemidiscretization) +function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::AbstractSemidiscretization) - #Initialize array of c - c_ts = c_Own(1.0, NumStages) + # Initialize array of c + c_ts = c_own(1.0, num_stages) - #Initialize the array of our solution - aUnknown = zeros(NumStageEvals) + # Initialize the array of our solution + a_unknown = zeros(num_stage_evals) - #Case of e = 3 - if NumStageEvals == 3 - aUnknown = [0, c_ts[2], 0.25] + # Special case of e = 3 + if num_stage_evals == 3 + a_unknown = [0, c_ts[2], 0.25] else - #Calculate MonCoeffs from polynomial Optimizer - ConsOrder = 3 - dtMax = 1.0 - dtEps = 1e-9 + # Calculate coefficients of the stability polynomial in monomial form + cons_order = 3 + dt_max = 1.0 + dt_eps = 1e-9 filter_threshold = 1e-12 - #Get EigVals + # Compute spectrum J = jacobian_ad_forward(semi) - EigVals = eigvals(J) - NumEigVals, EigVals = filter_Eigvals(EigVals, filter_threshold) + eig_vals = eigvals(J) + num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - MonCoeffs, dt, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) - MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + mon_coeffs, dt, _ = bisection(cons_order, num_eig_vals, num_stages, dt_max, dt_eps, eig_vals) + mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) - #Define the objective_function + # Define the objective_function function objective_function(x) - return F_own(x, NumStages, NumStageEvals, MonCoeffs) + return f_own(x, num_stages, num_stage_evals, mon_coeffs) end - #call nlsolver to solve until the answer is not NaN or negative values + # Call nlsolver to solve repeatedly until the result is not NaN or negative values is_sol_valid = false - while is_sol_valid == false - #Initialize initial guess - x0 = 0.1 .* rand(NumStageEvals) + while !is_sol_valid + # Initialize initial guess + x0 = 0.1 .* rand(num_stage_evals) x0[1] = 0.0 x0[2] = c_ts[2] - sol = nlsolve(objective_function, method = :trust_region, x0, - ftol = 1e-15, iterations = 10^4, xtol = 1e-13) + sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 1e-15, iterations = 10^4, xtol = 1e-13) - aUnknown = sol.zero + a_unknown = sol.zero - #First check the values of a[i, i-1] which are acquired by nonlinear solver and c[i] - a[i, i-1]. - is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) + # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 + is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- a_unknown[3:end]) end end - println("aUnknown") - println(aUnknown[3:end]) #To debug + println("a_unknown") + println(a_unknown[3:end]) # To debug - AMatrix = zeros(NumStages - 2, 2) - AMatrix[:, 1] = c_ts[3:end] - AMatrix[:, 1] -= aUnknown[3:end] - AMatrix[:, 2] = aUnknown[3:end] + a_matrix = zeros(num_stages - 2, 2) + a_matrix[:, 1] = c_ts[3:end] + a_matrix[:, 1] -= a_unknown[3:end] + a_matrix[:, 2] = a_unknown[3:end] - return AMatrix, c_ts + return a_matrix, c_ts end -function ComputePERK3_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, cS2::Float64) +function compute_PERK3_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, c_s2::Float64) - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(NumStages) - for k in 2:NumStages-2 - c[k] = cS2 * (k - 1)/(NumStages - 3) # Equidistant timestep distribution (similar to PERK2) - end + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages-2 + c[k] = c_s2 * (k - 1)/(num_stages - 3) # Equidistant timestep distribution (similar to PERK2) + end - # Original third order proposed PERK - #= - c[NumStages - 1] = 1.0/3.0 - c[NumStages] = 1.0 - =# - # Own third order PERK based on SSPRK33 - c[NumStages - 1] = 1.0 - c[NumStages] = 0.5 - - println("Timestep-split: "); display(c); println("\n") + # Original third order proposed PERK + #= + c[num_stages - 1] = 1.0/3.0 + c[num_stages] = 1.0 + =# + # Own third order PERK based on SSPRK33 + c[num_stages - 1] = 1.0 + c[num_stages] = 0.5 + + println("Timestep-split: "); display(c); println("\n") - # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) - CoeffsMax = NumStages - 2 + # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) + coeffs_max = num_stages - 2 - AMatrix = zeros(CoeffsMax, 2) - AMatrix[:, 1] = c[3:end] + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] - PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * "_" * string(NumStages) * ".txt" - NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) - @assert NumMonCoeffs == CoeffsMax + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * string(num_stages) * ".txt" + num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) + @assert num_mon_coeffs == coeffs_max - AMatrix[:, 1] -= A - AMatrix[:, 2] = A + a_matrix[:, 1] -= A + a_matrix[:, 2] = A - println("A matrix: "); display(AMatrix); println() + println("A matrix: "); display(AMatrix); println() - return AMatrix, c + return a_matrix, c end """ @@ -163,31 +162,31 @@ CarpenterKennedy2N{54, 43} methods. """ mutable struct PERK3 <: PERKSingle - const NumStages::Int + const num_stages::Int - AMatrix::Matrix{Float64} - c::Vector{Float64} + a_matrix::Matrix{Float64} + c::Vector{Float64} - # Constructor for previously computed A Coeffs - function PERK3(NumStages_::Int, BasePathMonCoeffs_::AbstractString, cS2_::Float64=1.0) + # Constructor for previously computed A Coeffs + function PERK3(num_stages_::Int, base_path_mon_coeffs_::AbstractString, c_s2_::Float64=1.0) + newPERK3 = new(num_stages_) - newPERK3 = new(NumStages_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, base_path_mon_coeffs_, c_s2_) + + return newPERK3 + end - newPERK3.AMatrix, newPERK3.c = - ComputePERK3_ButcherTableau(NumStages_, BasePathMonCoeffs_, cS2_) - return newPERK3 - end + # Constructor that computes Butcher matrix A coefficients from a semidiscretization + function PERK3(num_stages_::Int, num_stage_evals_ ::Int, semi_::AbstractSemidiscretization) - #Constructor that compute A coefficients from semidiscretization - function PERK3(NumStages_::Int, NumStageEvals_ ::Int, semi_::AbstractSemidiscretization) + newPERK3 = new(num_stages_) - newPERK3 = new(NumStages_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, num_stage_evals_, semi_) - newPERK3.AMatrix, newPERK3.c = - ComputePERK3_ButcherTableau(NumStages_, NumStageEvals_, semi_) + display(newPERK3.a_matrix) - return newPERK3 - end + return newPERK3 + end end # struct PERK3 diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 1ad29cd2f0c..acfb9128f71 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -3,139 +3,143 @@ using Convex using ECOS const MOI = Convex.MOI -function filter_Eigvals(EigVals::Array{Complex{Float64}}, threshold:: Float64) +function filter_eigvals(eig_vals::Array{Complex{Float64}}, threshold:: Float64) + filtered_eigvals_counter = 0 - #FilteredEigVals is not an Array but the code seems to work fine regardless - FilteredEigVals = Complex{Float64}[] - for EigVal in EigVals - if abs(EigVal) < threshold + filtered_eig_vals = Complex{Float64}[] + + for eig_val in eig_vals + if abs(eig_val) < threshold filtered_eigvals_counter += 1 else - push!(FilteredEigVals, EigVal) + push!(filtered_eig_vals, eig_val) end end println("Total number of $filtered_eigvals_counter eigenvalues has not been recorded because they were smaller than $threshold") - return length(FilteredEigVals), FilteredEigVals + return length(filtered_eig_vals), filtered_eig_vals end -function ReadInEigVals(Path_toEvalFile::AbstractString) +function read_in_eig_vals(path_to_eval_file::AbstractString) + + # Declare and set to some value + num_eig_vals = -1 + + open(path_to_eval_file, "r") do eval_file + num_eig_vals = countlines(eval_file) + end - NumEigVals = -1 # Declare and set to some value - open(Path_toEvalFile, "r") do EvalFile - NumEigVals = countlines(EvalFile) - end + eig_vals = Array{Complex{Float64}}(undef, num_eig_vals) - EigVals = Array{Complex{Float64}}(undef, NumEigVals) + line_index = 0 - LineIndex = 0 - open(Path_toEvalFile, "r") do EvalFile - # read till end of file - while !eof(EvalFile) + open(path_to_eval_file, "r") do eval_file + # Read till end of file + while !eof(eval_file) - # read a new / next line for every iteration - LineContent = readline(EvalFile) + # Read a new / next line for every iteration + line_content = readline(eval_file) - EigVals[LineIndex + 1] = parse(Complex{Float64}, LineContent) + eig_vals[line_index + 1] = parse(Complex{Float64}, line_content) - LineIndex += 1 + line_index += 1 + end end - end - return NumEigVals, EigVals + return num_eig_vals, eig_vals end -function Polynoms(ConsOrder::Int, NumStageEvals::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) - - for i in 1:NumEigVals - pnoms[i] = 1.0 - end +function polynoms(cons_order::Int, num_stage_evals::Int, num_eig_vals::Int, normalized_powered_eigvals_scaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) + for i in 1:num_eig_vals + pnoms[i] = 1.0 + end - for k in 1:ConsOrder - for i in 1:NumEigVals - pnoms[i] += normalized_powered_EigvalsScaled[i,k] + for k in 1:cons_order + for i in 1:num_eig_vals + pnoms[i] += normalized_powered_eigvals_scaled[i,k] + end end - end - for k in ConsOrder + 1:NumStageEvals - pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] - end + for k in cons_order + 1:num_stage_evals + pnoms += gamma[k - cons_order] * normalized_powered_eigvals_scaled[:,k] + end - return maximum(abs(pnoms)) + return maximum(abs(pnoms)) end -function Bisection(ConsOrder::Int, NumEigVals::Int, NumStageEvals::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) +function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_max::Float64, dt_eps::Float64, eig_vals::Array{Complex{Float64}}) - dtMin = 0.0 + dt_min = 0.0 - dt = -1.0 + dt = -1.0 - AbsP = -1.0 + abs_p = -1.0 - pnoms = ones(Complex{Float64}, NumEigVals, 1) + pnoms = ones(Complex{Float64}, num_eig_vals, 1) - gamma = Variable(NumStageEvals - ConsOrder) # Init datastructure for results + # Init datastructure for results + gamma = Variable(num_stage_evals - cons_order) - normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStageEvals) + normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - for j in 1:NumStageEvals - fac_j = factorial(j) - for i in 1:NumEigVals - normalized_powered_Eigvals[i, j] = EigVals[i]^j / fac_j + 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 - end - normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStageEvals) + normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - while dtMax - dtMin > dtEps - dt = 0.5 * (dtMax + dtMin) + while dt_max - dt_min > dt_eps + dt = 0.5 * (dt_max + dt_min) - for k in 1:NumStageEvals - dt_k = dt^k - for i in 1:NumEigVals - normalized_powered_EigvalsScaled[i,k] = dt_k * normalized_powered_Eigvals[i,k] - end - end + 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 gamm0 in (potentially) next iteration - problem = minimize(Polynoms(ConsOrder, NumStageEvals, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) - - Convex.solve!( - problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - #"eps" => 1e9, # 1e-13 - "feastol" => 1e-9, # 1e-9 - "abstol" => 1e-9, # 1e-9 - "reltol" => 1e-9, # 1e-9 - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent_solver = false - ) - - AbsP = problem.optval - - println("Current MaxAbsP: ", AbsP, "\nCurrent dt: ", dt, "\n") - - if AbsP < 1.0 - dtMin = dt - else - dtMax = dt + # Use last optimal values for gamm0 in (potentially) next iteration + problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma)) + + Convex.solve!( + problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + #"eps" => 1e9, # 1e-13 + "feastol" => 1e-9, # 1e-9 + "abstol" => 1e-9, # 1e-9 + "reltol" => 1e-9, # 1e-9 + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent_solver = false + ) + + abs_p = problem.optval + + println("Current MaxAbsP: ", abs_p, "\nCurrent dt: ", dt, "\n") + + if abs_p < 1.0 + dt_min = dt + else + dt_max = dt + end end - end - return evaluate(gamma), AbsP, dt + return evaluate(gamma), abs_p, dt end -function undo_normalization(ConsOrder::Int, NumStageEvals::Int, gammaOpt) - for k in ConsOrder + 1:NumStageEvals - gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) - end - return gammaOpt +function undo_normalization!(cons_order::Int, num_stage_evals::Int, gamma_opt) + for k in cons_order + 1:num_stage_evals + gamma_opt[k - cons_order] = gamma_opt[k - cons_order] / factorial(k) + end + return gamma_opt end \ No newline at end of file From 483015bdafdca28ff5b9e445a496f5eb5239e488 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 15 Feb 2024 20:52:25 +0100 Subject: [PATCH 10/64] spelling fix --- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index acfb9128f71..534ee883909 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -103,7 +103,7 @@ function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_ end end - # Use last optimal values for gamm0 in (potentially) next iteration + # Use last optimal values for gamma in (potentially) next iteration problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma)) Convex.solve!( From 03f0af14105dae821f9a1d18849c3b9af61f7321 Mon Sep 17 00:00:00 2001 From: Warisa Date: Sun, 3 Mar 2024 19:15:46 +0100 Subject: [PATCH 11/64] Change names and spacing according to style guide --- .../methods_PERK2.jl | 578 +++++++++--------- .../methods_PERK3.jl | 387 ++++++------ 2 files changed, 497 insertions(+), 468 deletions(-) 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 f8e061c0872..4e113d8720d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -8,106 +8,112 @@ abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end -function ComputeACoeffs(NumStageEvals::Int, - SE_Factors::Vector{Float64}, MonCoeffs::Vector{Float64}) - ACoeffs = MonCoeffs - - for stage in 1:NumStageEvals - 2 - ACoeffs[stage] /= SE_Factors[stage] - for prev_stage in 1:stage-1 - ACoeffs[stage] /= ACoeffs[prev_stage] +function compute_a_coeffs(num_stage_evals::Int, + se_factors::Vector{Float64}, mon_coeffs::Vector{Float64}) + a_coeffs = mon_coeffs + + for stage in 1:(num_stage_evals - 2) + a_coeffs[stage] /= se_factors[stage] + for prev_stage in 1:(stage - 1) + a_coeffs[stage] /= a_coeffs[prev_stage] + end end - end - return reverse(ACoeffs) + return reverse(a_coeffs) end -function ComputePERK2_ButcherTableau(NumStages::Int, semi::AbstractSemidiscretization, bS::Float64, cEnd::Float64) - - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(NumStages) - for k in 2:NumStages - c[k] = cEnd * (k - 1)/(NumStages - 1) - end - println("Timestep-split: "); display(c); println("\n") - SE_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) - CoeffsMax = NumStages - 2 - - AMatrix = zeros(CoeffsMax, 2) - AMatrix[:, 1] = c[3:end] - - - ConsOrder = 2 - filter_thres = 1e-12 - dtMax = 1.0 - dtEps = 1e-9 - - J = jacobian_ad_forward(semi) - EigVals = eigvals(J) - - NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) - - MonCoeffs, PWorstCase, dtOpt = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) - MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) - - NumMonCoeffs = length(MonCoeffs) - @assert NumMonCoeffs == CoeffsMax - A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) - - - #= - # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) - PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * ".txt" - NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) - @assert NumMonCoeffs == CoeffsMax - =# - - AMatrix[:, 1] -= A - AMatrix[:, 2] = A - - println("A matrix: "); display(AMatrix); println() - - return AMatrix, c +function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscretization, + b_s::Float64, c_end::Float64) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages + c[k] = c_end * (k - 1) / (num_stages - 1) + end + println("Timestep-split: ") + display(c) + println("\n") + se_factors = b_s * 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] + + cons_order = 2 + filter_thres = 1e-12 + dtmax = 1.0 + dt_eps = 1e-9 + + eig_vals = eigvals(jacobian_ad_forward(semi)) + + num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) + + mon_coeffs, p_worst_case, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, + eig_vals) + mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) + + num_mon_coeffs = length(mon_coeffs) + @assert num_mon_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) + + #= + # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * ".txt" + num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) + @assert num_mon_coeffs == coeffs_max + =# + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + println("A matrix: ") + display(a_matrix) + println() + + return a_matrix, c end -function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, bS::Float64, cEnd::Float64) - - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(NumStages) - for k in 2:NumStages - c[k] = cEnd * (k - 1)/(NumStages - 1) - end - println("Timestep-split: "); display(c); println("\n") - SE_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) - CoeffsMax = NumStages - 2 - - AMatrix = zeros(CoeffsMax, 2) - AMatrix[:, 1] = c[3:end] - - - PathMonCoeffs = BasePathMonCoeffs * "gamma_" * string(NumStages) * ".txt" - NumMonCoeffs, MonCoeffs = read_file(PathMonCoeffs, Float64) - @assert NumMonCoeffs == CoeffsMax - A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) - - - #= - # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) - PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * ".txt" - NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) - @assert NumMonCoeffs == CoeffsMax - =# - - AMatrix[:, 1] -= A - AMatrix[:, 2] = A - - println("A matrix: "); display(AMatrix); println() - - return AMatrix, c +function compute_PERK2_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, + b_s::Float64, c_end::Float64) + + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages + c[k] = c_end * (k - 1) / (num_stages - 1) + end + println("Timestep-split: ") + display(c) + println("\n") + se_factors = b_s * 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_mon_coeffs = base_path_mon_coeffs * "gamma_" * string(num_stages) * ".txt" + num_mon_coeffs, mon_coeffs = read_file(path_mon_coeffs, Float64) + @assert num_mon_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) + + #= + # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * ".txt" + num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) + @assert num_mon_coeffs == coeffs_max + =# + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + println("A matrix: ") + display(a_matrix) + println() + + return a_matrix, c end """ @@ -122,245 +128,249 @@ CarpenterKennedy2N{54, 43} methods. """ mutable struct PERK2 <: PERKSingle - const NumStages::Int - - AMatrix::Matrix{Float64} - c::Vector{Float64} - bS::Float64 - b1::Float64 - cEnd::Float64 - - #Constructor that read the coefficients from the file - function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5) - - newPERK2 = new(NumStages_) - - newPERK2.AMatrix, newPERK2.c = - ComputePERK2_ButcherTableau(NumStages_, BasePathMonCoeffs_, bS_, cEnd_) - - newPERK2.b1 = one(bS_) - bS_ - newPERK2.bS = bS_ - newPERK2.cEnd = cEnd_ - return newPERK2 - end - - #Constructor that calculate the coefficients with polynomial optimizer - function PERK2(NumStages_::Int, semi_::AbstractSemidiscretization, bS_::Float64=1.0, cEnd_::Float64=0.5) + const num_stages::Int + + a_matrix::Matrix{Float64} + c::Vector{Float64} + b_s::Float64 + b1::Float64 + c_end::Float64 + + #Constructor that read the coefficients from the file + function PERK2(num_stages_::Int, base_path_mon_coeffs_::AbstractString, b_s_::Float64 = 1.0, + c_end_::Float64 = 0.5) + newPERK2 = new(num_stages_) + + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, + base_path_mon_coeffs_, b_s_, + c_end_) + + newPERK2.b1 = one(b_s_) - b_s_ + newPERK2.b_s = b_s_ + newPERK2.c_end = c_end_ + return newPERK2 + end - newPERK2 = new(NumStages_) + #Constructor that calculate the coefficients with polynomial optimizer + function PERK2(num_stages_::Int, semi_::AbstractSemidiscretization, b_s_::Float64 = 1.0, + c_end_::Float64 = 0.5) + newPERK2 = new(num_stages_) - newPERK2.AMatrix, newPERK2.c = - ComputePERK2_ButcherTableau(NumStages_, semi_, bS_, cEnd_) + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, semi_, b_s_, + c_end_) - newPERK2.b1 = one(bS_) - bS_ - newPERK2.bS = bS_ - newPERK2.cEnd = cEnd_ - return newPERK2 - end + newPERK2.b1 = one(b_s_) - b_s_ + newPERK2.b_s = b_s_ + newPERK2.c_end = c_end_ + return newPERK2 + end end # struct PERK2 - # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 -mutable struct PERK_IntegratorOptions{Callback} - callback::Callback # callbacks; used in Trixi - adaptive::Bool # whether the algorithm is adaptive; ignored - dtmax::Float64 # ignored - maxiters::Int # maximal numer of time steps - tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored +mutable struct PERKIntegratorOptions{Callback} + callback::Callback # callbacks; used in Trixi + adaptive::Bool # whether the algorithm is adaptive; ignored + dtmax::Float64 # ignored + maxiters::Int # maximal numer of time steps + tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored end -function PERK_IntegratorOptions(callback, tspan; maxiters=typemax(Int), kwargs...) - PERK_IntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, [last(tspan)]) +function PERKIntegratorOptions(callback, tspan; maxiters = typemax(Int), kwargs...) + PERKIntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, [last(tspan)]) end -abstract type PERK_Integrator end -abstract type PERKSingle_Integrator <: PERK_Integrator end +abstract type PERKIntegrator end +abstract type PERKSingleIntegrator <: PERKIntegrator 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 PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK_IntegratorOptions} <: PERKSingle_Integrator - 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::PERK_IntegratorOptions - finalstep::Bool # added for convenience - # PERK2 stages: - k1::uType - k_higher::uType - t_stage::RealT +mutable struct PERK2Integrator{RealT <: Real, uType, Params, Sol, F, Alg, + PERKIntegratorOptions} <: PERKSingleIntegrator + 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::PERKIntegratorOptions + finalstep::Bool # added for convenience + # PERK2 stages: + k1::uType + k_higher::uType + t_stage::RealT end # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) -function Base.getproperty(integrator::PERK_Integrator, field::Symbol) - if field === :stats - return (naccept = getfield(integrator, :iter),) - end - # general fallback - return getfield(integrator, field) +function Base.getproperty(integrator::PERKIntegrator, field::Symbol) + if field === :stats + return (naccept = getfield(integrator, :iter),) + end + # general fallback + return getfield(integrator, field) end # Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 function solve(ode::ODEProblem, alg::PERK2; - dt, callback=nothing, kwargs...) - - u0 = copy(ode.u0) - du = zero(u0) #previously: similar(u0) - u_tmp = zero(u0) - - # PERK2 stages - k1 = zero(u0) - k_higher = zero(u0) - - t0 = first(ode.tspan) - iter = 0 - - integrator = PERK2_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, - (prob=ode,), ode.f, alg, - PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher, t0) - - # 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) + dt, callback = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(u0) #previously: similar(u0) + u_tmp = zero(u0) + + # PERK2 stages + k1 = zero(u0) + k_higher = zero(u0) + + t0 = first(ode.tspan) + iter = 0 + + integrator = PERK2Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, + (prob = ode,), ode.f, alg, + PERKIntegratorOptions(callback, ode.tspan; kwargs...), + false, + k1, k_higher, t0) + + # 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 - elseif !isnothing(callback) - error("unsupported") - end - solve!(integrator) + solve!(integrator) end -function solve!(integrator::PERK2_Integrator) - @unpack prob = integrator.sol - @unpack alg = integrator - t_end = last(prob.tspan) - callbacks = integrator.opts.callback +function solve!(integrator::PERK2Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback - integrator.finalstep = false + integrator.finalstep = false - @trixi_timeit timer() "main loop" while !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 - - # k2 - integrator.t_stage = integrator.t + alg.c[2] * integrator.dt - - # Construct current state - @threaded for i in eachindex(integrator.du) - integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] - end - - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) - - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt - end - - # Higher stages - for stage = 3:alg.NumStages - integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt - - # Construct current state - @threaded for i in eachindex(integrator.du) - integrator.u_tmp[i] = integrator.u[i] + alg.AMatrix[stage - 2, 1] * integrator.k1[i] + - alg.AMatrix[stage - 2, 2] * integrator.k_higher[i] + @trixi_timeit timer() "main loop" while !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") end - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + # 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 - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt + @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 + + # k2 + integrator.t_stage = integrator.t + alg.c[2] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @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) + integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt + + # 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_stage) + + @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] += integrator.k_higher[i] + integrator.u[i] += alg.b1 * integrator.k1[i] + + alg.b_s * integrator.k_higher[i] + end + end # PERK2 step + + 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 - end - - @threaded for i in eachindex(integrator.u) - #integrator.u[i] += integrator.k_higher[i] - integrator.u[i] += alg.b1 * integrator.k1[i] + alg.bS * integrator.k_higher[i] - end - end # PERK2 step - - 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) + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) end - end - end + end # "main loop" timer - # respect maximum number of iterations - if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep - @warn "Interrupted. Larger maxiters is needed." - terminate!(integrator) - end - end # "main loop" timer - - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) end # get a cache where the RHS can be stored -get_du(integrator::PERK_Integrator) = integrator.du -get_tmp_cache(integrator::PERK_Integrator) = (integrator.u_tmp,) +get_du(integrator::PERKIntegrator) = integrator.du +get_tmp_cache(integrator::PERKIntegrator) = (integrator.u_tmp,) # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u -u_modified!(integrator::PERK_Integrator, ::Bool) = false +u_modified!(integrator::PERKIntegrator, ::Bool) = false # used by adaptive timestepping algorithms in DiffEq -function set_proposed_dt!(integrator::PERK_Integrator, dt) - integrator.dt = dt +function set_proposed_dt!(integrator::PERKIntegrator, dt) + integrator.dt = dt end -function get_proposed_dt(integrator::PERK_Integrator) - return integrator.dt +function get_proposed_dt(integrator::PERKIntegrator) + return integrator.dt end # stop the time integration -function terminate!(integrator::PERK_Integrator) - integrator.finalstep = true - empty!(integrator.opts.tstops) +function terminate!(integrator::PERKIntegrator) + integrator.finalstep = true + empty!(integrator.opts.tstops) end # used for AMR (Adaptive Mesh Refinement) -function Base.resize!(integrator::PERK2_Integrator, new_size) - resize!(integrator.u, new_size) - resize!(integrator.du, new_size) - resize!(integrator.u_tmp, new_size) +function Base.resize!(integrator::PERK2Integrator, 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) + resize!(integrator.k1, new_size) + resize!(integrator.k_higher, new_size) end -end # @muladd \ No newline at end of file +end \ No newline at end of file diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 5901d3b5e65..5fd1dd19dc5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -18,8 +18,8 @@ function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) term1 *= a_unknown[num_stage_evals - 1 - j] term2 *= a_unknown[num_stage_evals - j] end - term1 *= c_ts[num_stages - 2 - i] * 1/6 - term2 *= c_ts[num_stages - 1 - i] * 4/6 + term1 *= c_ts[num_stages - 2 - i] * 1 / 6 + term2 *= c_ts[num_stages - 1 - i] * 4 / 6 c_eq[i] = mon_coeffs[i] - (term1 + term2) end @@ -30,11 +30,12 @@ function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) for j in 1:i term2 *= a_unknown[num_stage_evals - j] end - term2 *= c_ts[num_stages - 1 - i] * 4/6 + term2 *= c_ts[num_stages - 1 - i] * 4 / 6 c_eq[i] = mon_coeffs[i] - term2 - c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - a_unknown[num_stage_evals - 1] + c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - + a_unknown[num_stage_evals - 1] return c_eq end @@ -43,17 +44,18 @@ function c_own(c_s2, num_stages) c_ts = zeros(num_stages) # Last timesteps as for SSPRK33 - c_ts[num_stages] = 0.5 + c_ts[num_stages] = 0.5 c_ts[num_stages - 1] = 1 # Linear increasing timestep for remainder - for i in 2:num_stages-2 - c_ts[i] = c_s2 * (i-1)/(num_stages - 3) + for i in 2:(num_stages - 2) + c_ts[i] = c_s2 * (i - 1) / (num_stages - 3) end return c_ts end -function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::AbstractSemidiscretization) +function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, + semi::AbstractSemidiscretization) # Initialize array of c c_ts = c_own(1.0, num_stages) @@ -68,7 +70,7 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::Abstra else # Calculate coefficients of the stability polynomial in monomial form cons_order = 3 - dt_max = 1.0 + dtmax = 1.0 dt_eps = 1e-9 filter_threshold = 1e-12 @@ -77,7 +79,8 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::Abstra eig_vals = eigvals(J) num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - mon_coeffs, dt, _ = bisection(cons_order, num_eig_vals, num_stages, dt_max, dt_eps, eig_vals) + mon_coeffs, dt, _ = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, + eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) # Define the objective_function @@ -93,12 +96,14 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::Abstra x0[1] = 0.0 x0[2] = c_ts[2] - sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 1e-15, iterations = 10^4, xtol = 1e-13) - + sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 1e-15, + iterations = 10^4, xtol = 1e-13) + a_unknown = sol.zero # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 - is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- a_unknown[3:end]) + is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && + all(x -> !isnan(x) && x >= 0, c_ts[3:end] .- a_unknown[3:end]) end end @@ -108,19 +113,20 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::Abstra a_matrix = zeros(num_stages - 2, 2) a_matrix[:, 1] = c_ts[3:end] a_matrix[:, 1] -= a_unknown[3:end] - a_matrix[:, 2] = a_unknown[3:end] - + a_matrix[:, 2] = a_unknown[3:end] + return a_matrix, c_ts end -function compute_PERK3_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, c_s2::Float64) +function compute_PERK3_butcher_tableau(num_stages::Int, + base_path_mon_coeffs::AbstractString, c_s2::Float64) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) - for k in 2:num_stages-2 - c[k] = c_s2 * (k - 1)/(num_stages - 3) # Equidistant timestep distribution (similar to PERK2) + for k in 2:(num_stages - 2) + c[k] = c_s2 * (k - 1) / (num_stages - 3) # Equidistant timestep distribution (similar to PERK2) end - + # Original third order proposed PERK #= c[num_stages - 1] = 1.0/3.0 @@ -128,24 +134,29 @@ function compute_PERK3_butcher_tableau(num_stages::Int, base_path_mon_coeffs::Ab =# # Own third order PERK based on SSPRK33 c[num_stages - 1] = 1.0 - c[num_stages] = 0.5 + c[num_stages] = 0.5 + + println("Timestep-split: ") + display(c) + println("\n") - println("Timestep-split: "); display(c); println("\n") - # - 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_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * string(num_stages) * ".txt" + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * + string(num_stages) * ".txt" num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) @assert num_mon_coeffs == coeffs_max a_matrix[:, 1] -= A - a_matrix[:, 2] = A - - println("A matrix: "); display(AMatrix); println() + a_matrix[:, 2] = A + + println("A matrix: ") + display(a_matrix) + println() return a_matrix, c end @@ -168,202 +179,210 @@ mutable struct PERK3 <: PERKSingle c::Vector{Float64} # Constructor for previously computed A Coeffs - function PERK3(num_stages_::Int, base_path_mon_coeffs_::AbstractString, c_s2_::Float64=1.0) + function PERK3(num_stages_::Int, base_path_mon_coeffs_::AbstractString, + c_s2_::Float64 = 1.0) newPERK3 = new(num_stages_) - newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, base_path_mon_coeffs_, c_s2_) - + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, + base_path_mon_coeffs_, + c_s2_) + return newPERK3 end # Constructor that computes Butcher matrix A coefficients from a semidiscretization - function PERK3(num_stages_::Int, num_stage_evals_ ::Int, semi_::AbstractSemidiscretization) - + function PERK3(num_stages_::Int, num_stage_evals_::Int, + semi_::AbstractSemidiscretization) newPERK3 = new(num_stages_) - newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, num_stage_evals_, semi_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, + num_stage_evals_, + semi_) display(newPERK3.a_matrix) return newPERK3 end - end # struct PERK3 - # 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 PERK3_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK_IntegratorOptions} - 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::PERK_IntegratorOptions - finalstep::Bool # added for convenience - # PERK stages: - k1::uType - k_higher::uType - k_S1::uType # Required for third order - t_stage::RealT +mutable struct PERK3Integrator{RealT <: Real, uType, Params, Sol, F, Alg, + PERKIntegratorOptions} + u::uType + du::uType + u_tmp::uType + t::RealT + dt::RealT # current time step + dt_cache::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::PERKIntegratorOptions + final_step::Bool # added for convenience + # PERK stages: + k1::uType + k_higher::uType + k_s1::uType # Required for third order + t_stage::RealT end # Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 function solve(ode::ODEProblem, alg::PERK3; - dt, callback=nothing, kwargs...) - - u0 = copy(ode.u0) - du = zero(u0) #previously: similar(u0) - u_tmp = zero(u0) - - # PERK stages - k1 = zero(u0) - k_higher = zero(u0) - k_S1 = zero(u0) - - t0 = first(ode.tspan) - iter = 0 - - integrator = PERK3_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, - (prob=ode,), ode.f, alg, - PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher, k_S1, t0) - - # 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) + dt, callback = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(u0) #previously: similar(u0) + u_tmp = zero(u0) + + # PERK stages + k1 = zero(u0) + k_higher = zero(u0) + k_s1 = zero(u0) + + t0 = first(ode.tspan) + iter = 0 + + integrator = PERK3Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, + (prob = ode,), ode.f, alg, + PERKIntegratorOptions(callback, ode.tspan; kwargs...), + false, + k1, k_higher, k_s1, t0) + + # 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 - elseif !isnothing(callback) - error("unsupported") - end - solve!(integrator) + solve!(integrator) end -function solve!(integrator::PERK3_Integrator) - @unpack prob = integrator.sol - @unpack alg = integrator - t_end = last(prob.tspan) - callbacks = integrator.opts.callback - - integrator.finalstep = false +function solve!(integrator::PERK3Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback - #@trixi_timeit timer() "main loop" while !integrator.finalstep - while !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 + integrator.final_step = false - @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 - - # k2 - integrator.t_stage = integrator.t + alg.c[2] * integrator.dt - - # Construct current state - @threaded for i in eachindex(integrator.du) - integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] - end - - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) - - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt - end - - if alg.NumStages == 3 - @threaded for i in eachindex(integrator.du) - integrator.k_S1[i] = integrator.k_higher[i] + #@trixi_timeit timer() "main loop" while !integrator.final_step + while !integrator.final_step + if isnan(integrator.dt) + error("time step size `dt` is NaN") end - end - - # Higher stages - for stage = 3:alg.NumStages - integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt - - # Construct current state - @threaded for i in eachindex(integrator.du) - integrator.u_tmp[i] = integrator.u[i] + alg.AMatrix[stage - 2, 1] * integrator.k1[i] + - alg.AMatrix[stage - 2, 2] * integrator.k_higher[i] - end - - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt + # 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 - # IDEA: Stop for loop at NumStages -1 to avoid if (maybe more performant?) - if stage == alg.NumStages - 1 - @threaded for i in eachindex(integrator.du) - integrator.k_S1[i] = integrator.k_higher[i] - 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 + + # k2 + integrator.t_stage = integrator.t + alg.c[2] * integrator.dt + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + if alg.num_stages == 3 + @threaded for i in eachindex(integrator.du) + integrator.k_s1[i] = integrator.k_higher[i] + end + end + + # Higher stages + for stage in 3:(alg.num_stages) + integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt + + # 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_stage) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + # IDEA: Stop for loop at num_stages -1 to avoid if (maybe more performant?) + if stage == alg.num_stages - 1 + @threaded for i in eachindex(integrator.du) + integrator.k_s1[i] = integrator.k_higher[i] + end + end + end + + @threaded for i in eachindex(integrator.u) + # Original proposed PERK3 + #integrator.u[i] += 0.75 * integrator.k_s1[i] + 0.25 * integrator.k_higher[i] + # Own PERK based on SSPRK33 + integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + + 4.0 * integrator.k_higher[i]) / 6.0 + end + end # PERK 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 - end - - @threaded for i in eachindex(integrator.u) - # Original proposed PERK3 - #integrator.u[i] += 0.75 * integrator.k_S1[i] + 0.25 * integrator.k_higher[i] - # Own PERK based on SSPRK33 - integrator.u[i] += (integrator.k1[i] + integrator.k_S1[i] + 4.0 * integrator.k_higher[i])/6.0 - end - end # PERK 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) + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.final_step + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) end - end - end + end # "main loop" timer - # respect maximum number of iterations - if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep - @warn "Interrupted. Larger maxiters is needed." - terminate!(integrator) - end - end # "main loop" timer - - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) end # used for AMR (Adaptive Mesh Refinement) -function Base.resize!(integrator::PERK3_Integrator, 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) - resize!(integrator.k_S1, new_size) +function Base.resize!(integrator::PERK3Integrator, 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) + resize!(integrator.k_s1, new_size) end -end # @muladd \ No newline at end of file +#end # @muladd +end \ No newline at end of file From e182a89b76739a405abcb6190b08527e18148cec Mon Sep 17 00:00:00 2001 From: Warisa Date: Sun, 3 Mar 2024 19:25:00 +0100 Subject: [PATCH 12/64] spelling correction --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 4e113d8720d..6188b9e28d6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -171,7 +171,7 @@ mutable struct PERKIntegratorOptions{Callback} callback::Callback # callbacks; used in Trixi adaptive::Bool # whether the algorithm is adaptive; ignored dtmax::Float64 # ignored - maxiters::Int # maximal numer of time steps + 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 From 08b70d6a01f75c22a586299162f913da6f89fe43 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Mon, 4 Mar 2024 21:23:56 +0100 Subject: [PATCH 13/64] Apply suggestions from code review Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) 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 6188b9e28d6..3f701d3984d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -9,7 +9,7 @@ abstract type PERK end abstract type PERKSingle <: PERK end function compute_a_coeffs(num_stage_evals::Int, - se_factors::Vector{Float64}, mon_coeffs::Vector{Float64}) + se_factors::Vector{Float64}, mon_coeffs::Vector{Float64}) a_coeffs = mon_coeffs for stage in 1:(num_stage_evals - 2) @@ -23,7 +23,7 @@ function compute_a_coeffs(num_stage_evals::Int, end function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscretization, - b_s::Float64, c_end::Float64) + b_s::Float64, c_end::Float64) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -51,7 +51,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) mon_coeffs, p_worst_case, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, - eig_vals) + eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) num_mon_coeffs = length(mon_coeffs) @@ -76,7 +76,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre end function compute_PERK2_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, - b_s::Float64, c_end::Float64) + b_s::Float64, c_end::Float64) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -142,8 +142,8 @@ mutable struct PERK2 <: PERKSingle newPERK2 = new(num_stages_) newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, - base_path_mon_coeffs_, b_s_, - c_end_) + base_path_mon_coeffs_, b_s_, + c_end_) newPERK2.b1 = one(b_s_) - b_s_ newPERK2.b_s = b_s_ @@ -156,8 +156,7 @@ mutable struct PERK2 <: PERKSingle c_end_::Float64 = 0.5) newPERK2 = new(num_stages_) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, semi_, b_s_, - c_end_) + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, semi_, b_s_, c_end_) newPERK2.b1 = one(b_s_) - b_s_ newPERK2.b_s = b_s_ From bdb2db53f4a96bf581cf959f6fecaaa2a75e11c7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Mar 2024 12:09:29 +0100 Subject: [PATCH 14/64] snake case --- src/auxiliary/auxiliary.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 1ee36238276..82ed7284e4c 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -346,17 +346,17 @@ function register_error_hints() return nothing end -function read_file(FilePath::AbstractString, DataType::Type=Float64) - @assert isfile(FilePath) "Couldn't find file" - Data = zeros(DataType, 0) - open(FilePath, "r") do File - while !eof(File) - LineContent = readline(File) - append!(Data, parse(DataType, LineContent)) +function read_file(file_path::AbstractString, data_type::Type=Float64) + @assert isfile(file_path) "Couldn't find file" + data = zeros(data_type, 0) + open(file_path, "r") do file + while !eof(file) + line_content = readline(file) + append!(data, parse(data_type, line_content)) end end - NumLines = length(Data) + num_lines = length(data) - return NumLines, Data + return num_lines, data end end # @muladd From df4733d483b9fd23c66c4ec4b93fbb5357c0a485 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Mar 2024 15:12:15 +0100 Subject: [PATCH 15/64] revisit perk single p2 p3 --- Project.toml | 4 + .../elixir_burgers_PERK3.jl | 53 ++++++ .../tree_1d_dgsem/elixir_advection_PERK2.jl | 16 +- .../methods_PERK2.jl | 90 ++++------ .../methods_PERK3.jl | 166 +++++++----------- .../polynomial_optimizer.jl | 29 +-- 6 files changed, 173 insertions(+), 185 deletions(-) create mode 100644 examples/structured_1d_dgsem/elixir_burgers_PERK3.jl diff --git a/Project.toml b/Project.toml index 1a0aa0103dc..dbe541ede11 100644 --- a/Project.toml +++ b/Project.toml @@ -6,9 +6,11 @@ version = "0.6.8-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +Convex = "f65535da-76fb-5f13-bab9-19810c17039a" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -19,8 +21,10 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" 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 00000000000..056f83a3e99 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_burgers_PERK3.jl @@ -0,0 +1,53 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the (inviscid) Burgers' equation + +equations = InviscidBurgersEquation1D() + +initial_condition = initial_condition_convergence_test + +# Create DG solver with polynomial degree = 3 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,) + +# Create curved mesh with 16 cells +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) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +ode_algorithm = PERK3(8, 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() \ 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 index 84ded3528b8..8f12bfb328f 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi -using Downloads: download ############################################################################### # semidiscretization of the linear advection equation @@ -38,7 +37,7 @@ summary_callback = SummaryCallback() analysis_callback = AnalysisCallback(semi, interval = 100) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 2.5) +stepsize_callback = StepsizeCallback(cfl = 2.5/4) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, analysis_callback, @@ -47,22 +46,11 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation -MonCoeffsFile = "gamma_6.txt" - -download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", -MonCoeffsFile) - -#ode_algorithm = PERK2(6, "./") - ode_algorithm = PERK2(6, semi) -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks 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() - -using Plots -plot(sol) \ No newline at end of file +summary_callback() \ No newline at end of file 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 3f701d3984d..ebff168145a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -4,12 +4,12 @@ # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin +# Abstract base type for both single/standalone and multi-level P-ERK time integration schemes abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end -function compute_a_coeffs(num_stage_evals::Int, - se_factors::Vector{Float64}, mon_coeffs::Vector{Float64}) +function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) a_coeffs = mon_coeffs for stage in 1:(num_stage_evals - 2) @@ -22,8 +22,7 @@ function compute_a_coeffs(num_stage_evals::Int, return reverse(a_coeffs) end -function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscretization, - b_s::Float64, c_end::Float64) +function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -33,7 +32,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre println("Timestep-split: ") display(c) println("\n") - se_factors = b_s * reverse(c[2:(end - 1)]) + se_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 @@ -50,7 +49,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) - mon_coeffs, p_worst_case, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, + mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) @@ -58,13 +57,6 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre @assert num_mon_coeffs == coeffs_max A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) - #= - # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) - path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * ".txt" - num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) - @assert num_mon_coeffs == coeffs_max - =# - a_matrix[:, 1] -= A a_matrix[:, 2] = A @@ -75,8 +67,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, semi::AbstractSemidiscre return a_matrix, c end -function compute_PERK2_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, - b_s::Float64, c_end::Float64) +function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -86,7 +77,7 @@ function compute_PERK2_butcher_tableau(num_stages::Int, base_path_mon_coeffs::Ab println("Timestep-split: ") display(c) println("\n") - se_factors = b_s * reverse(c[2:(end - 1)]) + se_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 @@ -99,13 +90,6 @@ function compute_PERK2_butcher_tableau(num_stages::Int, base_path_mon_coeffs::Ab @assert num_mon_coeffs == coeffs_max A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) - #= - # TODO: Not sure if I not rather want to read-in values (especially those from Many Stage C++ Optim) - path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * ".txt" - num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) - @assert num_mon_coeffs == coeffs_max - =# - a_matrix[:, 1] -= A a_matrix[:, 2] = A @@ -122,9 +106,6 @@ end The following structures and methods provide a minimal implementation of the paired explicit Runge-Kutta method (https://doi.org/10.1016/j.jcp.2019.05.014) optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). - -This is using the same interface as OrdinaryDiffEq.jl, copied from file "methods_2N.jl" for the -CarpenterKennedy2N{54, 43} methods. """ mutable struct PERK2 <: PERKSingle @@ -132,35 +113,34 @@ mutable struct PERK2 <: PERKSingle a_matrix::Matrix{Float64} c::Vector{Float64} - b_s::Float64 + bS::Float64 b1::Float64 c_end::Float64 #Constructor that read the coefficients from the file - function PERK2(num_stages_::Int, base_path_mon_coeffs_::AbstractString, b_s_::Float64 = 1.0, - c_end_::Float64 = 0.5) - newPERK2 = new(num_stages_) + function PERK2(num_stages, base_path_mon_coeffs::AbstractString, bS = 1.0, c_end = 0.5) + newPERK2 = new(num_stages) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, - base_path_mon_coeffs_, b_s_, - c_end_) + newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, + base_path_mon_coeffs, + bS, c_end) - newPERK2.b1 = one(b_s_) - b_s_ - newPERK2.b_s = b_s_ - newPERK2.c_end = c_end_ + newPERK2.b1 = one(bS) - bS + newPERK2.bS = bS + newPERK2.c_end = c_end return newPERK2 end #Constructor that calculate the coefficients with polynomial optimizer - function PERK2(num_stages_::Int, semi_::AbstractSemidiscretization, b_s_::Float64 = 1.0, - c_end_::Float64 = 0.5) - newPERK2 = new(num_stages_) + function PERK2(num_stages, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5) + newPERK2 = new(num_stages) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages_, semi_, b_s_, c_end_) + newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, semi, + bS, c_end) - newPERK2.b1 = one(b_s_) - b_s_ - newPERK2.b_s = b_s_ - newPERK2.c_end = c_end_ + newPERK2.b1 = one(bS) - bS + newPERK2.bS = bS + newPERK2.c_end = c_end return newPERK2 end end # struct PERK2 @@ -203,7 +183,6 @@ mutable struct PERK2Integrator{RealT <: Real, uType, Params, Sol, F, Alg, # PERK2 stages: k1::uType k_higher::uType - t_stage::RealT end # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) @@ -219,7 +198,7 @@ end function solve(ode::ODEProblem, alg::PERK2; dt, callback = nothing, kwargs...) u0 = copy(ode.u0) - du = zero(u0) #previously: similar(u0) + du = zero(u0) u_tmp = zero(u0) # PERK2 stages @@ -233,7 +212,7 @@ function solve(ode::ODEProblem, alg::PERK2; (prob = ode,), ode.f, alg, PERKIntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher, t0) + k1, k_higher) # initialize callbacks if callback isa CallbackSet @@ -271,21 +250,18 @@ function solve!(integrator::PERK2Integrator) end @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin - # k1: + # 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 - # k2 - integrator.t_stage = integrator.t + alg.c[2] * integrator.dt - # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] end - - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + # 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 @@ -293,8 +269,6 @@ function solve!(integrator::PERK2Integrator) # Higher stages for stage in 3:(alg.num_stages) - integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt - # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + @@ -302,7 +276,7 @@ function solve!(integrator::PERK2Integrator) alg.a_matrix[stage - 2, 2] * integrator.k_higher[i] end - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + 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 @@ -310,9 +284,7 @@ function solve!(integrator::PERK2Integrator) end @threaded for i in eachindex(integrator.u) - #integrator.u[i] += integrator.k_higher[i] - integrator.u[i] += alg.b1 * integrator.k1[i] + - alg.b_s * integrator.k_higher[i] + integrator.u[i] += alg.b1 * integrator.k1[i] + alg.bS * integrator.k_higher[i] end end # PERK2 step @@ -372,4 +344,4 @@ function Base.resize!(integrator::PERK2Integrator, new_size) resize!(integrator.k_higher, new_size) end -end \ No newline at end of file +end # @muladd \ No newline at end of file diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 5fd1dd19dc5..69e119f3a27 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -2,14 +2,30 @@ # 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 NLsolve +using NLsolve: nlsolve + @muladd begin -function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) - c_s2 = 1 # c_s2 = c_ts(S-2) - c_ts = c_own(c_s2, num_stages) # ts = timestep +# Initialize Butcher array abscissae c for PERK3 based on SSPRK33 base method +function c_PERK3_SSP33(num_stages, c_S2) + c = zeros(num_stages) + + # Last timesteps as for SSPRK33 + c[num_stages] = 0.5 + c[num_stages - 1] = 1 + + # Linear increasing timestep for remainder + for i in 2:(num_stages - 2) + c[i] = c_S2 * (i - 1) / (num_stages - 3) + end - c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_s2 is equal to 1 + return c +end + +function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, mon_coeffs, c_S2) + c_ts = c_PERK3_SSP33(num_stages, c_S2) # ts = timestep + + c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_S2 is equal to 1 # Both terms should be present for i in 1:(num_stage_evals - 4) term1 = a_unknown[num_stage_evals - 1] @@ -18,8 +34,8 @@ function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) term1 *= a_unknown[num_stage_evals - 1 - j] term2 *= a_unknown[num_stage_evals - j] end - term1 *= c_ts[num_stages - 2 - i] * 1 / 6 - term2 *= c_ts[num_stages - 1 - i] * 4 / 6 + term1 *= c_ts[num_stages - 2 - i] * 1/6 + term2 *= c_ts[num_stages - 1 - i] * 4/6 c_eq[i] = mon_coeffs[i] - (term1 + term2) end @@ -34,38 +50,22 @@ function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) c_eq[i] = mon_coeffs[i] - term2 - c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - - a_unknown[num_stage_evals - 1] + c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - a_unknown[num_stage_evals - 1] return c_eq end -function c_own(c_s2, num_stages) - c_ts = zeros(num_stages) - - # Last timesteps as for SSPRK33 - c_ts[num_stages] = 0.5 - c_ts[num_stages - 1] = 1 - - # Linear increasing timestep for remainder - for i in 2:(num_stages - 2) - c_ts[i] = c_s2 * (i - 1) / (num_stages - 3) - end - return c_ts -end - -function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, - semi::AbstractSemidiscretization) +function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, c_S2) # Initialize array of c - c_ts = c_own(1.0, num_stages) + c = c_PERK3_SSP33(num_stages, c_S2) # Initialize the array of our solution - a_unknown = zeros(num_stage_evals) + a_unknown = zeros(num_stages) # Special case of e = 3 - if num_stage_evals == 3 - a_unknown = [0, c_ts[2], 0.25] + if num_stages == 3 + a_unknown = [0, c[2], 0.25] else # Calculate coefficients of the stability polynomial in monomial form @@ -79,31 +79,30 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, eig_vals = eigvals(J) num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - mon_coeffs, dt, _ = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, - eig_vals) + mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) # Define the objective_function function objective_function(x) - return f_own(x, num_stages, num_stage_evals, mon_coeffs) + return PERK3_Butcher_tableau_objective(x, num_stages, num_stages, mon_coeffs, c_S2) end # Call nlsolver to solve repeatedly until the result is not NaN or negative values is_sol_valid = false while !is_sol_valid # Initialize initial guess - x0 = 0.1 .* rand(num_stage_evals) + x0 = 0.1 .* rand(num_stages) x0[1] = 0.0 - x0[2] = c_ts[2] + x0[2] = c[2] - sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 1e-15, + sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 4e-16, iterations = 10^4, xtol = 1e-13) a_unknown = sol.zero # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && - all(x -> !isnan(x) && x >= 0, c_ts[3:end] .- a_unknown[3:end]) + all(x -> !isnan(x) && x >= 0, c[3:end] .- a_unknown[3:end]) end end @@ -111,30 +110,17 @@ function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, println(a_unknown[3:end]) # To debug a_matrix = zeros(num_stages - 2, 2) - a_matrix[:, 1] = c_ts[3:end] + a_matrix[:, 1] = c[3:end] a_matrix[:, 1] -= a_unknown[3:end] a_matrix[:, 2] = a_unknown[3:end] - return a_matrix, c_ts + return a_matrix, c end -function compute_PERK3_butcher_tableau(num_stages::Int, - base_path_mon_coeffs::AbstractString, c_s2::Float64) - - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(num_stages) - for k in 2:(num_stages - 2) - c[k] = c_s2 * (k - 1) / (num_stages - 3) # Equidistant timestep distribution (similar to PERK2) - end +function compute_PERK3_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, c_S2) - # Original third order proposed PERK - #= - c[num_stages - 1] = 1.0/3.0 - c[num_stages] = 1.0 - =# - # Own third order PERK based on SSPRK33 - c[num_stages - 1] = 1.0 - c[num_stages] = 0.5 + # Initialize array of c + c = c_PERK3_SSP33(num_stages, c_S2) println("Timestep-split: ") display(c) @@ -146,8 +132,7 @@ function compute_PERK3_butcher_tableau(num_stages::Int, a_matrix = zeros(coeffs_max, 2) a_matrix[:, 1] = c[3:end] - path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * - string(num_stages) * ".txt" + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * string(num_stages) * ".txt" num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) @assert num_mon_coeffs == coeffs_max @@ -167,9 +152,6 @@ end The following structures and methods provide a minimal implementation of the third order paired explicit Runge-Kutta method (https://www.sciencedirect.com/science/article/pii/S0021999122005320) optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). - -This is using the same interface as OrdinaryDiffEq.jl, copied from file "methods_2N.jl" for the -CarpenterKennedy2N{54, 43} methods. """ mutable struct PERK3 <: PERKSingle @@ -179,25 +161,23 @@ mutable struct PERK3 <: PERKSingle c::Vector{Float64} # Constructor for previously computed A Coeffs - function PERK3(num_stages_::Int, base_path_mon_coeffs_::AbstractString, - c_s2_::Float64 = 1.0) - newPERK3 = new(num_stages_) + function PERK3(num_stages, base_path_mon_coeffs::AbstractString, c_S2::Float64 = 1.0) + newPERK3 = new(num_stages) - newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, - base_path_mon_coeffs_, - c_s2_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, + base_path_mon_coeffs, + c_S2) return newPERK3 end # Constructor that computes Butcher matrix A coefficients from a semidiscretization - function PERK3(num_stages_::Int, num_stage_evals_::Int, - semi_::AbstractSemidiscretization) - newPERK3 = new(num_stages_) + function PERK3(num_stages, semi::AbstractSemidiscretization, c_S2::Float64 = 1.0) + newPERK3 = new(num_stages) - newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, - num_stage_evals_, - semi_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, + semi, + c_S2) display(newPERK3.a_matrix) @@ -210,32 +190,31 @@ end # struct PERK3 # https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1 # which are used in Trixi. mutable struct PERK3Integrator{RealT <: Real, uType, Params, Sol, F, Alg, - PERKIntegratorOptions} + PERKIntegratorOptions} <: PERKSingleIntegrator u::uType du::uType u_tmp::uType t::RealT dt::RealT # current time step - dt_cache::RealT # ignored + 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::PERKIntegratorOptions - final_step::Bool # added for convenience + finalstep::Bool # added for convenience # PERK stages: k1::uType k_higher::uType k_s1::uType # Required for third order - t_stage::RealT end # Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 function solve(ode::ODEProblem, alg::PERK3; dt, callback = nothing, kwargs...) u0 = copy(ode.u0) - du = zero(u0) #previously: similar(u0) + du = zero(u0) u_tmp = zero(u0) # PERK stages @@ -250,7 +229,7 @@ function solve(ode::ODEProblem, alg::PERK3; (prob = ode,), ode.f, alg, PERKIntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher, k_s1, t0) + k1, k_higher, k_s1) # initialize callbacks if callback isa CallbackSet @@ -273,10 +252,10 @@ function solve!(integrator::PERK3Integrator) t_end = last(prob.tspan) callbacks = integrator.opts.callback - integrator.final_step = false + integrator.finalstep = false - #@trixi_timeit timer() "main loop" while !integrator.final_step - while !integrator.final_step + #@trixi_timeit timer() "main loop" while !integrator.finalstep + while !integrator.finalstep if isnan(integrator.dt) error("time step size `dt` is NaN") end @@ -289,21 +268,18 @@ function solve!(integrator::PERK3Integrator) end @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin - # k1: + # 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 - # k2 - integrator.t_stage = integrator.t + alg.c[2] * integrator.dt - # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] end - - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + # 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 @@ -316,9 +292,7 @@ function solve!(integrator::PERK3Integrator) end # Higher stages - for stage in 3:(alg.num_stages) - integrator.t_stage = integrator.t + alg.c[stage] * integrator.dt - + for stage in 3:alg.num_stages # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + @@ -326,7 +300,7 @@ function solve!(integrator::PERK3Integrator) alg.a_matrix[stage - 2, 2] * integrator.k_higher[i] end - integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t_stage) + 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 @@ -341,11 +315,8 @@ function solve!(integrator::PERK3Integrator) end @threaded for i in eachindex(integrator.u) - # Original proposed PERK3 - #integrator.u[i] += 0.75 * integrator.k_s1[i] + 0.25 * integrator.k_higher[i] - # Own PERK based on SSPRK33 - integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + - 4.0 * integrator.k_higher[i]) / 6.0 + # "Own" PERK based on SSPRK33 + integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + 4.0 * integrator.k_higher[i]) / 6.0 end end # PERK step timer @@ -362,7 +333,7 @@ function solve!(integrator::PERK3Integrator) end # respect maximum number of iterations - if integrator.iter >= integrator.opts.maxiters && !integrator.final_step + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep @warn "Interrupted. Larger maxiters is needed." terminate!(integrator) end @@ -384,5 +355,4 @@ function Base.resize!(integrator::PERK3Integrator, new_size) resize!(integrator.k_s1, new_size) end -#end # @muladd -end \ No newline at end of file +end # @muladd \ No newline at end of file diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index a4ae0f00ba1..1051988e45c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -1,9 +1,9 @@ -using LinearAlgebra -using Convex +using LinearAlgebra: eigvals using ECOS +using Convex const MOI = Convex.MOI -function filter_eigvals(eig_vals::Array{Complex{Float64}}, threshold:: Float64) +function filter_eigvals(eig_vals, threshold) filtered_eigvals_counter = 0 filtered_eig_vals = Complex{Float64}[] @@ -16,12 +16,12 @@ function filter_eigvals(eig_vals::Array{Complex{Float64}}, threshold:: Float64) end end - println("Total number of $filtered_eigvals_counter eigenvalues has not been recorded because they were smaller than $threshold") + println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they are in magnitude smaller than $threshold \n") return length(filtered_eig_vals), filtered_eig_vals end -function read_in_eig_vals(path_to_eval_file::AbstractString) +function read_in_eig_vals(path_to_eval_file) # Declare and set to some value num_eig_vals = -1 @@ -50,7 +50,7 @@ function read_in_eig_vals(path_to_eval_file::AbstractString) return num_eig_vals, eig_vals end -function polynoms(cons_order::Int, num_stage_evals::Int, num_eig_vals::Int, normalized_powered_eigvals_scaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) +function polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma::Variable) for i in 1:num_eig_vals pnoms[i] = 1.0 end @@ -69,12 +69,9 @@ function polynoms(cons_order::Int, num_stage_evals::Int, num_eig_vals::Int, norm end -function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_max::Float64, dt_eps::Float64, eig_vals::Array{Complex{Float64}}) - +function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, eig_vals) dt_min = 0.0 - dt = -1.0 - abs_p = -1.0 pnoms = ones(Complex{Float64}, num_eig_vals, 1) @@ -93,6 +90,8 @@ function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_ normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) + println("Start optimization of stability polynomial \n") + while dt_max - dt_min > dt_eps dt = 0.5 * (dt_max + dt_min) @@ -120,12 +119,12 @@ function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_ "reltol_inacc" => 5e-5, "nitref" => 9, "maxit" => 100, - "verbose" => 3); silent_solver = false + "verbose" => 3); silent_solver = true ) abs_p = problem.optval - println("Current MaxAbsP: ", abs_p, "\nCurrent dt: ", dt, "\n") + println("MaxAbsP: ", abs_p, "\ndt: ", dt, "\n") if abs_p < 1.0 dt_min = dt @@ -134,10 +133,12 @@ function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_ end end - return evaluate(gamma), abs_p, dt + println("Concluded stability polynomial optimization \n") + + return evaluate(gamma), dt end -function undo_normalization!(cons_order::Int, num_stage_evals::Int, gamma_opt) +function undo_normalization!(cons_order, num_stage_evals, gamma_opt) for k in cons_order + 1:num_stage_evals gamma_opt[k - cons_order] = gamma_opt[k - cons_order] / factorial(k) end From ee9e56978d51c0c994b479feacce99f7f5212689 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Mar 2024 15:13:29 +0100 Subject: [PATCH 16/64] fmt --- .../paired_explicit_runge_kutta.jl | 3 +- .../polynomial_optimizer.jl | 78 +++++++++---------- 2 files changed, 40 insertions(+), 41 deletions(-) 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 9bb3d0b817c..b8d8839b9bf 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 @@ -8,5 +8,4 @@ include("methods_PERK2.jl") include("methods_PERK3.jl") include("polynomial_optimizer.jl") - -end # @muladd \ No newline at end of file +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 index 1051988e45c..e8d6847fb7c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -4,7 +4,6 @@ using Convex const MOI = Convex.MOI function filter_eigvals(eig_vals, threshold) - filtered_eigvals_counter = 0 filtered_eig_vals = Complex{Float64}[] @@ -36,11 +35,11 @@ function read_in_eig_vals(path_to_eval_file) open(path_to_eval_file, "r") do eval_file # Read till end of file - while !eof(eval_file) - + while !eof(eval_file) + # Read a new / next line for every iteration - line_content = readline(eval_file) - + line_content = readline(eval_file) + eig_vals[line_index + 1] = parse(Complex{Float64}, line_content) line_index += 1 @@ -50,32 +49,32 @@ function read_in_eig_vals(path_to_eval_file) return num_eig_vals, eig_vals end -function polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma::Variable) +function polynoms(cons_order, num_stage_evals, num_eig_vals, + normalized_powered_eigvals_scaled, pnoms, gamma::Variable) for i in 1:num_eig_vals pnoms[i] = 1.0 end for k in 1:cons_order for i in 1:num_eig_vals - pnoms[i] += normalized_powered_eigvals_scaled[i,k] + pnoms[i] += normalized_powered_eigvals_scaled[i, k] end end - for k in cons_order + 1:num_stage_evals - pnoms += gamma[k - cons_order] * normalized_powered_eigvals_scaled[:,k] + for k in (cons_order + 1):num_stage_evals + pnoms += gamma[k - cons_order] * normalized_powered_eigvals_scaled[:, k] end - + return maximum(abs(pnoms)) end - function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, eig_vals) dt_min = 0.0 - dt = -1.0 - abs_p = -1.0 + dt = -1.0 + abs_p = -1.0 pnoms = ones(Complex{Float64}, num_eig_vals, 1) - + # Init datastructure for results gamma = Variable(num_stage_evals - cons_order) @@ -84,43 +83,44 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei 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 + normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j end end - normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) + normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, + num_stage_evals) println("Start optimization of stability polynomial \n") while dt_max - dt_min > dt_eps - dt = 0.5 * (dt_max + dt_min) + dt = 0.5 * (dt_max + dt_min) 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] + 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(polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma)) - - Convex.solve!( - problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - #"eps" => 1e13, # 1e13 - "feastol" => 1e-9, # 1e-9 - "abstol" => 1e-9, # 1e-9 - "reltol" => 1e-9, # 1e-9 - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent_solver = true - ) + problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, + normalized_powered_eigvals_scaled, pnoms, gamma)) + + Convex.solve!(problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + #"eps" => 1e13, # 1e13 + "feastol" => 1e-9, # 1e-9 + "abstol" => 1e-9, # 1e-9 + "reltol" => 1e-9, # 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 @@ -135,12 +135,12 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei println("Concluded stability polynomial optimization \n") - return evaluate(gamma), dt + return evaluate(gamma), dt end function undo_normalization!(cons_order, num_stage_evals, gamma_opt) - for k in cons_order + 1:num_stage_evals + for k in (cons_order + 1):num_stage_evals gamma_opt[k - cons_order] = gamma_opt[k - cons_order] / factorial(k) end return gamma_opt -end \ No newline at end of file +end From 32cbb1e2af93f66181dd792606ac07412f720fb0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Mar 2024 15:20:08 +0100 Subject: [PATCH 17/64] fmt --- src/auxiliary/auxiliary.jl | 14 ++--- .../methods_PERK2.jl | 48 +++++++++------- .../methods_PERK3.jl | 55 ++++++++++++------- .../polynomial_optimizer.jl | 42 ++++++-------- 4 files changed, 87 insertions(+), 72 deletions(-) diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 82ed7284e4c..f9d85a3c831 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -346,17 +346,17 @@ function register_error_hints() return nothing end -function read_file(file_path::AbstractString, data_type::Type=Float64) +function read_file(file_path::AbstractString, data_type::Type = Float64) @assert isfile(file_path) "Couldn't find file" data = zeros(data_type, 0) open(file_path, "r") do file - while !eof(file) - line_content = readline(file) - append!(data, parse(data_type, line_content)) - end + while !eof(file) + line_content = readline(file) + append!(data, parse(data_type, line_content)) + end end num_lines = length(data) - + return num_lines, data - end +end end # @muladd 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 ebff168145a..8243a2c7362 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -3,6 +3,7 @@ # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin +#! format: noindent # Abstract base type for both single/standalone and multi-level P-ERK time integration schemes abstract type PERK end @@ -22,7 +23,8 @@ function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) return reverse(a_coeffs) end -function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, bS, c_end) +function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, + bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -50,7 +52,7 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, - eig_vals) + eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) num_mon_coeffs = length(mon_coeffs) @@ -67,7 +69,8 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat return a_matrix, c end -function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, bS, c_end) +function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, + bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) @@ -118,11 +121,12 @@ mutable struct PERK2 <: PERKSingle c_end::Float64 #Constructor that read the coefficients from the file - function PERK2(num_stages, base_path_mon_coeffs::AbstractString, bS = 1.0, c_end = 0.5) + function PERK2(num_stages, base_path_mon_coeffs::AbstractString, bS = 1.0, + c_end = 0.5) newPERK2 = new(num_stages) newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, - base_path_mon_coeffs, + base_path_mon_coeffs, bS, c_end) newPERK2.b1 = one(bS) - bS @@ -135,7 +139,8 @@ mutable struct PERK2 <: PERKSingle function PERK2(num_stages, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5) newPERK2 = new(num_stages) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, semi, + newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, + semi, bS, c_end) newPERK2.b1 = one(bS) - bS @@ -155,7 +160,8 @@ mutable struct PERKIntegratorOptions{Callback} end function PERKIntegratorOptions(callback, tspan; maxiters = typemax(Int), kwargs...) - PERKIntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, [last(tspan)]) + PERKIntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, + [last(tspan)]) end abstract type PERKIntegrator end @@ -166,7 +172,7 @@ abstract type PERKSingleIntegrator <: PERKIntegrator end # https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1 # which are used in Trixi. mutable struct PERK2Integrator{RealT <: Real, uType, Params, Sol, F, Alg, - PERKIntegratorOptions} <: PERKSingleIntegrator + PERKIntegratorOptions} <: PERKSingleIntegrator u::uType du::uType u_tmp::uType @@ -209,10 +215,10 @@ function solve(ode::ODEProblem, alg::PERK2; iter = 0 integrator = PERK2Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, - (prob = ode,), ode.f, alg, - PERKIntegratorOptions(callback, ode.tspan; kwargs...), - false, - k1, k_higher) + (prob = ode,), ode.f, alg, + PERKIntegratorOptions(callback, ode.tspan; kwargs...), + false, + k1, k_higher) # initialize callbacks if callback isa CallbackSet @@ -261,7 +267,8 @@ function solve!(integrator::PERK2Integrator) 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) + 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 @@ -272,11 +279,14 @@ function solve!(integrator::PERK2Integrator) # 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] + 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) + 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 @@ -284,7 +294,8 @@ function solve!(integrator::PERK2Integrator) end @threaded for i in eachindex(integrator.u) - integrator.u[i] += alg.b1 * integrator.k1[i] + alg.bS * integrator.k_higher[i] + integrator.u[i] += alg.b1 * integrator.k1[i] + + alg.bS * integrator.k_higher[i] end end # PERK2 step @@ -343,5 +354,4 @@ function Base.resize!(integrator::PERK2Integrator, new_size) resize!(integrator.k1, new_size) resize!(integrator.k_higher, new_size) end - -end # @muladd \ No newline at end of file +end # @muladd diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 69e119f3a27..e4203408d00 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -5,6 +5,7 @@ using NLsolve: nlsolve @muladd begin +#! format: noindent # Initialize Butcher array abscissae c for PERK3 based on SSPRK33 base method function c_PERK3_SSP33(num_stages, c_S2) @@ -22,7 +23,8 @@ function c_PERK3_SSP33(num_stages, c_S2) return c end -function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, mon_coeffs, c_S2) +function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, + mon_coeffs, c_S2) c_ts = c_PERK3_SSP33(num_stages, c_S2) # ts = timestep c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_S2 is equal to 1 @@ -34,8 +36,8 @@ function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, term1 *= a_unknown[num_stage_evals - 1 - j] term2 *= a_unknown[num_stage_evals - j] end - term1 *= c_ts[num_stages - 2 - i] * 1/6 - term2 *= c_ts[num_stages - 1 - i] * 4/6 + term1 *= c_ts[num_stages - 2 - i] * 1 / 6 + term2 *= c_ts[num_stages - 1 - i] * 4 / 6 c_eq[i] = mon_coeffs[i] - (term1 + term2) end @@ -49,13 +51,14 @@ function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, term2 *= c_ts[num_stages - 1 - i] * 4 / 6 c_eq[i] = mon_coeffs[i] - term2 - - c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - a_unknown[num_stage_evals - 1] + c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - + a_unknown[num_stage_evals - 1] return c_eq end -function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, c_S2) +function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, + c_S2) # Initialize array of c c = c_PERK3_SSP33(num_stages, c_S2) @@ -79,12 +82,14 @@ function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat eig_vals = eigvals(J) num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, eig_vals) + mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, + dt_eps, eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) # Define the objective_function function objective_function(x) - return PERK3_Butcher_tableau_objective(x, num_stages, num_stages, mon_coeffs, c_S2) + return PERK3_Butcher_tableau_objective(x, num_stages, num_stages, + mon_coeffs, c_S2) end # Call nlsolver to solve repeatedly until the result is not NaN or negative values @@ -117,7 +122,8 @@ function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat return a_matrix, c end -function compute_PERK3_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, c_S2) +function compute_PERK3_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, + c_S2) # Initialize array of c c = c_PERK3_SSP33(num_stages, c_S2) @@ -132,7 +138,8 @@ function compute_PERK3_Butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac a_matrix = zeros(coeffs_max, 2) a_matrix[:, 1] = c[3:end] - path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * string(num_stages) * ".txt" + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * + string(num_stages) * ".txt" num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) @assert num_mon_coeffs == coeffs_max @@ -161,7 +168,8 @@ mutable struct PERK3 <: PERKSingle c::Vector{Float64} # Constructor for previously computed A Coeffs - function PERK3(num_stages, base_path_mon_coeffs::AbstractString, c_S2::Float64 = 1.0) + function PERK3(num_stages, base_path_mon_coeffs::AbstractString, + c_S2 = 1.0) newPERK3 = new(num_stages) newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, @@ -172,11 +180,12 @@ mutable struct PERK3 <: PERKSingle end # Constructor that computes Butcher matrix A coefficients from a semidiscretization - function PERK3(num_stages, semi::AbstractSemidiscretization, c_S2::Float64 = 1.0) + function PERK3(num_stages, semi::AbstractSemidiscretization, + c_S2 = 1.0) newPERK3 = new(num_stages) newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, - semi, + semi, c_S2) display(newPERK3.a_matrix) @@ -279,7 +288,8 @@ function solve!(integrator::PERK3Integrator) 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) + 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 @@ -292,15 +302,18 @@ function solve!(integrator::PERK3Integrator) end # Higher stages - for stage in 3:alg.num_stages + for stage in 3:(alg.num_stages) # 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] + 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) + 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 @@ -316,7 +329,8 @@ function solve!(integrator::PERK3Integrator) @threaded for i in eachindex(integrator.u) # "Own" PERK based on SSPRK33 - integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + 4.0 * integrator.k_higher[i]) / 6.0 + integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + + 4.0 * integrator.k_higher[i]) / 6.0 end end # PERK step timer @@ -354,5 +368,4 @@ function Base.resize!(integrator::PERK3Integrator, new_size) resize!(integrator.k_higher, new_size) resize!(integrator.k_s1, new_size) end - -end # @muladd \ No newline at end of file +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 index e8d6847fb7c..9b4d59d42cb 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -3,45 +3,20 @@ using ECOS using Convex const MOI = Convex.MOI -function filter_eigvals(eig_vals, threshold) - filtered_eigvals_counter = 0 - filtered_eig_vals = Complex{Float64}[] - - for eig_val in eig_vals - if abs(eig_val) < threshold - filtered_eigvals_counter += 1 - else - push!(filtered_eig_vals, eig_val) - end - end - - println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they are in magnitude smaller than $threshold \n") - - return length(filtered_eig_vals), filtered_eig_vals -end - function read_in_eig_vals(path_to_eval_file) - # Declare and set to some value num_eig_vals = -1 - open(path_to_eval_file, "r") do eval_file num_eig_vals = countlines(eval_file) end - eig_vals = Array{Complex{Float64}}(undef, num_eig_vals) - line_index = 0 - open(path_to_eval_file, "r") do eval_file # Read till end of file while !eof(eval_file) - # Read a new / next line for every iteration line_content = readline(eval_file) - eig_vals[line_index + 1] = parse(Complex{Float64}, line_content) - line_index += 1 end end @@ -49,6 +24,23 @@ function read_in_eig_vals(path_to_eval_file) return num_eig_vals, eig_vals end +function filter_eigvals(eig_vals, threshold) + filtered_eigvals_counter = 0 + filtered_eig_vals = Complex{Float64}[] + + for eig_val in eig_vals + if abs(eig_val) < threshold + filtered_eigvals_counter += 1 + else + push!(filtered_eig_vals, eig_val) + end + end + + println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they are in magnitude smaller than $threshold \n") + + return length(filtered_eig_vals), filtered_eig_vals +end + function polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma::Variable) for i in 1:num_eig_vals From e3ec50cb1bfd630505b20ab4fbb6197636d18c18 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Mar 2024 15:25:51 +0100 Subject: [PATCH 18/64] semantic ordering --- .../paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 b8d8839b9bf..9ac1941f987 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 @@ -5,7 +5,8 @@ @muladd begin #! format: noindent +include("polynomial_optimizer.jl") + include("methods_PERK2.jl") include("methods_PERK3.jl") -include("polynomial_optimizer.jl") end # @muladd From fbcd6e3047f22f44582706b55ad6200448cbe31c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Mar 2024 08:53:25 +0100 Subject: [PATCH 19/64] add literature --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 6 +++++- .../paired_explicit_runge_kutta/methods_PERK3.jl | 6 +++++- 2 files changed, 10 insertions(+), 2 deletions(-) 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 8243a2c7362..59d25e05892 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -107,8 +107,12 @@ end PERK2() The following structures and methods provide a minimal implementation of -the paired explicit Runge-Kutta method (https://doi.org/10.1016/j.jcp.2019.05.014) +the second-order paired explicit Runge-Kutta method optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +- 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 PERK2 <: PERKSingle diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index e4203408d00..bed2ca3c5e4 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -157,8 +157,12 @@ end PERK3() The following structures and methods provide a minimal implementation of -the third order paired explicit Runge-Kutta method (https://www.sciencedirect.com/science/article/pii/S0021999122005320) +the third-order paired explicit Runge-Kutta method optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +- 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) """ mutable struct PERK3 <: PERKSingle From b1b6d7cbaadaeeae500aa4d5bdb47360e8627856 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 6 Mar 2024 18:05:41 +0100 Subject: [PATCH 20/64] Strip code of p = 3 --- .../methods_PERK3.jl | 375 ------------------ .../paired_explicit_runge_kutta.jl | 1 - 2 files changed, 376 deletions(-) delete mode 100644 src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl deleted file mode 100644 index bed2ca3c5e4..00000000000 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ /dev/null @@ -1,375 +0,0 @@ -# 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 NLsolve: nlsolve - -@muladd begin -#! format: noindent - -# Initialize Butcher array abscissae c for PERK3 based on SSPRK33 base method -function c_PERK3_SSP33(num_stages, c_S2) - c = zeros(num_stages) - - # Last timesteps as for SSPRK33 - c[num_stages] = 0.5 - c[num_stages - 1] = 1 - - # Linear increasing timestep for remainder - for i in 2:(num_stages - 2) - c[i] = c_S2 * (i - 1) / (num_stages - 3) - end - - return c -end - -function PERK3_Butcher_tableau_objective(a_unknown, num_stages, num_stage_evals, - mon_coeffs, c_S2) - c_ts = c_PERK3_SSP33(num_stages, c_S2) # ts = timestep - - c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_S2 is equal to 1 - # Both terms should be present - for i in 1:(num_stage_evals - 4) - term1 = a_unknown[num_stage_evals - 1] - term2 = a_unknown[num_stage_evals] - for j in 1:i - term1 *= a_unknown[num_stage_evals - 1 - j] - term2 *= a_unknown[num_stage_evals - j] - end - term1 *= c_ts[num_stages - 2 - i] * 1 / 6 - term2 *= c_ts[num_stages - 1 - i] * 4 / 6 - - c_eq[i] = mon_coeffs[i] - (term1 + term2) - end - - # Highest coefficient: Only one term present - i = num_stage_evals - 3 - term2 = a_unknown[num_stage_evals] - for j in 1:i - term2 *= a_unknown[num_stage_evals - j] - end - term2 *= c_ts[num_stages - 1 - i] * 4 / 6 - - c_eq[i] = mon_coeffs[i] - term2 - c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - - a_unknown[num_stage_evals - 1] - - return c_eq -end - -function compute_PERK3_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, - c_S2) - - # Initialize array of c - c = c_PERK3_SSP33(num_stages, c_S2) - - # Initialize the array of our solution - a_unknown = zeros(num_stages) - - # Special case of e = 3 - if num_stages == 3 - a_unknown = [0, c[2], 0.25] - - else - # Calculate coefficients of the stability polynomial in monomial form - cons_order = 3 - dtmax = 1.0 - dt_eps = 1e-9 - filter_threshold = 1e-12 - - # Compute spectrum - J = jacobian_ad_forward(semi) - eig_vals = eigvals(J) - num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - - mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, - dt_eps, eig_vals) - mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) - - # Define the objective_function - function objective_function(x) - return PERK3_Butcher_tableau_objective(x, num_stages, num_stages, - mon_coeffs, c_S2) - end - - # Call nlsolver to solve repeatedly until the result is not NaN or negative values - is_sol_valid = false - while !is_sol_valid - # Initialize initial guess - x0 = 0.1 .* rand(num_stages) - x0[1] = 0.0 - x0[2] = c[2] - - sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 4e-16, - iterations = 10^4, xtol = 1e-13) - - a_unknown = sol.zero - - # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 - is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && - all(x -> !isnan(x) && x >= 0, c[3:end] .- a_unknown[3:end]) - end - end - - println("a_unknown") - println(a_unknown[3:end]) # To debug - - a_matrix = zeros(num_stages - 2, 2) - a_matrix[:, 1] = c[3:end] - a_matrix[:, 1] -= a_unknown[3:end] - a_matrix[:, 2] = a_unknown[3:end] - - return a_matrix, c -end - -function compute_PERK3_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, - c_S2) - - # Initialize array of c - c = c_PERK3_SSP33(num_stages, c_S2) - - println("Timestep-split: ") - display(c) - println("\n") - - # - 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_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * - string(num_stages) * ".txt" - num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) - @assert num_mon_coeffs == coeffs_max - - a_matrix[:, 1] -= A - a_matrix[:, 2] = A - - println("A matrix: ") - display(a_matrix) - println() - - return a_matrix, c -end - -""" - PERK3() - -The following structures and methods provide a minimal implementation of -the third-order paired explicit Runge-Kutta method -optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). - -- 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) -""" - -mutable struct PERK3 <: PERKSingle - const num_stages::Int - - a_matrix::Matrix{Float64} - c::Vector{Float64} - - # Constructor for previously computed A Coeffs - function PERK3(num_stages, base_path_mon_coeffs::AbstractString, - c_S2 = 1.0) - newPERK3 = new(num_stages) - - newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, - base_path_mon_coeffs, - c_S2) - - return newPERK3 - end - - # Constructor that computes Butcher matrix A coefficients from a semidiscretization - function PERK3(num_stages, semi::AbstractSemidiscretization, - c_S2 = 1.0) - newPERK3 = new(num_stages) - - newPERK3.a_matrix, newPERK3.c = compute_PERK3_Butcher_tableau(num_stages, - semi, - c_S2) - - display(newPERK3.a_matrix) - - return newPERK3 - end -end # struct PERK3 - -# 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 PERK3Integrator{RealT <: Real, uType, Params, Sol, F, Alg, - PERKIntegratorOptions} <: PERKSingleIntegrator - 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::PERKIntegratorOptions - finalstep::Bool # added for convenience - # PERK stages: - k1::uType - k_higher::uType - k_s1::uType # Required for third order -end - -# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 -function solve(ode::ODEProblem, alg::PERK3; - dt, callback = nothing, kwargs...) - u0 = copy(ode.u0) - du = zero(u0) - u_tmp = zero(u0) - - # PERK stages - k1 = zero(u0) - k_higher = zero(u0) - k_s1 = zero(u0) - - t0 = first(ode.tspan) - iter = 0 - - integrator = PERK3Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p, - (prob = ode,), ode.f, alg, - PERKIntegratorOptions(callback, ode.tspan; kwargs...), - false, - k1, k_higher, k_s1) - - # 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 - - solve!(integrator) -end - -function solve!(integrator::PERK3Integrator) - @unpack prob = integrator.sol - @unpack alg = integrator - t_end = last(prob.tspan) - callbacks = integrator.opts.callback - - integrator.finalstep = false - - #@trixi_timeit timer() "main loop" while !integrator.finalstep - while !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.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 - - if alg.num_stages == 3 - @threaded for i in eachindex(integrator.du) - integrator.k_s1[i] = integrator.k_higher[i] - end - end - - # Higher stages - for stage in 3:(alg.num_stages) - # 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 - - # IDEA: Stop for loop at num_stages -1 to avoid if (maybe more performant?) - if stage == alg.num_stages - 1 - @threaded for i in eachindex(integrator.du) - integrator.k_s1[i] = integrator.k_higher[i] - end - end - end - - @threaded for i in eachindex(integrator.u) - # "Own" PERK based on SSPRK33 - integrator.u[i] += (integrator.k1[i] + integrator.k_s1[i] + - 4.0 * integrator.k_higher[i]) / 6.0 - end - end # PERK 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 # "main loop" timer - - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) -end - -# used for AMR (Adaptive Mesh Refinement) -function Base.resize!(integrator::PERK3Integrator, 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) - resize!(integrator.k_s1, 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 9ac1941f987..df9b1b6c7cf 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 @@ -8,5 +8,4 @@ include("polynomial_optimizer.jl") include("methods_PERK2.jl") -include("methods_PERK3.jl") end # @muladd From 33bb8da25e0eabd9e87ee511249f7e9dd6584515 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 7 Mar 2024 15:08:25 +0100 Subject: [PATCH 21/64] Add the line to show the error of the elixir --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 8f12bfb328f..1cd128de831 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -53,4 +53,7 @@ sol = Trixi.solve(ode, ode_algorithm, save_everystep = false, callback = callbacks); # Print the timer summary -summary_callback() \ No newline at end of file +summary_callback() + +# Show the error of the elixir +analysis_callback(sol) \ No newline at end of file From 5f763bbacbef00bc424fc33eca7a51aa8b875d0d Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 8 Mar 2024 10:37:40 +0100 Subject: [PATCH 22/64] Make adjustments to a test file and delete example of PERK3 --- .../elixir_burgers_PERK3.jl | 53 ------------------- test/test_tree_1d_advection.jl | 14 +++++ 2 files changed, 14 insertions(+), 53 deletions(-) delete mode 100644 examples/structured_1d_dgsem/elixir_burgers_PERK3.jl diff --git a/examples/structured_1d_dgsem/elixir_burgers_PERK3.jl b/examples/structured_1d_dgsem/elixir_burgers_PERK3.jl deleted file mode 100644 index 056f83a3e99..00000000000 --- a/examples/structured_1d_dgsem/elixir_burgers_PERK3.jl +++ /dev/null @@ -1,53 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the (inviscid) Burgers' equation - -equations = InviscidBurgersEquation1D() - -initial_condition = initial_condition_convergence_test - -# Create DG solver with polynomial degree = 3 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,) - -# Create curved mesh with 16 cells -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) - -stepsize_callback = StepsizeCallback(cfl = 3.0) - -callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, - stepsize_callback) - -############################################################################### -# run the simulation - -ode_algorithm = PERK3(8, 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() \ No newline at end of file diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index a580a3b5600..d00bc82deeb 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.011662300515980219], + linf=[0.01647256923710194]) + # 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)) < 1000 + end +end end end # module From 9149b7b97d3d0f8dd0f42f3e140512a614eaa9b5 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:39:04 +0100 Subject: [PATCH 23/64] Update Project.toml Co-authored-by: Daniel Doehring --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index d55d741435f..e51a245ef31 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,6 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" From 6c04331866b372c1cfdc4b052177a7bd1f2772de Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:39:21 +0100 Subject: [PATCH 24/64] Update examples/tree_1d_dgsem/elixir_advection_PERK2.jl Co-authored-by: Daniel Doehring --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 1cd128de831..e30eb44370d 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -37,7 +37,7 @@ summary_callback = SummaryCallback() analysis_callback = AnalysisCallback(semi, interval = 100) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 2.5/4) +stepsize_callback = StepsizeCallback(cfl = 2.5) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, analysis_callback, From 0bafb7f9fccc60ab3c190da3a1d621ab68951081 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:39:52 +0100 Subject: [PATCH 25/64] Add comments in an example of PERK2 Co-authored-by: Daniel Doehring --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index e30eb44370d..680ed2da916 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -46,6 +46,8 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation +# Construct second order P-ERK method with 6 stages for +# given simulation setup ode_algorithm = PERK2(6, semi) sol = Trixi.solve(ode, ode_algorithm, From 64926fbb2abf7f6cba485cea47bf93951ee11fea Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:40:02 +0100 Subject: [PATCH 26/64] Update src/Trixi.jl Co-authored-by: Daniel Doehring --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 15c95963674..334aca03208 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -276,7 +276,7 @@ export trixi_include, examples_dir, get_examples, default_example, export ode_norm, ode_unstable_check -export PERK2, PERK3 +export PERK2 export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure From 6ce99e61e21f4336a51a0f86c6a2a9e0a4837446 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:40:18 +0100 Subject: [PATCH 27/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 3 --- 1 file changed, 3 deletions(-) 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 59d25e05892..8b4fe77a6a3 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -31,9 +31,6 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat for k in 2:num_stages c[k] = c_end * (k - 1) / (num_stages - 1) end - println("Timestep-split: ") - display(c) - println("\n") se_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) From 50cf650a866275c8b332f5f91c45600874977219 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:42:15 +0100 Subject: [PATCH 28/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 3 --- 1 file changed, 3 deletions(-) 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 8b4fe77a6a3..c6bc27fa7ad 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -59,9 +59,6 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat a_matrix[:, 1] -= A a_matrix[:, 2] = A - println("A matrix: ") - display(a_matrix) - println() return a_matrix, c end From f48cf85770a9b1d9ba85c8b7ef496fd2ed97ce28 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:42:34 +0100 Subject: [PATCH 29/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 3 --- 1 file changed, 3 deletions(-) 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 c6bc27fa7ad..0a4f5ffd7f3 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -71,9 +71,6 @@ function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac for k in 2:num_stages c[k] = c_end * (k - 1) / (num_stages - 1) end - println("Timestep-split: ") - display(c) - println("\n") se_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) From e7060e56c3a8ad45079d7f275566137a980865b9 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:43:33 +0100 Subject: [PATCH 30/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 3 --- 1 file changed, 3 deletions(-) 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 0a4f5ffd7f3..68f1a7303b4 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -87,9 +87,6 @@ function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac a_matrix[:, 1] -= A a_matrix[:, 2] = A - println("A matrix: ") - display(a_matrix) - println() return a_matrix, c end From 88e81682d8f94a1098fc2646a9cd8ef15c14479d Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 8 Mar 2024 13:28:02 +0100 Subject: [PATCH 31/64] add tspan as a parameter and adjust the elixir accordingly --- .../tree_1d_dgsem/elixir_advection_PERK2.jl | 7 ++-- .../methods_PERK2.jl | 15 ++++---- .../polynomial_optimizer.jl | 35 ++++--------------- 3 files changed, 18 insertions(+), 39 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 680ed2da916..85ce6167f1a 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -46,9 +46,12 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation +# Define tspan to calculate maximum time step allowed for the bisection algorithm used in calculate polynomial coefficients in ODE algorithm +tspan = [0.0, 1.0] + # Construct second order P-ERK method with 6 stages for # given simulation setup -ode_algorithm = PERK2(6, semi) +ode_algorithm = PERK2(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 @@ -58,4 +61,4 @@ sol = Trixi.solve(ode, ode_algorithm, summary_callback() # Show the error of the elixir -analysis_callback(sol) \ No newline at end of file +analysis_callback(sol) 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 68f1a7303b4..50183e4573a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -23,7 +23,7 @@ function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) return reverse(a_coeffs) end -function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, +function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, tspan, bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) @@ -41,14 +41,14 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat cons_order = 2 filter_thres = 1e-12 - dtmax = 1.0 - dt_eps = 1e-9 + dtmax = tspan[2] - tspan[1] + dteps = 1e-9 eig_vals = eigvals(jacobian_ad_forward(semi)) num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) - mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps, + mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dteps, eig_vals) mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) @@ -59,7 +59,6 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat a_matrix[:, 1] -= A a_matrix[:, 2] = A - return a_matrix, c end @@ -87,7 +86,6 @@ function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac a_matrix[:, 1] -= A a_matrix[:, 2] = A - return a_matrix, c end @@ -128,11 +126,12 @@ mutable struct PERK2 <: PERKSingle end #Constructor that calculate the coefficients with polynomial optimizer - function PERK2(num_stages, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5) + function PERK2(num_stages, tspan, semi::AbstractSemidiscretization, bS = 1.0, + c_end = 0.5) newPERK2 = new(num_stages) newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, - semi, + semi, tspan, bS, c_end) newPERK2.b1 = one(bS) - bS diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 9b4d59d42cb..09f6a612db6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -3,27 +3,6 @@ using ECOS using Convex const MOI = Convex.MOI -function read_in_eig_vals(path_to_eval_file) - # Declare and set to some value - num_eig_vals = -1 - open(path_to_eval_file, "r") do eval_file - num_eig_vals = countlines(eval_file) - end - eig_vals = Array{Complex{Float64}}(undef, num_eig_vals) - line_index = 0 - open(path_to_eval_file, "r") do eval_file - # Read till end of file - while !eof(eval_file) - # Read a new / next line for every iteration - line_content = readline(eval_file) - eig_vals[line_index + 1] = parse(Complex{Float64}, line_content) - line_index += 1 - end - end - - return num_eig_vals, eig_vals -end - function filter_eigvals(eig_vals, threshold) filtered_eigvals_counter = 0 filtered_eig_vals = Complex{Float64}[] @@ -60,8 +39,8 @@ function polynoms(cons_order, num_stage_evals, num_eig_vals, return maximum(abs(pnoms)) end -function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, eig_vals) - dt_min = 0.0 +function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals) + dtmin = 0.0 dt = -1.0 abs_p = -1.0 @@ -84,8 +63,8 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei println("Start optimization of stability polynomial \n") - while dt_max - dt_min > dt_eps - dt = 0.5 * (dt_max + dt_min) + while dtmax - dtmin > dteps + dt = 0.5 * (dtmax + dtmin) for k in 1:num_stage_evals dt_k = dt^k @@ -116,12 +95,10 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei abs_p = problem.optval - println("MaxAbsP: ", abs_p, "\ndt: ", dt, "\n") - if abs_p < 1.0 - dt_min = dt + dtmin = dt else - dt_max = dt + dtmax = dt end end From bb938d9d9f64db6336eebdc2be5e778cc537855f Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 8 Mar 2024 15:56:27 +0100 Subject: [PATCH 32/64] add an additional constructor and modify the function compute_PERK2_butcher_tableau --- .../methods_PERK2.jl | 32 +++++++++++++------ .../polynomial_optimizer.jl | 2 ++ 2 files changed, 25 insertions(+), 9 deletions(-) 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 50183e4573a..8c0cf83d159 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -23,7 +23,7 @@ function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) return reverse(a_coeffs) end -function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, tspan, +function compute_PERK2_butcher_tableau(num_stages, eig_vals, tspan, bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) @@ -44,8 +44,6 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat dtmax = tspan[2] - tspan[1] dteps = 1e-9 - eig_vals = eigvals(jacobian_ad_forward(semi)) - num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres) mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dteps, @@ -62,7 +60,7 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat return a_matrix, c end -function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, +function compute_PERK2_butcher_tableau(num_stages, base_path_mon_coeffs::AbstractString, bS, c_end) # c Vector form Butcher Tableau (defines timestep per stage) @@ -110,12 +108,12 @@ mutable struct PERK2 <: PERKSingle b1::Float64 c_end::Float64 - #Constructor that read the coefficients from the file + # Constructor that read the coefficients from the file function PERK2(num_stages, base_path_mon_coeffs::AbstractString, bS = 1.0, c_end = 0.5) newPERK2 = new(num_stages) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages, base_path_mon_coeffs, bS, c_end) @@ -125,13 +123,29 @@ mutable struct PERK2 <: PERKSingle return newPERK2 end - #Constructor that calculate the coefficients with polynomial optimizer + # Constructor that calculate the coefficients with polynomial optimizer from a semidiscretization function PERK2(num_stages, tspan, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5) + eig_vals = eigvals(jacobian_ad_forward(semi)) + newPERK2 = new(num_stages) + + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages, + eig_vals, tspan, + bS, c_end) + + newPERK2.b1 = one(bS) - bS + newPERK2.bS = bS + newPERK2.c_end = c_end + return newPERK2 + end + + # Constructor that calculate the coefficients with polynomial optimizer from a list of eigenvalues + function PERK2(num_stages, tspan, eig_vals::Vector{ComplexF64}, bS = 1.0, + c_end = 0.5) newPERK2 = new(num_stages) - newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages, - semi, tspan, + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages, + eig_vals, tspan, bS, c_end) newPERK2.b1 = one(bS) - bS diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 09f6a612db6..f1be279533a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -102,6 +102,8 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_ end end + println("MaxAbsP: ", abs_p, "\ndt: ", dt, "\n") + println("Concluded stability polynomial optimization \n") return evaluate(gamma), dt From 38f4fcb8432c1cec33cc07e0bb8c874cfaff153c Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Mon, 11 Mar 2024 17:59:45 +0100 Subject: [PATCH 33/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 8c0cf83d159..02818a2e8f7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -11,7 +11,7 @@ abstract type PERK end abstract type PERKSingle <: PERK end function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) - a_coeffs = mon_coeffs + a_coeffs = copy(mon_coeffs) for stage in 1:(num_stage_evals - 2) a_coeffs[stage] /= se_factors[stage] From 7619b2a4168a340c836dc8db054e10c4d8a84e9c Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 12 Mar 2024 09:17:05 +0100 Subject: [PATCH 34/64] update filter function in polynomial optimizer, add literature, and change se_factors to bc_factors --- .../methods_PERK2.jl | 12 ++++---- .../polynomial_optimizer.jl | 29 +++++++++++++++---- 2 files changed, 29 insertions(+), 12 deletions(-) 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 8c0cf83d159..d763989a3a0 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -10,11 +10,11 @@ abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end -function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs) +function compute_a_coeffs(num_stage_evals, bc_factors, mon_coeffs) a_coeffs = mon_coeffs for stage in 1:(num_stage_evals - 2) - a_coeffs[stage] /= se_factors[stage] + a_coeffs[stage] /= bc_factors[stage] for prev_stage in 1:(stage - 1) a_coeffs[stage] /= a_coeffs[prev_stage] end @@ -31,7 +31,7 @@ function compute_PERK2_butcher_tableau(num_stages, eig_vals, tspan, for k in 2:num_stages c[k] = c_end * (k - 1) / (num_stages - 1) end - se_factors = bS * reverse(c[2:(end - 1)]) + bc_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 @@ -52,7 +52,7 @@ function compute_PERK2_butcher_tableau(num_stages, eig_vals, tspan, num_mon_coeffs = length(mon_coeffs) @assert num_mon_coeffs == coeffs_max - A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) + A = compute_a_coeffs(num_stages, bc_factors, mon_coeffs) a_matrix[:, 1] -= A a_matrix[:, 2] = A @@ -68,7 +68,7 @@ function compute_PERK2_butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac for k in 2:num_stages c[k] = c_end * (k - 1) / (num_stages - 1) end - se_factors = bS * reverse(c[2:(end - 1)]) + bc_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 @@ -79,7 +79,7 @@ function compute_PERK2_butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac path_mon_coeffs = base_path_mon_coeffs * "gamma_" * string(num_stages) * ".txt" num_mon_coeffs, mon_coeffs = read_file(path_mon_coeffs, Float64) @assert num_mon_coeffs == coeffs_max - A = compute_a_coeffs(num_stages, se_factors, mon_coeffs) + A = compute_a_coeffs(num_stages, bc_factors, mon_coeffs) a_matrix[:, 1] -= A a_matrix[:, 2] = A diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index f1be279533a..c380bdcd557 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -8,20 +8,23 @@ function filter_eigvals(eig_vals, threshold) filtered_eig_vals = Complex{Float64}[] for eig_val in eig_vals - if abs(eig_val) < threshold + # Filter out eigenvalues with positive real parts, those with negative imaginary parts due to eigenvalues' symmetry + # around real axis, or the eigenvalues that are too small. + if real(eig_val) > 0 || imag(eig_val) < 0 || abs(eig_val) < threshold filtered_eigvals_counter += 1 else push!(filtered_eig_vals, eig_val) end end - println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they are in magnitude smaller than $threshold \n") + println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they either are in magnitude smaller than $threshold, \n + have positive real parts, or have negative imaginary parts") return length(filtered_eig_vals), filtered_eig_vals end -function polynoms(cons_order, num_stage_evals, num_eig_vals, - normalized_powered_eigvals_scaled, pnoms, gamma::Variable) +function stability_polynomials(cons_order, num_stage_evals, num_eig_vals, + normalized_powered_eigvals_scaled, pnoms, gamma::Variable) for i in 1:num_eig_vals pnoms[i] = 1.0 end @@ -39,6 +42,19 @@ function polynoms(cons_order, num_stage_evals, num_eig_vals, return maximum(abs(pnoms)) end +""" + bisection() + +The following structures and methods provide a simplified implementation to +discover optimally stable polynomial approximations of the exponential function. +These are designed for the one-step 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) +""" + function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals) dtmin = 0.0 dt = -1.0 @@ -75,8 +91,9 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_ end # Use last optimal values for gamma in (potentially) next iteration - problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, - normalized_powered_eigvals_scaled, pnoms, gamma)) + problem = minimize(stability_polynomials(cons_order, num_stage_evals, num_eig_vals, + normalized_powered_eigvals_scaled, pnoms, + gamma)) Convex.solve!(problem, # Parameters taken from default values for EiCOS From f2d9fa8083e6f78ba495289b4814ee82fcddeac3 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Tue, 12 Mar 2024 13:37:35 +0100 Subject: [PATCH 35/64] Apply suggestions from code review Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 6 +++--- .../polynomial_optimizer.jl | 13 ++++++++----- 2 files changed, 11 insertions(+), 8 deletions(-) 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 0c8f7406fe8..6dd9957a9b2 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -108,7 +108,7 @@ mutable struct PERK2 <: PERKSingle b1::Float64 c_end::Float64 - # Constructor that read the coefficients from the file + # Constructor that reads the coefficients from a file function PERK2(num_stages, base_path_mon_coeffs::AbstractString, bS = 1.0, c_end = 0.5) newPERK2 = new(num_stages) @@ -123,7 +123,7 @@ mutable struct PERK2 <: PERKSingle return newPERK2 end - # Constructor that calculate the coefficients with polynomial optimizer from a semidiscretization + # Constructor that calculates the coefficients with polynomial optimizer from a semidiscretization function PERK2(num_stages, tspan, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5) eig_vals = eigvals(jacobian_ad_forward(semi)) @@ -139,7 +139,7 @@ mutable struct PERK2 <: PERKSingle return newPERK2 end - # Constructor that calculate the coefficients with polynomial optimizer from a list of eigenvalues + # Constructor that calculates the coefficients with polynomial optimizer from a list of eigenvalues function PERK2(num_stages, tspan, eig_vals::Vector{ComplexF64}, bS = 1.0, c_end = 0.5) newPERK2 = new(num_stages) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index c380bdcd557..e60912254a7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -25,20 +25,24 @@ end function stability_polynomials(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma::Variable) + # Initialize with zero'th order (z^0) coefficient for i in 1:num_eig_vals pnoms[i] = 1.0 end + # First `cons_order` terms of the exponential for k in 1:cons_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 (cons_order + 1):num_stage_evals pnoms += gamma[k - cons_order] * normalized_powered_eigvals_scaled[:, k] end + # For optimization only the maximum is relevant return maximum(abs(pnoms)) end @@ -46,8 +50,8 @@ end bisection() The following structures and methods provide a simplified implementation to -discover optimally stable polynomial approximations of the exponential function. -These are designed for the one-step integration of initial value ordinary +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). @@ -60,6 +64,7 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_ dt = -1.0 abs_p = -1.0 + # Construct stability polynomial for each eigenvalue pnoms = ones(Complex{Float64}, num_eig_vals, 1) # Init datastructure for results @@ -82,6 +87,7 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_ 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 @@ -118,9 +124,6 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_ dtmax = dt end end - - println("MaxAbsP: ", abs_p, "\ndt: ", dt, "\n") - println("Concluded stability polynomial optimization \n") return evaluate(gamma), dt From ccf4260992e8d1169d0665560b4a2ae021625efe Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 12 Mar 2024 14:10:55 +0100 Subject: [PATCH 36/64] change tspan to have (,) instead of[,]. filter_eigenvalues minor adjustments --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 85ce6167f1a..f8c20adbdb9 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -47,7 +47,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, # run the simulation # Define tspan to calculate maximum time step allowed for the bisection algorithm used in calculate polynomial coefficients in ODE algorithm -tspan = [0.0, 1.0] +tspan = (0.0, 1.0) # Construct second order P-ERK method with 6 stages for # given simulation setup diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index e60912254a7..34bb4c9e791 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -17,8 +17,9 @@ function filter_eigvals(eig_vals, threshold) end end - println("$filtered_eigvals_counter eigenvalue(s) are not passed on because they either are in magnitude smaller than $threshold, \n - have positive real parts, or have negative imaginary parts") + println("$filtered_eigvals_counter 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") return length(filtered_eig_vals), filtered_eig_vals end From 6cb05722ab46a15fceabf0f8ddca25ce28aaa907 Mon Sep 17 00:00:00 2001 From: warisa-r <81345089+warisa-r@users.noreply.github.com> Date: Thu, 14 Mar 2024 11:22:58 +0100 Subject: [PATCH 37/64] Apply suggestions from code review Co-authored-by: Daniel Doehring --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 -- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- test/test_tree_1d_advection.jl | 4 ++-- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index f8c20adbdb9..86936209360 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -60,5 +60,3 @@ sol = Trixi.solve(ode, ode_algorithm, # Print the timer summary summary_callback() -# Show the error of the elixir -analysis_callback(sol) 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 6dd9957a9b2..d4a66af7fdd 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -280,7 +280,7 @@ function solve!(integrator::PERK2Integrator) end # Higher stages - for stage in 3:(alg.num_stages) + for stage in 3:alg.num_stages # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index d00bc82deeb..39ead7707e3 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -84,8 +84,8 @@ end @trixi_testset "elixir_advection_PERK2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_PERK2.jl"), - l2=[0.011662300515980219], - linf=[0.01647256923710194]) + l2=[0.0006911538643023453], + linf=[0.0009859411767757509]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 1f4e9d0a336d7857e4e9438abf067f35b2ba9085 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 14 Mar 2024 18:06:32 +0100 Subject: [PATCH 38/64] Update src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl Co-authored-by: Hendrik Ranocha --- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 34bb4c9e791..aa9c8b8f695 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -59,7 +59,6 @@ and partial differential equations. 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) """ - function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals) dtmin = 0.0 dt = -1.0 From 61da622e85ec9f73f5c34bf5b1a13ffeb7bfda39 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 09:47:10 +0100 Subject: [PATCH 39/64] Apply suggestions from code review --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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 d4a66af7fdd..141cdbbefcb 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -91,12 +91,12 @@ end PERK2() The following structures and methods provide a minimal implementation of -the second-order paired explicit Runge-Kutta method +the second-order paired explicit Runge-Kutta (PERK) method optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). -- 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) +- 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 PERK2 <: PERKSingle @@ -104,8 +104,8 @@ mutable struct PERK2 <: PERKSingle a_matrix::Matrix{Float64} c::Vector{Float64} - bS::Float64 b1::Float64 + bS::Float64 c_end::Float64 # Constructor that reads the coefficients from a file From 16ac4d55e25f9e8d72cd27a8b4633cfe36aea143 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 09:48:01 +0100 Subject: [PATCH 40/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 141cdbbefcb..ce8e38b0de5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -5,7 +5,8 @@ @muladd begin #! format: noindent -# Abstract base type for both single/standalone and multi-level P-ERK time integration schemes +# Abstract base type for both single/standalone and multi-level +P-ERK (Paired-Explicit Runge-Kutta) time integration schemes abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end From b3c12ee5dc8d76b30da6956c73f9922854f14901 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 09:48:19 +0100 Subject: [PATCH 41/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 ce8e38b0de5..96b22d9b45f 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -6,7 +6,7 @@ #! format: noindent # Abstract base type for both single/standalone and multi-level -P-ERK (Paired-Explicit Runge-Kutta) time integration schemes +# P-ERK (Paired-Explicit Runge-Kutta) time integration schemes abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end From cb0ef290c2baddca0df8edd5f40b65d45adfb030 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 09:48:39 +0100 Subject: [PATCH 42/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 96b22d9b45f..936056ec749 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -6,7 +6,7 @@ #! format: noindent # Abstract base type for both single/standalone and multi-level -# P-ERK (Paired-Explicit Runge-Kutta) time integration schemes +# PERK (Paired-Explicit Runge-Kutta) time integration schemes abstract type PERK end # Abstract base type for single/standalone P-ERK time integration schemes abstract type PERKSingle <: PERK end From a1506a0b6cb5584d278738ecf5be347d32ae5cf4 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 09:49:19 +0100 Subject: [PATCH 43/64] Apply suggestions from code review --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 86936209360..bbc446139eb 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -49,7 +49,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, # Define tspan to calculate maximum time step allowed for the bisection algorithm used in calculate polynomial coefficients in ODE algorithm tspan = (0.0, 1.0) -# Construct second order P-ERK method with 6 stages for +# Construct second order PERK method with 6 stages for # given simulation setup ode_algorithm = PERK2(6, tspan, semi) 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 936056ec749..df586b2284f 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -8,7 +8,7 @@ # Abstract base type for both single/standalone and multi-level # PERK (Paired-Explicit Runge-Kutta) time integration schemes abstract type PERK end -# Abstract base type for single/standalone P-ERK time integration schemes +# Abstract base type for single/standalone PERK time integration schemes abstract type PERKSingle <: PERK end function compute_a_coeffs(num_stage_evals, bc_factors, mon_coeffs) From 6fc11d585ddbf4525466aec5bb18f55cc9433bf6 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 15 Mar 2024 11:48:22 +0100 Subject: [PATCH 44/64] use readdlm instead of read_file --- src/auxiliary/auxiliary.jl | 14 -------------- .../paired_explicit_runge_kutta/methods_PERK2.jl | 9 +++++++-- 2 files changed, 7 insertions(+), 16 deletions(-) diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 3331a17cd25..92da9a5ba8b 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -346,20 +346,6 @@ function register_error_hints() return nothing end -function read_file(file_path::AbstractString, data_type::Type = Float64) - @assert isfile(file_path) "Couldn't find file" - data = zeros(data_type, 0) - open(file_path, "r") do file - while !eof(file) - line_content = readline(file) - append!(data, parse(data_type, line_content)) - end - end - num_lines = length(data) - - return num_lines, data -end - """ Trixi.download(src_url, file_path) 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 df586b2284f..26de1ecfbe0 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,6 +2,8 @@ # 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 + @muladd begin #! format: noindent @@ -78,7 +80,10 @@ function compute_PERK2_butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac a_matrix[:, 1] = c[3:end] path_mon_coeffs = base_path_mon_coeffs * "gamma_" * string(num_stages) * ".txt" - num_mon_coeffs, mon_coeffs = read_file(path_mon_coeffs, Float64) + @assert isfile(path_mon_coeffs) "Couldn't find file" + mon_coeffs = readdlm(path_mon_coeffs, Float64) + num_mon_coeffs = size(mon_coeffs, 1) + @assert num_mon_coeffs == coeffs_max A = compute_a_coeffs(num_stages, bc_factors, mon_coeffs) @@ -281,7 +286,7 @@ function solve!(integrator::PERK2Integrator) end # Higher stages - for stage in 3:alg.num_stages + for stage in 3:(alg.num_stages) # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + From ef1f5083bf149b4417f70e5dfe9a8fe7f99e338e Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 15 Mar 2024 11:57:33 +0100 Subject: [PATCH 45/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 26de1ecfbe0..f78e45cbe77 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,7 +2,7 @@ # 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 +using DelimitedFiles: readdlm @muladd begin #! format: noindent From 5bb554a4cbddd32bb446b5d6dc91e5160dcb6549 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 12:12:31 +0100 Subject: [PATCH 46/64] add DelimFIles --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 188e38b6093..960c7a0724c 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" 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" From 60cd20ecddafa4e683ffc824c1f15fa7a8d01880 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 12:13:19 +0100 Subject: [PATCH 47/64] fmt --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index bbc446139eb..cd1a005cdc4 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -59,4 +59,3 @@ sol = Trixi.solve(ode, ode_algorithm, # Print the timer summary summary_callback() - From 8e850ebaa1aae4389ac164f53257f3d635b6a3de Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 12:43:13 +0100 Subject: [PATCH 48/64] Unit tests --- test/test_unit.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/test/test_unit.jl b/test/test_unit.jl index 1907a281718..5f1035268f4 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -3,6 +3,8 @@ module TestUnit using Test using Trixi +using DelimitedFiles: readdlm + include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists @@ -1575,6 +1577,37 @@ end @test mesh.boundary_faces[:entire_boundary] == [1, 2] end + +@testset "PERK Single p2 Constructors" begin + path_coeff_file = joinpath(examples_dir(), "tree_1d_dgsem/") + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", + joinpath(path_coeff_file, "gamma_6.txt")) + + ode_algorithm = PERK2(6, path_coeff_file) + + @test ode_algorithm.a_matrix == [0.12405417889682908 0.07594582110317093 + 0.16178873711001726 0.13821126288998273 + 0.16692313960864164 0.2330768603913584 + 0.12281292901258256 0.37718707098741744] + + 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 = PERK2(12, tspan, vec(eig_vals)) + + @test ode_algorithm.a_matrix == [0.06460030228718522 0.026308788621905697 + 0.09476728053612687 0.04159635582750947 + 0.12338555995285924 0.058432621865322575 + 0.14992047611512016 0.0773522511576071 + 0.17346509036190638 0.09926218236536634 + 0.19265658773758804 0.12552523044423014 + 0.20526524445560954 0.1583711191807541 + 0.2073723615813174 0.20171854750959173 + 0.19137613735041836 0.26316931719503617 + 0.13943312463511776 0.36056687536488224] +end end end #module From 08faa0a2d726178e0a444768bb086b213018d65b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 13:46:32 +0100 Subject: [PATCH 49/64] compat --- Project.toml | 4 +++- test/Project.toml | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bae1b15b16a..7953eb1f8f8 100644 --- a/Project.toml +++ b/Project.toml @@ -23,7 +23,6 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" @@ -62,10 +61,13 @@ TrixiMakieExt = "Makie" [compat] CodeTracking = "1.0.5" ConstructionBase = "1.3" +Convex = "0.15.4" DataStructures = "0.18.15" +DelimitedFiles = "1.9" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" +ECOS = "1.1.12" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.24" diff --git a/test/Project.toml b/test/Project.toml index 1a042dab44f..d6821862adb 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -25,6 +26,7 @@ Plots = "1.19" Printf = "1" Random = "1" Test = "1" +DelimitedFiles = "1.9" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false From e17f7ad05b93482cfd31fc302234a25f9dc603d9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 13:48:59 +0100 Subject: [PATCH 50/64] ecos compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7953eb1f8f8..f0bf95eb883 100644 --- a/Project.toml +++ b/Project.toml @@ -67,7 +67,7 @@ DelimitedFiles = "1.9" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" -ECOS = "1.1.12" +ECOS = "1.1.2" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.24" From f3220e7fdd0640be3b882a859db7fc89b5495f2b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 13:55:22 +0100 Subject: [PATCH 51/64] compat --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f0bf95eb883..b9073b7fe21 100644 --- a/Project.toml +++ b/Project.toml @@ -63,7 +63,7 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" -DelimitedFiles = "1.9" +DelimitedFiles = "1.8.5" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index d6821862adb..3c16e5b4ab8 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,7 +26,7 @@ Plots = "1.19" Printf = "1" Random = "1" Test = "1" -DelimitedFiles = "1.9" +DelimitedFiles = "1.8.5" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false From 17096feb563b5b05aba96dde9f0415e3999d6daf Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 14:00:00 +0100 Subject: [PATCH 52/64] compat --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index b9073b7fe21..cca34e6a650 100644 --- a/Project.toml +++ b/Project.toml @@ -63,7 +63,7 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" -DelimitedFiles = "1.8.5" +DelimitedFiles = "1.9.0" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index 3c16e5b4ab8..91b34461493 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,7 +26,7 @@ Plots = "1.19" Printf = "1" Random = "1" Test = "1" -DelimitedFiles = "1.8.5" +DelimitedFiles = "1.9.0" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false From 01d04f2943403f58f76f33e64b65a14c617ee3f4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 14:03:45 +0100 Subject: [PATCH 53/64] compat --- test/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index 91b34461493..539d96d9eee 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,7 +26,6 @@ Plots = "1.19" Printf = "1" Random = "1" Test = "1" -DelimitedFiles = "1.9.0" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false From 44f413a2a5cbb54f3f0f319033c461e70cd0f3d5 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 14:05:58 +0100 Subject: [PATCH 54/64] compat --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 539d96d9eee..3878f69c4c6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" +DelimitedFiles = "1.9.0" Downloads = "1" FFMPEG = "0.4" ForwardDiff = "0.10.24" From 3280fc1ed0a790d5858bf57d3ddc4f17bf046647 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 15 Mar 2024 15:04:06 +0100 Subject: [PATCH 55/64] callbacks --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 f78e45cbe77..8c2a69405d4 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -315,7 +315,7 @@ function solve!(integrator::PERK2Integrator) # handle callbacks if callbacks isa CallbackSet - for cb in callbacks.discrete_callbacks + foreach(callbacks.discrete_callbacks) do cb if cb.condition(integrator.u, integrator.t, integrator) cb.affect!(integrator) end From 4a3deda2ce6abd712033cae9e29b96c0a6d963a0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 17 Mar 2024 19:35:34 +0100 Subject: [PATCH 56/64] increase allowed allocs --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 6 +++++- .../paired_explicit_runge_kutta/methods_PERK2.jl | 1 + test/test_tree_1d_advection.jl | 2 +- 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index cd1a005cdc4..7510e722b89 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -39,8 +39,12 @@ analysis_callback = AnalysisCallback(semi, interval = 100) # 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 = 1) + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, +callbacks = CallbackSet(summary_callback, + alive_callback, + analysis_callback, stepsize_callback) ############################################################################### 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 8c2a69405d4..3258485ee7f 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -319,6 +319,7 @@ function solve!(integrator::PERK2Integrator) if cb.condition(integrator.u, integrator.t, integrator) cb.affect!(integrator) end + return nothing end end diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 39ead7707e3..0869e8b3a81 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -92,7 +92,7 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 end end end From a1da64e7e5eea1a078703fd27ba66f2227c925c3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 09:50:08 +0100 Subject: [PATCH 57/64] fmt --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 7510e722b89..86c5d545062 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -42,7 +42,7 @@ stepsize_callback = StepsizeCallback(cfl = 2.5) alive_callback = AliveCallback(alive_interval = 1) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, stepsize_callback) From e804a63b6e962e8186a6a6cbb33fea5fdca16476 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 13:07:05 +0100 Subject: [PATCH 58/64] timer for step callbacks --- examples/tree_1d_dgsem/elixir_advection_PERK2.jl | 11 +++++------ .../paired_explicit_runge_kutta/methods_PERK2.jl | 14 ++++++++------ 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 86c5d545062..e27bc81960a 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -27,7 +27,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergen # ODE solvers, callbacks etc. # Create ODE problem with time span from 0.0 to 1.0 -ode = semidiscretize(semi, (0.0, 1.0)); +tspan = (0.0, 1.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 @@ -50,11 +51,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -# Define tspan to calculate maximum time step allowed for the bisection algorithm used in calculate polynomial coefficients in ODE algorithm -tspan = (0.0, 1.0) - -# Construct second order PERK method with 6 stages for -# given simulation setup +# Construct second order PERK 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 = PERK2(6, tspan, semi) sol = Trixi.solve(ode, ode_algorithm, 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 3258485ee7f..01e75d39da7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -313,13 +313,15 @@ function solve!(integrator::PERK2Integrator) integrator.iter += 1 integrator.t += integrator.dt - # handle callbacks - if callbacks isa CallbackSet - foreach(callbacks.discrete_callbacks) do cb - if cb.condition(integrator.u, integrator.t, integrator) - cb.affect!(integrator) + @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 - return nothing end end From 214fa6a9cbbdf636f7763393d7e32bf910a10e05 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 14:33:33 +0100 Subject: [PATCH 59/64] deps compat --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index cca34e6a650..b9073b7fe21 100644 --- a/Project.toml +++ b/Project.toml @@ -63,7 +63,7 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" -DelimitedFiles = "1.9.0" +DelimitedFiles = "1.8.5" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index 3878f69c4c6..515c27d3a3d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,7 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" -DelimitedFiles = "1.9.0" +DelimitedFiles = "1.8.5" Downloads = "1" FFMPEG = "0.4" ForwardDiff = "0.10.24" From 443a1625167445f6efe1a876c9219bc248076175 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 14:38:24 +0100 Subject: [PATCH 60/64] remove del files compat --- Project.toml | 1 - test/Project.toml | 1 - 2 files changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index b9073b7fe21..f4b164bb11e 100644 --- a/Project.toml +++ b/Project.toml @@ -63,7 +63,6 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" -DelimitedFiles = "1.8.5" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index 515c27d3a3d..539d96d9eee 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,7 +16,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" -DelimitedFiles = "1.8.5" Downloads = "1" FFMPEG = "0.4" ForwardDiff = "0.10.24" From 247946bc72420dd06d75519159b7afece0ef7644 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 16:07:38 +0100 Subject: [PATCH 61/64] skip delimitedfiles in downgrade compat --- .github/workflows/Downgrade.yml | 2 +- Project.toml | 1 + test/Project.toml | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index dd5d8ee7e32..974000eb430 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/Project.toml b/Project.toml index f4b164bb11e..f0bf95eb883 100644 --- a/Project.toml +++ b/Project.toml @@ -63,6 +63,7 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" +DelimitedFiles = "1.9" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index 539d96d9eee..2588a1fab43 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" +DelimitedFiles = "1.9" Downloads = "1" FFMPEG = "0.4" ForwardDiff = "0.10.24" From 1b3d4f64c5d2c1663c315eef1a984506822d67c5 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 18 Mar 2024 16:10:31 +0100 Subject: [PATCH 62/64] v1 --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f0bf95eb883..06b946fe8c9 100644 --- a/Project.toml +++ b/Project.toml @@ -63,7 +63,7 @@ CodeTracking = "1.0.5" ConstructionBase = "1.3" Convex = "0.15.4" DataStructures = "0.18.15" -DelimitedFiles = "1.9" +DelimitedFiles = "1" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" diff --git a/test/Project.toml b/test/Project.toml index 2588a1fab43..3222ae00ac5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,7 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" -DelimitedFiles = "1.9" +DelimitedFiles = "1" Downloads = "1" FFMPEG = "0.4" ForwardDiff = "0.10.24" From 1b94f8326e6c4c5e37e37573d59511f45d0acd62 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 20 Mar 2024 17:04:48 +0100 Subject: [PATCH 63/64] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 01e75d39da7..9d12ccd0c90 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -286,7 +286,7 @@ function solve!(integrator::PERK2Integrator) end # Higher stages - for stage in 3:(alg.num_stages) + for stage in 3:alg.num_stages # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + From 0d410e4a9ca4f5c610be82c1737c4cafd3386a98 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 12 Apr 2024 15:06:26 +0200 Subject: [PATCH 64/64] modular imp of the integrator --- .../methods_PERK2.jl | 152 ++++++++++-------- 1 file changed, 85 insertions(+), 67 deletions(-) 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 9d12ccd0c90..cfa1d12086e 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -200,6 +200,7 @@ mutable struct PERK2Integrator{RealT <: Real, uType, Params, Sol, F, Alg, # PERK2 stages: k1::uType k_higher::uType + t_stage::RealT end # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) @@ -211,9 +212,8 @@ function Base.getproperty(integrator::PERKIntegrator, field::Symbol) return getfield(integrator, field) end -# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 -function solve(ode::ODEProblem, alg::PERK2; - dt, callback = nothing, kwargs...) +function init(ode::ODEProblem, alg::PERK2; + dt, callback = nothing, kwargs...) u0 = copy(ode.u0) du = zero(u0) u_tmp = zero(u0) @@ -229,7 +229,7 @@ function solve(ode::ODEProblem, alg::PERK2; (prob = ode,), ode.f, alg, PERKIntegratorOptions(callback, ode.tspan; kwargs...), false, - k1, k_higher) + k1, k_higher, t0) # initialize callbacks if callback isa CallbackSet @@ -243,10 +243,33 @@ function solve(ode::ODEProblem, alg::PERK2; error("unsupported") end - solve!(integrator) + return integrator +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::PERK2; + dt, callback = nothing, kwargs...) + integrator = init(ode, alg, dt = dt, callback = callback; kwargs...) + + # Start actual solve + solve_steps!(integrator) end -function solve!(integrator::PERK2Integrator) +function solve_steps!(integrator::PERK2Integrator) + @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::PERK2Integrator) @unpack prob = integrator.sol @unpack alg = integrator t_end = last(prob.tspan) @@ -254,87 +277,82 @@ function solve!(integrator::PERK2Integrator) integrator.finalstep = false - @trixi_timeit timer() "main loop" while !integrator.finalstep - if isnan(integrator.dt) - error("time step size `dt` is NaN") + @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 - # 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) + # 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) - @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 + @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.du) - integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + 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 - # k2 + integrator.f(integrator.du, integrator.u_tmp, prob.p, - integrator.t + alg.c[2] * integrator.dt) + 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 - # Higher stages - for stage in 3:alg.num_stages - # 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.u) + integrator.u[i] += alg.b1 * integrator.k1[i] + + alg.bS * integrator.k_higher[i] + end + end # PERK2 step - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt - end - end + integrator.iter += 1 + integrator.t += integrator.dt - @threaded for i in eachindex(integrator.u) - integrator.u[i] += alg.b1 * integrator.k1[i] + - alg.bS * integrator.k_higher[i] - end - end # PERK2 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 + @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 # "main loop" timer - - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) + # 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