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 27df49ed4fa..04c9f2f3d21 100644 --- a/Project.toml +++ b/Project.toml @@ -6,10 +6,13 @@ version = "0.7.5-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" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -57,10 +60,13 @@ TrixiMakieExt = "Makie" [compat] CodeTracking = "1.0.5" ConstructionBase = "1.3" +Convex = "0.15.4" DataStructures = "0.18.15" +DelimitedFiles = "1" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" Downloads = "1.6" +ECOS = "1.1.2" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.24" 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..e27bc81960a --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -0,0 +1,64 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = 1.0 +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = -1.0 # minimum coordinate +coordinates_max = 1.0 # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 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 +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) + +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, + alive_callback, + analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# 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, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 9375c80d77e..00aa9583081 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -276,6 +276,8 @@ export trixi_include, examples_dir, get_examples, default_example, export ode_norm, ode_unstable_check +export PERK2 + export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP 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..cfa1d12086e --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -0,0 +1,389 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +using DelimitedFiles: readdlm + +@muladd begin +#! format: noindent + +# 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 PERK time integration schemes +abstract type PERKSingle <: PERK end + +function compute_a_coeffs(num_stage_evals, bc_factors, mon_coeffs) + a_coeffs = copy(mon_coeffs) + + for stage in 1:(num_stage_evals - 2) + a_coeffs[stage] /= bc_factors[stage] + for prev_stage in 1:(stage - 1) + a_coeffs[stage] /= a_coeffs[prev_stage] + end + end + + return reverse(a_coeffs) +end + +function compute_PERK2_butcher_tableau(num_stages, eig_vals, tspan, + bS, c_end) + + # 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 + 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 + + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] + + cons_order = 2 + filter_thres = 1e-12 + dtmax = tspan[2] - tspan[1] + dteps = 1e-9 + + 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, + 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, bc_factors, mon_coeffs) + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + return a_matrix, 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) + for k in 2:num_stages + c[k] = c_end * (k - 1) / (num_stages - 1) + end + 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 + + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] + + path_mon_coeffs = base_path_mon_coeffs * "gamma_" * string(num_stages) * ".txt" + @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) + + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + + return a_matrix, c +end + +""" + PERK2() + +The following structures and methods provide a minimal implementation of +the second-order paired explicit Runge-Kutta (PERK) method +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). + +- Brian Vermeire (2019). + Paired explicit Runge-Kutta schemes for stiff systems of equations + [DOI: 10.1016/j.jcp.2019.05.014](https://doi.org/10.1016/j.jcp.2019.05.014) +""" + +mutable struct PERK2 <: PERKSingle + const num_stages::Int + + a_matrix::Matrix{Float64} + c::Vector{Float64} + b1::Float64 + bS::Float64 + c_end::Float64 + + # 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) + + newPERK2.a_matrix, newPERK2.c = compute_PERK2_butcher_tableau(num_stages, + base_path_mon_coeffs, + bS, c_end) + + newPERK2.b1 = one(bS) - bS + newPERK2.bS = bS + newPERK2.c_end = c_end + return newPERK2 + end + + # 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)) + 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 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) + + 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 +end # struct PERK2 + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 +mutable struct PERKIntegratorOptions{Callback} + callback::Callback # callbacks; used in Trixi + adaptive::Bool # whether the algorithm is adaptive; ignored + dtmax::Float64 # ignored + maxiters::Int # maximal number of time steps + tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored +end + +function PERKIntegratorOptions(callback, tspan; maxiters = typemax(Int), kwargs...) + PERKIntegratorOptions{typeof(callback)}(callback, false, Inf, maxiters, + [last(tspan)]) +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 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::PERKIntegrator, field::Symbol) + if field === :stats + return (naccept = getfield(integrator, :iter),) + end + # general fallback + return getfield(integrator, field) +end + +function init(ode::ODEProblem, alg::PERK2; + dt, callback = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(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 + + 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_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) + callbacks = integrator.opts.callback + + integrator.finalstep = false + + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end + + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end + + @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin + # k1 + integrator.f(integrator.du, integrator.u, prob.p, integrator.t) + @threaded for i in eachindex(integrator.du) + integrator.k1[i] = integrator.du[i] * integrator.dt + end + + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] + end + # k2 + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[2] * integrator.dt) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + # Higher stages + for stage in 3:(alg.num_stages) + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.a_matrix[stage - 2, 1] * + integrator.k1[i] + + alg.a_matrix[stage - 2, 2] * + integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[stage] * integrator.dt) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + end + + @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 + end + end + end + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end +end + +# get a cache where the RHS can be stored +get_du(integrator::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::PERKIntegrator, ::Bool) = false + +# used by adaptive timestepping algorithms in DiffEq +function set_proposed_dt!(integrator::PERKIntegrator, dt) + integrator.dt = dt +end + +function get_proposed_dt(integrator::PERKIntegrator) + return integrator.dt +end + +# stop the time integration +function terminate!(integrator::PERKIntegrator) + integrator.finalstep = true + empty!(integrator.opts.tstops) +end + +# used for AMR (Adaptive Mesh Refinement) +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) +end +end # @muladd diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl new file mode 100644 index 00000000000..df9b1b6c7cf --- /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("polynomial_optimizer.jl") + +include("methods_PERK2.jl") +end # @muladd diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl new file mode 100644 index 00000000000..aa9c8b8f695 --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -0,0 +1,137 @@ +using LinearAlgebra: eigvals +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 + # 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 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 + +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 + +""" + bisection() + +The following structures and methods provide a simplified implementation to +discover optimal stability polynomial for a given set of `eig_vals` +These are designed for the one-step (i.e., Runge-Kutta methods) integration of initial value ordinary +and partial differential equations. + +- Ketcheson and Ahmadia (2012). +Optimal stability polynomials for numerical integration of initial value problems +[DOI: 10.2140/camcos.2012.7.247](https://doi.org/10.2140/camcos.2012.7.247) +""" +function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals) + dtmin = 0.0 + dt = -1.0 + abs_p = -1.0 + + # Construct stability polynomial for each eigenvalue + pnoms = ones(Complex{Float64}, num_eig_vals, 1) + + # Init datastructure for results + gamma = Variable(num_stage_evals - cons_order) + + normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) + + for j in 1:num_stage_evals + fac_j = factorial(j) + for i in 1:num_eig_vals + normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j + end + end + + normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, + num_stage_evals) + + println("Start optimization of stability polynomial \n") + + while dtmax - dtmin > dteps + dt = 0.5 * (dtmax + dtmin) + + # Compute stability polynomial for current timestep + for k in 1:num_stage_evals + dt_k = dt^k + for i in 1:num_eig_vals + normalized_powered_eigvals_scaled[i, k] = dt_k * + normalized_powered_eigvals[i, k] + end + end + + # Use last optimal values for gamma in (potentially) next iteration + problem = minimize(stability_polynomials(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 + + if abs_p < 1.0 + dtmin = dt + else + dtmax = dt + end + end + println("Concluded stability polynomial optimization \n") + + return evaluate(gamma), dt +end + +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 + return gamma_opt +end diff --git a/src/time_integration/time_integration.jl b/src/time_integration/time_integration.jl index c1e53527121..d19a1fcc37c 100644 --- a/src/time_integration/time_integration.jl +++ b/src/time_integration/time_integration.jl @@ -16,4 +16,5 @@ end include("methods_2N.jl") include("methods_3Sstar.jl") include("methods_SSP.jl") +include("paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl") end # @muladd diff --git a/test/Project.toml b/test/Project.toml index 1491d7a5c5f..5b795891391 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" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" @@ -16,6 +17,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CairoMakie = "0.10" +DelimitedFiles = "1" Downloads = "1" ExplicitImports = "1.0.1" FFMPEG = "0.4" diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index a580a3b5600..0869e8b3a81 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.0006911538643023453], + linf=[0.0009859411767757509]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 03a78f6918a..66904d6c55e 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 @@ -1576,6 +1578,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