From 4f958718c40ca3a9d7c38c2acfc0ae0df7aa0fb0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 5 Nov 2024 17:27:02 +0100 Subject: [PATCH 01/61] PERK4 with imported coefficients --- .../elixir_euler_vortex_perk4.jl | 115 ++++++++ .../methods_PERK2.jl | 27 +- .../methods_PERK3.jl | 16 +- .../methods_PERK4.jl | 262 ++++++++++++++++++ .../paired_explicit_runge_kutta.jl | 5 +- 5 files changed, 402 insertions(+), 23 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl create mode 100644 src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl new file mode 100644 index 00000000000..0f918bed16c --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -0,0 +1,115 @@ + +using OrdinaryDiffEq # Requiref for `CallbackSet` +using Trixi + +# Ratio of specific heats +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc) + +""" + initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D) + +The classical isentropic vortex test case as presented in Section 5.1 of + +- 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) + https://spectrum.library.concordia.ca/id/eprint/985444/1/Paired-explicit-Runge-Kutta-schemes-for-stiff-sy_2019_Journal-of-Computation.pdf +""" +function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D) + # Evaluate error after full domain traversion + if t == T_end + t = 0 + end + + # initial center of the vortex + inicenter = SVector(0.0, 0.0) + # strength of the vortex + S = 13.5 + # Radius of vortex + R = 1.5 + # Free-stream Mach + M = 0.4 + # base flow + v1 = 1.0 + v2 = 1.0 + vel = SVector(v1, v2) + + center = inicenter + vel * t # advection of center + center = x - center # distance to centerpoint + center = SVector(center[2], -center[1]) + r2 = center[1]^2 + center[2]^2 + + f = (1 - r2) / (2 * R^2) + + rho = (1 - (S * M / pi)^2 * (gamma - 1) * exp(2 * f) / 8)^(1 / (gamma - 1)) + + du = S / (2 * π * R) * exp(f) # vel. perturbation + vel = vel + du * center + v1, v2 = vel + + p = rho^gamma / (gamma * M^2) + prim = SVector(rho, v1, v2, p) + return prim2cons(prim, equations) +end +initial_condition = initial_condition_isentropic_vortex + +EdgeLength = 20 + +N_passes = 1 +T_end = EdgeLength * N_passes +tspan = (0.0, T_end) + +coordinates_min = (-EdgeLength / 2, -EdgeLength / 2) +coordinates_max = (EdgeLength / 2, EdgeLength / 2) + +cells_per_dimension = (32, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_integrals = (entropy,)) + +# Note quite large CFL number +stepsize_callback = StepsizeCallback(cfl = 9.1) + +callbacks = CallbackSet(summary_callback, + stepsize_callback, + analysis_callback) + +############################################################################### +# set up time integration algorithm + +num_stages = 19 + +current_directory = @__DIR__ +coefficient_file = "a_" * string(num_stages) * ".txt" + +# Download the optimized PERK4 coefficients +Trixi.download("https://gist.githubusercontent.com/DanielDoehring/84f266ff61f0a69a0127cec64056275e/raw/1a66adbe1b425d33daf502311ecbdd4b191b89cc/a_19.txt", + joinpath(current_directory, coefficient_file)) + +ode_algorithm = Trixi.PairedExplicitRK4(num_stages, current_directory) + +############################################################################### +# run the simulation + +sol = Trixi.solve(ode, ode_algorithm, + dt = 42.0, # Not used + save_everystep = false, callback = callbacks); + +# Clean up the downloaded file +rm(joinpath(current_directory, coefficient_file)) + +summary_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 0ff648556c4..d196c70db1f 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -128,7 +128,8 @@ optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solve 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) -Note: To use this integrator, the user must import the `Convex` and `ECOS` packages. +Note: To use this integrator, the user must import the `Convex` and `ECOS` packages +unless the coefficients are provided in a "gamma_.txt" file. """ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle const num_stages::Int @@ -240,6 +241,20 @@ function init(ode::ODEProblem, alg::PairedExplicitRK2; return integrator end +# Function that computes the first stage of a general P-ERK method +# (in fact every explicit Runge-Kutta method) +@inline function k1!(integrator::AbstractPairedExplicitRKIntegrator, p, c) + integrator.f(integrator.du, integrator.u, 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.u) + integrator.u_tmp[i] = integrator.u[i] + c[2] * integrator.k1[i] + end +end + function step!(integrator::PairedExplicitRK2Integrator) @unpack prob = integrator.sol @unpack alg = integrator @@ -261,16 +276,8 @@ function step!(integrator::PairedExplicitRK2Integrator) 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 + k1!(integrator, prob.p, alg.c) - # Construct current state - @threaded for i in eachindex(integrator.u) - integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i] - end # k2 integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t + alg.c[2] * integrator.dt) 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 384f4b408ca..5b1396d8204 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -132,7 +132,8 @@ While the changes to SSPRK33 base-scheme are described in Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively Refined Meshes for Compressible Fluid Dynamics [DOI: 10.1016/j.jcp.2024.113223](https://doi.org/10.1016/j.jcp.2024.113223) -Note: To use this integrator, the user must import the `Convex`, `ECOS`, and `NLsolve` packages. +Note: To use this integrator, the user must import the `Convex`, `ECOS`, and `NLsolve` packages +unless the A-matrix coefficients are provided in a "a_.txt" file. """ mutable struct PairedExplicitRK3 <: AbstractPairedExplicitRKSingle const num_stages::Int # S @@ -225,8 +226,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK3; # initialize callbacks if callback isa CallbackSet for cb in callback.continuous_callbacks - throw(ArgumentError("Continuous callbacks are unsupported with paired explicit Runge- - Kutta methods.")) + throw(ArgumentError("Continuous callbacks are unsupported with paired explicit Runge-Kutta methods.")) end for cb in callback.discrete_callbacks cb.initialize(cb, integrator.u, integrator.t, integrator) @@ -257,16 +257,8 @@ function step!(integrator::PairedExplicitRK3Integrator) 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 + k1!(integrator, prob.p, alg.c) - # 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) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl new file mode 100644 index 00000000000..358030ac84b --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -0,0 +1,262 @@ +# 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 + +function compute_PairedExplicitRK4_butcher_tableau(num_stages, + base_path_a_coeffs::AbstractString; + c_const = 1.0) # Default value for best internal stability + c = c_const * ones(num_stages) # Use same abscissae for free coefficients + c[1] = 0.0 + + cS3 = c_const + c[num_stages - 3] = cS3 + c[num_stages - 2] = 0.479274057836310 + c[num_stages - 1] = sqrt(3) / 6 + 0.5 + c[num_stages] = -sqrt(3) / 6 + 0.5 + + num_coeffs_max = num_stages - 5 + + a_matrix = zeros(num_coeffs_max, 2) + a_matrix[:, 1] = c[3:(num_stages - 3)] + + path_a_coeffs = joinpath(base_path_a_coeffs, + "a_" * string(num_stages) * ".txt") + + @assert isfile(path_a_coeffs) "Couldn't find file $path_a_coeffs" + a_coeffs = readdlm(path_a_coeffs, Float64) + num_a_coeffs = size(a_coeffs, 1) + + @assert num_a_coeffs == num_coeffs_max + if num_coeffs_max > 0 + a_matrix[:, 1] -= a_coeffs + a_matrix[:, 2] = a_coeffs + end + + # Constant/non-optimized part of the Butcher matrix + a_matrix_constant = [0.479274057836310-0.114851811257441 / cS3 0.114851811257441/cS3 + 0.1397682537005989 0.648906880894214 + 0.1830127018922191 0.028312163512968] + + return a_matrix, a_matrix_constant, c +end + +@doc raw""" + PairedExplicitRK4(num_stages, base_path_a_coeffs::AbstractString, dt_opt = nothing; + c_const = 1.0f0) + + Parameters: + - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method. + - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in + the Butcher tableau of the Runge Kutta method. + The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks. + - `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. + In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used. + - `c_const` (`Float64`, optional): Value of abscissae $c$ in the Butcher tableau for the optimized coefficients. + Default is 1.0. + +The following structures and methods provide an implementation of +the fourth-order paired explicit Runge-Kutta (P-ERK) method +optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). +The method has been proposed in +- D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024). + Fourth-Order Paired-Explicit Runge-Kutta Methods + [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470) +""" +mutable struct PairedExplicitRK4 <: AbstractPairedExplicitRKSingle + const num_stages::Int # S + + a_matrix::Matrix{Float64} + # This part of the Butcher array matrix A is constant for all PERK methods, i.e., + # regardless of the optimized coefficients. + a_matrix_constant::Matrix{Float64} + c::Vector{Float64} + dt_opt::Union{Float64, Nothing} +end # struct PairedExplicitRK4 + +# Constructor for previously computed A Coeffs +function PairedExplicitRK4(num_stages, base_path_a_coeffs::AbstractString, + dt_opt = nothing; + c_const = 1.0f0) + a_matrix, a_matrix_constant, c = compute_PairedExplicitRK4_butcher_tableau(num_stages, + base_path_a_coeffs; + c_const) + + return PairedExplicitRK4(num_stages, a_matrix, a_matrix_constant, c, dt_opt) +end + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77 +# This implements the interface components described at +# https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1 +# which are used in Trixi.jl. +mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F, Alg, + PairedExplicitRKOptions} <: + AbstractPairedExplicitRKSingleIntegrator + u::uType + du::uType + u_tmp::uType + t::RealT + tdir::RealT + dt::RealT # current time step + dtcache::RealT # manually set time step + iter::Int # current number of time steps (iteration) + p::Params # will be the semidiscretization from Trixi + sol::Sol # faked + f::F + alg::Alg # This is our own class written above; Abbreviation for ALGorithm + opts::PairedExplicitRKOptions + finalstep::Bool # added for convenience + dtchangeable::Bool + force_stepfail::Bool + # PairedExplicitRK stages: + k1::uType + k_higher::uType +end + +function init(ode::ODEProblem, alg::PairedExplicitRK4; + dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...) + u0 = copy(ode.u0) + du = zero(u0) + u_tmp = zero(u0) + + # PairedExplicitRK stages + k1 = zero(u0) + k_higher = zero(u0) + + t0 = first(ode.tspan) + tdir = sign(ode.tspan[end] - ode.tspan[1]) + iter = 0 + + integrator = PairedExplicitRK4Integrator(u0, du, u_tmp, t0, tdir, dt, dt, iter, + ode.p, + (prob = ode,), ode.f, alg, + PairedExplicitRKOptions(callback, + ode.tspan; + kwargs...), + false, true, false, + k1, k_higher) + + # initialize callbacks + if callback isa CallbackSet + for cb in callback.continuous_callbacks + throw(ArgumentError("Continuous callbacks are unsupported with paired explicit Runge-Kutta methods.")) + end + for cb in callback.discrete_callbacks + cb.initialize(cb, integrator.u, integrator.t, integrator) + end + end + + return integrator +end + +@inline function last_three_stages!(integrator::PairedExplicitRK4Integrator, alg, p) + for stage in 1:2 + @threaded for u_ind in eachindex(integrator.u) + integrator.u_tmp[u_ind] = integrator.u[u_ind] + + alg.a_matrix_constant[stage, 1] * + integrator.k1[u_ind] + + alg.a_matrix_constant[stage, 2] * + integrator.k_higher[u_ind] + end + + integrator.f(integrator.du, integrator.u_tmp, p, + integrator.t + + alg.c[alg.num_stages - 3 + stage] * integrator.dt) + + @threaded for u_ind in eachindex(integrator.du) + integrator.k_higher[u_ind] = integrator.du[u_ind] * integrator.dt + end + end + + # Last stage + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.a_matrix_constant[3, 1] * integrator.k1[i] + + alg.a_matrix_constant[3, 2] * integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, p, + integrator.t + alg.c[alg.num_stages] * integrator.dt) + + @threaded for u_ind in eachindex(integrator.u) + # Note that 'k_higher' carries the values of K_{S-1} + # and that we construct 'K_S' "in-place" from 'integrator.du' + integrator.u[u_ind] += 0.5 * (integrator.k_higher[u_ind] + + integrator.du[u_ind] * integrator.dt) + end +end + +function step!(integrator::PairedExplicitRK4Integrator) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end + + modify_dt_for_tstops!(integrator) + + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end + + @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin + k1!(integrator, prob.p, alg.c) + + # 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 until "constant" stages + for stage in 3:(alg.num_stages - 3) + # Construct current state + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.a_matrix[stage - 2, 1] * + integrator.k1[i] + + alg.a_matrix[stage - 2, 2] * + integrator.k_higher[i] + end + + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + alg.c[stage] * integrator.dt) + + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + end + + last_three_stages!(integrator, alg, prob.p) + end # PairedExplicitRK step timer + + integrator.iter += 1 + integrator.t += integrator.dt + + # handle callbacks + if callbacks isa CallbackSet + for cb in callbacks.discrete_callbacks + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) + end + end + end + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end +end +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 c606326738f..7183f0ff823 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 @@ -9,7 +9,7 @@ include("polynomial_optimizer.jl") # Abstract base type for both single/standalone and multi-level -# PERK (Paired-Explicit Runge-Kutta) time integration schemes +# PERK (Paired Explicit Runge-Kutta) time integration schemes abstract type AbstractPairedExplicitRK end # Abstract base type for single/standalone PERK time integration schemes abstract type AbstractPairedExplicitRKSingle <: AbstractPairedExplicitRK end @@ -178,5 +178,8 @@ function solve_a_butcher_coeffs_unknown! end # Basic implementation of the second-order paired explicit Runge-Kutta (PERK) method include("methods_PERK2.jl") +# Slightly customized implementation of the third-order PERK method include("methods_PERK3.jl") +# Basic implementation of the fourth-order PERK method +include("methods_PERK4.jl") end # @muladd From 3b84504cf436079313a5ebb21b4089b8e70a04d5 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Nov 2024 14:46:51 +0100 Subject: [PATCH 02/61] examples, more constructors --- .../elixir_hypdiff_nonperiodic_perk4.jl | 65 +++++++++ ext/TrixiConvexECOSExt.jl | 138 +++++++++++++++++- .../methods_PERK2.jl | 15 +- .../methods_PERK3.jl | 11 +- .../methods_PERK4.jl | 107 ++++++++++++-- .../polynomial_optimizer.jl | 2 + test/test_structured_2d.jl | 28 ++++ test/test_tree_1d_hypdiff.jl | 19 +++ test/test_unit.jl | 21 +++ 9 files changed, 379 insertions(+), 27 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl new file mode 100644 index 00000000000..676cc2a9da6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl @@ -0,0 +1,65 @@ + +using OrdinaryDiffEq +using Trixi + +# Convex and ECOS are imported because they are used for finding the optimal time step and optimal +# monomial coefficients in the stability polynomial of P-ERK time integrators. +using Convex, ECOS + +############################################################################### +# semidiscretization of the hyperbolic diffusion equations + +equations = HyperbolicDiffusionEquations1D() + +initial_condition = initial_condition_poisson_nonperiodic + +boundary_conditions = boundary_condition_poisson_nonperiodic + +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) + +coordinates_min = 0.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 30_000, + periodicity = false) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions, + source_terms = source_terms_poisson_nonperiodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +resid_tol = 5.0e-12 +steady_state_callback = SteadyStateCallback(abstol = resid_tol, reltol = 0.0) + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +# Construct third order paired explicit Runge-Kutta method with 8 stages for given simulation setup. +# Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used +# in calculating the polynomial coefficients in the ODE algorithm. +ode_algorithm = Trixi.PairedExplicitRK4(11, tspan, semi) + +cfl_number = Trixi.calculate_cfl(ode_algorithm, ode) +stepsize_callback = StepsizeCallback(cfl = 0.9 * cfl_number) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = Trixi.solve(ode, ode_algorithm, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 8251fe3eed9..c992dc55878 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -15,7 +15,8 @@ end using LinearAlgebra: eigvals # Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd +using Trixi: Trixi, undo_normalization!, undo_normalization_PERK4!, + bisect_stability_polynomial, bisect_stability_polynomial_PERK4, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). # Since these FMAs can increase the performance of many numerical algorithms, @@ -28,10 +29,14 @@ using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd # relative to consistency order. function Trixi.undo_normalization!(gamma_opt, consistency_order, num_stage_evals) for k in (consistency_order + 1):num_stage_evals - gamma_opt[k - consistency_order] = gamma_opt[k - consistency_order] / - factorial(k) + gamma_opt[k - consistency_order] /= factorial(k) + end +end + +function Trixi.undo_normalization_PERK4!(gamma_opt, num_stage_evals) + for k in 1:(num_stage_evals - 5) + gamma_opt[k] /= factorial(k + 4) end - return gamma_opt end # Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency @@ -63,6 +68,44 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals, return maximum(abs(pnoms)) end +# Specialized form of the stability polynomials for fourth-order PERK schemes. +function stability_polynomials_PERK4!(pnoms, num_stage_evals, + normalized_powered_eigvals, + gamma, dt, cS3) + num_eig_vals = length(pnoms) + + # Constants arising from the particular form of Butcher tableau chosen for the 4th order PERK methods + k1 = 0.001055026310046423 / cS3 + k2 = 0.03726406530405851 / cS3 + # Note: `cS3` = c_{S-3} is in principle free, while the other abscissae are fixed to 1.0 + + # Initialize with zero'th order (z^0) coefficient + for i in 1:num_eig_vals + pnoms[i] = 1.0 + end + + # First `consistency_order` = 4 terms of the exponential + for k in 1:4 + for i in 1:num_eig_vals + pnoms[i] += dt^k * normalized_powered_eigvals[i, k] + end + end + + # "Fixed" term due to choice of the PERK4 Butcher tableau + # Required to un-do the normalization of the eigenvalues here + pnoms += k1 * dt^5 * normalized_powered_eigvals[:, 5] * factorial(5) + + # Contribution from free coefficients + for k in 1:(num_stage_evals - 5) + pnoms += (k2 * dt^(k + 4) * normalized_powered_eigvals[:, k + 4] * gamma[k] + + k1 * dt^(k + 5) * normalized_powered_eigvals[:, k + 5] * gamma[k] * + (k + 5)) + end + + # For optimization only the maximum is relevant + return maximum(abs(pnoms)) +end + #= The following structures and methods provide a simplified implementation to discover optimal stability polynomial for a given set of `eig_vals` @@ -158,6 +201,93 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, gamma_opt = [gamma_opt] end + undo_normalization!(gamma_opt, consistency_order, num_stage_evals) + + return gamma_opt, dt +end + +# Specialized routine for PERK4. +# For details, see Section 4 in +# - D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024). +# Fourth-Order Paired-Explicit Runge-Kutta Methods +# [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470) +function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, + num_stage_evals, + dtmax, dteps, eig_vals, cS3; + verbose = false) + dtmin = 0.0 + dt = -1.0 + abs_p = -1.0 + + # Construct stability polynomial for each eigenvalue + pnoms = ones(Complex{Float64}, num_eig_vals, 1) + + # Init datastructure for monomial coefficients + gamma = Variable(num_stage_evals - 5) + + 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 + + if verbose + println("Start optimization of stability polynomial \n") + end + + # Bisection on timestep + while dtmax - dtmin > dteps + dt = 0.5 * (dtmax + dtmin) + + # Use last optimal values for gamma in (potentially) next iteration + problem = minimize(stability_polynomials_PERK4!(pnoms, + num_stage_evals, + normalized_powered_eigvals, + gamma, dt, cS3)) + + solve!(problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + "feastol" => 1e-9, + "abstol" => 1e-9, + "reltol" => 1e-9, + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent = true) + + abs_p = problem.optval + + if abs_p < 1 + dtmin = dt + else + dtmax = dt + end + end + + if verbose + println("Concluded stability polynomial optimization \n") + end + + gamma_opt = [] + + if num_stage_evals > 5 + gamma_opt = evaluate(gamma) + + # Catch case S = 6 (only one opt. variable) + if isa(gamma_opt, Number) + gamma_opt = [gamma_opt] + end + + undo_normalization_PERK4!(gamma_opt, num_stage_evals) + end + return gamma_opt, dt 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 d196c70db1f..8ce2bb7b4ed 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -49,8 +49,6 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, dtmax, dteps, eig_vals; verbose) - monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, - num_stages) num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max @@ -112,13 +110,13 @@ end In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used. - `tspan`: Time span of the simulation. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the + - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the equation has been semidiscretized. - `verbose` (`Bool`, optional): Verbosity flag, default is false. - - `bS` (`Float64`, optional): Value of b in the Butcher tableau at b_s, when - s is the number of stages, default is 1.0. - - `cS` (`Float64`, optional): Value of c in the Butcher tableau at c_s, when - s is the number of stages, default is 0.5. + - `bS` (`Float64`, optional): Value of $b_S$ in the Butcher tableau, where + $S$ is the number of stages. Default is 1.0. + - `cS` (`Float64`, optional): Value of $c_S$ in the Butcher tableau, where + $S$ is the number of stages. Default is 0.5. The following structures and methods provide a minimal implementation of the second-order paired explicit Runge-Kutta (PERK) method @@ -146,6 +144,7 @@ end # struct PairedExplicitRK2 function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing, bS = 1.0, cS = 0.5) + @assert num_stages>=2 "PERK2 requires at least two stages" # If the user has the monomial coefficients, they also must have the optimal time step a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs, @@ -159,6 +158,7 @@ end function PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) + @assert num_stages>=2 "PERK2 requires at least two stages" eig_vals = eigvals(jacobian_ad_forward(semi)) return PairedExplicitRK2(num_stages, tspan, eig_vals; verbose, bS, cS) @@ -169,6 +169,7 @@ end function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, bS = 1.0, cS = 0.5) + @assert num_stages>=2 "PERK2 requires at least two stages" a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, bS, cS; 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 5b1396d8204..c6dd328b10e 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -48,8 +48,6 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, num_eig_vals, num_stages, dtmax, dteps, eig_vals; verbose) - monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, - num_stages) # Solve the nonlinear system of equations from monomial coefficient and # Butcher array abscissae c to find Butcher matrix A @@ -114,11 +112,11 @@ end In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used. - `tspan`: Time span of the simulation. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the + - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the equation has been semidiscretized. - `verbose` (`Bool`, optional): Verbosity flag, default is false. - - `cS2` (`Float64`, optional): Value of c in the Butcher tableau at c_{s-2}, when - s is the number of stages, default is 1.0f0. + - `cS2` (`Float64`, optional): Value of $c_{S-2}$ in the Butcher tableau, where + $S$ is the number of stages. Default is 1.0f0. The following structures and methods provide an implementation of the third-order paired explicit Runge-Kutta (P-ERK) method @@ -147,6 +145,7 @@ end # struct PairedExplicitRK3 function PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt = nothing; cS2 = 1.0f0) + @assert num_stages>=3 "PERK3 requires at least three stages" a_matrix, c = compute_PairedExplicitRK3_butcher_tableau(num_stages, base_path_a_coeffs; cS2) @@ -157,6 +156,7 @@ end # Constructor that computes Butcher matrix A coefficients from a semidiscretization function PairedExplicitRK3(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, cS2 = 1.0f0) + @assert num_stages>=3 "PERK3 requires at least three stages" eig_vals = eigvals(jacobian_ad_forward(semi)) return PairedExplicitRK3(num_stages, tspan, eig_vals; verbose, cS2) @@ -165,6 +165,7 @@ end # Constructor that calculates the coefficients with polynomial optimizer from a list of eigenvalues function PairedExplicitRK3(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, cS2 = 1.0f0) + @assert num_stages>=3 "PERK3 requires at least three stages" a_matrix, c, dt_opt = compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, eig_vals; diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 358030ac84b..dcce0409924 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -5,13 +5,70 @@ @muladd begin #! format: noindent +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 4 +# using a list of eigenvalues +function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, + eig_vals::Vector{ComplexF64}; + verbose = false, cS3) + c = ones(num_stages) # Best internal stability properties + c[1] = 0.0 + + c[num_stages - 3] = cS3 + c[num_stages - 2] = 0.479274057836310 + c[num_stages - 1] = sqrt(3) / 6 + 0.5 + c[num_stages] = -sqrt(3) / 6 + 0.5 + + num_coeffs_max = num_stages - 5 + + a_matrix = nothing + dt_opt = nothing + if num_stages > 5 + # Calculate coefficients of the stability polynomial in monomial form + dtmax = tspan[2] - tspan[1] + dteps = 1.0f-9 + + num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + + monomial_coeffs, dt_opt = bisect_stability_polynomial_PERK4(num_eig_vals, + num_stages, + dtmax, dteps, + eig_vals, cS3; + verbose) + + a_unknown = zeros(num_coeffs_max) + a_unknown[1] = monomial_coeffs[1] + l = 2 + for _ in 5:(num_stages - 2) + a_unknown[l] = monomial_coeffs[l] / monomial_coeffs[l - 1] + l += 1 + end + reverse!(a_unknown) + + num_coeffs_max = num_stages - 5 + + a_matrix = zeros(num_coeffs_max, 2) + a_matrix[:, 1] = c[3:(num_stages - 3)] + + a_matrix[:, 1] -= a_unknown + a_matrix[:, 2] = a_unknown + end + + # Constant/non-optimized part of the Butcher matrix + a_matrix_constant = [0.479274057836310-0.114851811257441 / cS3 0.114851811257441/cS3 + 0.1397682537005989 0.648906880894214 + 0.1830127018922191 0.028312163512968] + + return a_matrix, a_matrix_constant, c, dt_opt +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 4 +# using provided values of coefficients a in A-matrix of Butcher tableau function compute_PairedExplicitRK4_butcher_tableau(num_stages, - base_path_a_coeffs::AbstractString; - c_const = 1.0) # Default value for best internal stability - c = c_const * ones(num_stages) # Use same abscissae for free coefficients + base_path_a_coeffs::AbstractString, + cS3) + c = ones(num_stages) # Best internal stability properties c[1] = 0.0 - cS3 = c_const c[num_stages - 3] = cS3 c[num_stages - 2] = 0.479274057836310 c[num_stages - 1] = sqrt(3) / 6 + 0.5 @@ -45,7 +102,9 @@ end @doc raw""" PairedExplicitRK4(num_stages, base_path_a_coeffs::AbstractString, dt_opt = nothing; - c_const = 1.0f0) + cS3 = 1.0f0) + PairedExplicitRK4(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, cS3 = 1.0f0) Parameters: - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method. @@ -54,8 +113,12 @@ end The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks. - `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used. - - `c_const` (`Float64`, optional): Value of abscissae $c$ in the Butcher tableau for the optimized coefficients. - Default is 1.0. + - `tspan`: Time span of the simulation. + - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. + - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the + equation has been semidiscretized. + - `cS3` (`Float64`, optional): Value of $c_{S-3}$ in the Butcher tableau, where + $S$ is the number of stages. Default is 1.0f0. The following structures and methods provide an implementation of the fourth-order paired explicit Runge-Kutta (P-ERK) method @@ -68,7 +131,7 @@ The method has been proposed in mutable struct PairedExplicitRK4 <: AbstractPairedExplicitRKSingle const num_stages::Int # S - a_matrix::Matrix{Float64} + a_matrix::Union{Matrix{Float64}, Nothing} # Nothing for S = 5 # This part of the Butcher array matrix A is constant for all PERK methods, i.e., # regardless of the optimized coefficients. a_matrix_constant::Matrix{Float64} @@ -79,11 +142,33 @@ end # struct PairedExplicitRK4 # Constructor for previously computed A Coeffs function PairedExplicitRK4(num_stages, base_path_a_coeffs::AbstractString, dt_opt = nothing; - c_const = 1.0f0) + cS3 = 1.0f0) # Default value for best internal stability + @assert num_stages>=5 "PERK4 requires at least five stages" a_matrix, a_matrix_constant, c = compute_PairedExplicitRK4_butcher_tableau(num_stages, - base_path_a_coeffs; - c_const) + base_path_a_coeffs, + cS3) + + return PairedExplicitRK4(num_stages, a_matrix, a_matrix_constant, c, dt_opt) +end + +# Constructor that computes Butcher matrix A coefficients from a semidiscretization +function PairedExplicitRK4(num_stages, tspan, semi::AbstractSemidiscretization; + verbose = false, cS3 = 1.0f0) + @assert num_stages>=5 "PERK4 requires at least five stages" + eig_vals = eigvals(jacobian_ad_forward(semi)) + + return PairedExplicitRK4(num_stages, tspan, eig_vals; verbose, cS3) +end +# Constructor that calculates the coefficients with polynomial optimizer from a list of eigenvalues +function PairedExplicitRK4(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, cS3 = 1.0f0) + @assert num_stages>=5 "PERK4 requires at least five stages" + a_matrix, a_matrix_constant, c, dt_opt = compute_PairedExplicitRK4_butcher_tableau(num_stages, + tspan, + eig_vals; + verbose, + cS3) return PairedExplicitRK4(num_stages, a_matrix, a_matrix_constant, c, dt_opt) end 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 bfd53ba2eaf..430016d1270 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -25,4 +25,6 @@ end # such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package # extension or by the Convex and ECOS-specific code loaded by Requires.jl function undo_normalization! end +function undo_normalization_PERK4! end function bisect_stability_polynomial end +function bisect_stability_polynomial_PERK4 end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index f78f127f434..7c3a97010cd 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -593,6 +593,34 @@ end end end +@trixi_testset "elixir_euler_vortex_perk4.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_perk4.jl"), + l2=[ + 0.0001846244731283424, + 0.00042537910268029285, + 0.0003724909264689687, + 0.0026689613797051493 + ], + linf=[ + 0.0025031072787504716, + 0.009266316022570331, + 0.009876399281272374, + 0.0306915591360557 + ]) + # 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) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end + @trixi_testset "elixir_euler_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), l2=[ diff --git a/test/test_tree_1d_hypdiff.jl b/test/test_tree_1d_hypdiff.jl index cd570c16708..b2515d20a01 100644 --- a/test/test_tree_1d_hypdiff.jl +++ b/test/test_tree_1d_hypdiff.jl @@ -25,6 +25,25 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") end end +@trixi_testset "elixir_hypdiff_nonperiodic_perk4.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic_perk4.jl"), + l2=[1.3655114989569954e-7, 1.020034502220849e-6], + linf=[7.173289721662535e-7, 4.507115265006689e-6], + atol=2.5e-13) + # 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) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end + @trixi_testset "elixir_hypdiff_harmonic_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_harmonic_nonperiodic.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index b0925f31cd8..7c59b1692e1 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1749,6 +1749,27 @@ end 0.31168238866709846 0.18831761133290154], atol = 1e-13) end +@testset "PERK Single p4 Constructors" begin + path_coeff_file = mktempdir() + Trixi.download("https://gist.githubusercontent.com/warisa-r/8d93f6a3ae0635e13b9f51ee32ab7fff/raw/54dc5b14be9288e186b745facb5bbcb04d1476f8/EigenvalueList_Refined2.txt", + joinpath(path_coeff_file, "spectrum.txt")) + + eig_vals = readdlm(joinpath(path_coeff_file, "spectrum.txt"), ComplexF64) + tspan = (0.0, 1.0) + ode_algorithm = Trixi.PairedExplicitRK4(14, tspan, vec(eig_vals)) + + @test isapprox(ode_algorithm.a_matrix, + [0.9935760056654522 0.006423994334547779 + 0.984991598524171 0.01500840147582901 + 0.9731962964227893 0.026803703577210732 + 0.9564649429772978 0.04353505702270225 + 0.9319644395990909 0.0680355604009091 + 0.8955295433261026 0.10447045667389748 + 0.8444449101874965 0.15555508981250354 + 0.7923618582881918 0.20763814171180817 + 0.7723161199637355 0.2276838800362645], atol = 1e-13) +end + @testset "Sutherlands Law" begin function mu(u, equations) T_ref = 291.15 From 749a5059264b9b6414ab37ec25d989221c7a3602 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Nov 2024 14:50:45 +0100 Subject: [PATCH 03/61] typos --- examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl index 676cc2a9da6..2bbf930932c 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl @@ -44,7 +44,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -# Construct third order paired explicit Runge-Kutta method with 8 stages for given simulation setup. +# Construct fourth order paired explicit Runge-Kutta method with 11 stages for given simulation setup. # Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used # in calculating the polynomial coefficients in the ODE algorithm. ode_algorithm = Trixi.PairedExplicitRK4(11, tspan, semi) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index dcce0409924..0a2c38878e1 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -105,6 +105,8 @@ end cS3 = 1.0f0) PairedExplicitRK4(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, cS3 = 1.0f0) + PairedExplicitRK4(num_stages, tspan, eig_vals::Vector{ComplexF64}; + verbose = false, cS3 = 1.0f0) Parameters: - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method. From 8ce78250d6997b4b109f26e0dda7b7186161eb15 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Nov 2024 14:54:55 +0100 Subject: [PATCH 04/61] comments --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 2 +- examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 0f918bed16c..ce4f0aceacf 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -1,5 +1,5 @@ -using OrdinaryDiffEq # Requiref for `CallbackSet` +using OrdinaryDiffEq # Required for `CallbackSet` using Trixi # Ratio of specific heats diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl index 2bbf930932c..b43ab3c416d 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl @@ -1,5 +1,5 @@ -using OrdinaryDiffEq +using OrdinaryDiffEq # Required for `CallbackSet` using Trixi # Convex and ECOS are imported because they are used for finding the optimal time step and optimal From 072693a1ae8dc21e5cb59fd055bfac7a19017178 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 8 Nov 2024 15:16:40 +0100 Subject: [PATCH 05/61] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl --- .../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 8ce2bb7b4ed..c4aab7441f3 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -249,9 +249,6 @@ end @threaded for i in eachindex(integrator.du) integrator.k1[i] = integrator.du[i] * integrator.dt - end - - @threaded for i in eachindex(integrator.u) integrator.u_tmp[i] = integrator.u[i] + c[2] * integrator.k1[i] end end From 5266a6c68476af295d9ea764e1a18a891a857bca Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 1 Dec 2024 16:04:31 +0100 Subject: [PATCH 06/61] Make `PairedExplicitRK2` constructors consistent --- NEWS.md | 8 ++++++++ .../paired_explicit_runge_kutta/methods_PERK2.jl | 4 ++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 2f71c45e57e..6c7c4fe8362 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.10 from v0.9.x + +#### Changed + +- The `PairedExplicitRK2` constructor with second argument `base_path_monomial_coeffs::AbstractString` requires + now `dt_opt`, `b1`, `cS` to be given as keyword arguments ([#2184]). + Previously, those where standard function parameters, in the same order as listed above. + ## Changes in the v0.9 lifecycle #### Added 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 2451680a505..d75ee041e3d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -98,7 +98,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; dt_opt = nothing, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -143,7 +143,7 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file -function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; dt_opt = nothing, bS = 1.0, cS = 0.5) # If the user has the monomial coefficients, they also must have the optimal time step From ba49773fbb091d1a30b08ffbc319799cbda0e165 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 1 Dec 2024 16:06:37 +0100 Subject: [PATCH 07/61] Make PERK2 constructor consistent --- NEWS.md | 8 -------- .../paired_explicit_runge_kutta/methods_PERK2.jl | 4 ++-- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6c7c4fe8362..2f71c45e57e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,14 +4,6 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -## Changes when updating to v0.10 from v0.9.x - -#### Changed - -- The `PairedExplicitRK2` constructor with second argument `base_path_monomial_coeffs::AbstractString` requires - now `dt_opt`, `b1`, `cS` to be given as keyword arguments ([#2184]). - Previously, those where standard function parameters, in the same order as listed above. - ## Changes in the v0.9 lifecycle #### Added 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 d75ee041e3d..2451680a505 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -98,7 +98,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; dt_opt = nothing, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -143,7 +143,7 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file -function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing, bS = 1.0, cS = 0.5) # If the user has the monomial coefficients, they also must have the optimal time step From c486b8babeb5c95795cd7c592edef08e4ebb1782 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 1 Dec 2024 16:09:01 +0100 Subject: [PATCH 08/61] changes --- NEWS.md | 8 ++++++++ .../paired_explicit_runge_kutta/methods_PERK2.jl | 4 ++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 2f71c45e57e..6c7c4fe8362 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.10 from v0.9.x + +#### Changed + +- The `PairedExplicitRK2` constructor with second argument `base_path_monomial_coeffs::AbstractString` requires + now `dt_opt`, `b1`, `cS` to be given as keyword arguments ([#2184]). + Previously, those where standard function parameters, in the same order as listed above. + ## Changes in the v0.9 lifecycle #### Added 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 2451680a505..d75ee041e3d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -98,7 +98,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; dt_opt = nothing, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -143,7 +143,7 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file -function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString; dt_opt = nothing, bS = 1.0, cS = 0.5) # If the user has the monomial coefficients, they also must have the optimal time step From cd2e86827ae6b6f7c1e2dec76c039cd0db32b641 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 2 Dec 2024 17:46:34 +0100 Subject: [PATCH 09/61] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 6c7c4fe8362..db06e55d249 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ for human readability. #### Changed - The `PairedExplicitRK2` constructor with second argument `base_path_monomial_coeffs::AbstractString` requires - now `dt_opt`, `b1`, `cS` to be given as keyword arguments ([#2184]). + now `dt_opt`, `bS`, `cS` to be given as keyword arguments ([#2184]). Previously, those where standard function parameters, in the same order as listed above. ## Changes in the v0.9 lifecycle From fb21c9eaa2c304281f9329ef8fb194fb7a960199 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 09:41:45 +0100 Subject: [PATCH 10/61] version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 87be80456c3..ac1af2c49eb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.4-DEV" +version = "0.9.8-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From fb029167ae68c271ff2b134189f64fe538d5658f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 09:49:14 +0100 Subject: [PATCH 11/61] cosmetics --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 4 ++-- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 3 files changed, 4 insertions(+), 4 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 1ed89a124e2..156785ae40c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -116,9 +116,9 @@ end equation has been semidiscretized. - `verbose` (`Bool`, optional): Verbosity flag, default is false. - `bS` (`Float64`, optional): Value of $b_S$ in the Butcher tableau, where - $S$ is the number of stages. Default is 1.0. + $S$ is the number of stages. Default is `1.0`. - `cS` (`Float64`, optional): Value of $c_S$ in the Butcher tableau, where - $S$ is the number of stages. Default is 0.5. + $S$ is the number of stages. Default is `0.5`. The following structures and methods provide a minimal implementation of the second-order paired explicit Runge-Kutta (PERK) method 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 853f5170337..ed59e157528 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -118,7 +118,7 @@ end equation has been semidiscretized. - `verbose` (`Bool`, optional): Verbosity flag, default is false. - `cS2` (`Float64`, optional): Value of $c_{S-2}$ in the Butcher tableau, where - $S$ is the number of stages. Default is 1.0f0. + $S$ is the number of stages. Default is `1.0f0`. The following structures and methods provide an implementation of the third-order paired explicit Runge-Kutta (PERK) method diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 7895585143d..bda60ed3c36 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -115,7 +115,7 @@ end - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the equation has been semidiscretized. - `cS3` (`Float64`, optional): Value of $c_{S-3}$ in the Butcher tableau, where - $S$ is the number of stages. Default is 1.0f0. + $S$ is the number of stages. Default is `1.0f0`. The following structures and methods provide an implementation of the fourth-order paired explicit Runge-Kutta (P-ERK) method From 3e231222fceac19b8aa9318d2d1ce193e8a6a2e8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 10:32:09 +0100 Subject: [PATCH 12/61] remove func --- ext/TrixiConvexECOSExt.jl | 33 +++++++++---------- .../methods_PERK2.jl | 1 - .../methods_PERK3.jl | 2 -- .../polynomial_optimizer.jl | 1 - 4 files changed, 15 insertions(+), 22 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 4e9f69f3713..fa0ed53319d 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -15,7 +15,7 @@ end using LinearAlgebra: eigvals # Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, undo_normalization!, undo_normalization_PERK4!, +using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, bisect_stability_polynomial_PERK4, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). @@ -27,15 +27,10 @@ using Trixi: Trixi, undo_normalization!, undo_normalization_PERK4!, # Undo normalization of stability polynomial coefficients by index factorial # relative to consistency order. -function Trixi.undo_normalization!(gamma_opt, consistency_order, num_stage_evals) - for k in 1:(num_stage_evals - consistency_order) - gamma_opt[k] /= factorial(k+consistency_order) - end -end - -function Trixi.undo_normalization_PERK4!(gamma_opt, num_stage_evals) - for k in 1:(num_stage_evals - 5) - gamma_opt[k] /= factorial(k + 4) +function Trixi.undo_normalization!(gamma_opt, num_stage_evals, + num_reduced_coeffs, fac_offset) + for k in 1:(num_stage_evals - num_reduced_coeffs) + gamma_opt[k] /= factorial(k + fac_offset) end end @@ -209,7 +204,8 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, gamma_opt = [gamma_opt] end - undo_normalization!(gamma_opt, consistency_order, num_stage_evals) + undo_normalization!(gamma_opt, num_stage_evals, + consistency_order, consistency_order) return gamma_opt, dt end @@ -283,18 +279,19 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, println("Concluded stability polynomial optimization \n") end - gamma_opt = [] if num_stage_evals > 5 gamma_opt = evaluate(gamma) + else + gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. + end - # Catch case S = 6 (only one opt. variable) - if isa(gamma_opt, Number) - gamma_opt = [gamma_opt] - end - - undo_normalization_PERK4!(gamma_opt, num_stage_evals) + # Catch case S = 6 (only one opt. variable) + if isa(gamma_opt, Number) + gamma_opt = [gamma_opt] end + undo_normalization!(gamma_opt, num_stage_evals, 5, 4) + return gamma_opt, dt 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 156785ae40c..0fe709d12af 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -51,7 +51,6 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, eig_vals; verbose) if num_stages != consistency_order - undo_normalization!(monomial_coeffs, consistency_order, num_stages) num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) 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 ed59e157528..6c402d8ebfa 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -49,8 +49,6 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, if num_stages == consistency_order a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau else - undo_normalization!(monomial_coeffs, consistency_order, num_stages) - # Solve the nonlinear system of equations from monomial coefficient and # Butcher array abscissae c to find Butcher matrix A # This function is extended in TrixiNLsolveExt.jl 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 430016d1270..6658f1ac4e6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -25,6 +25,5 @@ end # such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package # extension or by the Convex and ECOS-specific code loaded by Requires.jl function undo_normalization! end -function undo_normalization_PERK4! end function bisect_stability_polynomial end function bisect_stability_polynomial_PERK4 end From 4cea3a9707bc7186aa6cd457fa8e344e153ac941 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 10:35:18 +0100 Subject: [PATCH 13/61] remove unused func --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 11 ----------- 1 file changed, 11 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 0fe709d12af..f06fb90c66b 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -241,17 +241,6 @@ function init(ode::ODEProblem, alg::PairedExplicitRK2; return integrator end -# Function that computes the first stage of a general P-ERK method -# (in fact every explicit Runge-Kutta method) -@inline function k1!(integrator::AbstractPairedExplicitRKIntegrator, p, c) - integrator.f(integrator.du, integrator.u, p, integrator.t) - - @threaded for i in eachindex(integrator.du) - integrator.k1[i] = integrator.du[i] * integrator.dt - integrator.u_tmp[i] = integrator.u[i] + c[2] * integrator.k1[i] - end -end - function step!(integrator::PairedExplicitRK2Integrator) @unpack prob = integrator.sol @unpack alg = integrator From 69a1e5e712f5a28363d156bf1fba3fe104dc1e50 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 10:55:20 +0100 Subject: [PATCH 14/61] make timestep computable perk4 s=5 --- ext/TrixiConvexECOSExt.jl | 6 ++++- .../methods_PERK2.jl | 12 +++++----- .../methods_PERK3.jl | 6 ++--- .../methods_PERK4.jl | 23 +++++++++---------- 4 files changed, 25 insertions(+), 22 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index fa0ed53319d..07647af0658 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -102,7 +102,11 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, end # For optimization only the maximum is relevant - return maximum(abs(pnoms)) + if num_stage_evals == 5 + return maximum(abs.(pnoms)) # If there is no variable to optimize, we need to use the broadcast operator. + else + return maximum(abs(pnoms)) + end end #= 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 f06fb90c66b..01964219eae 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -32,9 +32,9 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, stage_scaling_factors = bS * reverse(c[2:(end - 1)]) # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) - coeffs_max = num_stages - 2 + num_coeffs_max = num_stages - 2 - a_matrix = zeros(2, coeffs_max) + a_matrix = zeros(2, num_coeffs_max) a_matrix[1, :] = c[3:end] consistency_order = 2 @@ -52,7 +52,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, if num_stages != consistency_order num_monomial_coeffs = length(monomial_coeffs) - @assert num_monomial_coeffs == coeffs_max + @assert num_monomial_coeffs == num_coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) a_matrix[1, :] -= A a_matrix[2, :] = A @@ -74,9 +74,9 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, stage_scaling_factors = bS * reverse(c[2:(end - 1)]) # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) - coeffs_max = num_stages - 2 + num_coeffs_max = num_stages - 2 - a_matrix = zeros(2, coeffs_max) + a_matrix = zeros(2, num_coeffs_max) a_matrix[1, :] = c[3:end] path_monomial_coeffs = joinpath(base_path_monomial_coeffs, @@ -86,7 +86,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, monomial_coeffs = readdlm(path_monomial_coeffs, Float64) num_monomial_coeffs = size(monomial_coeffs, 1) - @assert num_monomial_coeffs == coeffs_max + @assert num_monomial_coeffs == num_coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) a_matrix[1, :] -= A 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 6c402d8ebfa..c4fb048c6d1 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -75,9 +75,9 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, c = compute_c_coeffs(num_stages, cS2) # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) - a_coeffs_max = num_stages - 2 + num_a_coeffs_max = num_stages - 2 - a_matrix = zeros(2, a_coeffs_max) + a_matrix = zeros(2, num_a_coeffs_max) a_matrix[1, :] = c[3:end] path_a_coeffs = joinpath(base_path_a_coeffs, @@ -87,7 +87,7 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, a_coeffs = readdlm(path_a_coeffs, Float64) num_a_coeffs = size(a_coeffs, 1) - @assert num_a_coeffs == a_coeffs_max + @assert num_a_coeffs == num_a_coeffs_max # Fill A-matrix in PERK style a_matrix[1, :] -= a_coeffs a_matrix[2, :] = a_coeffs diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index bda60ed3c36..7b9c48b8ac6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -19,22 +19,21 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, c[num_stages] = -sqrt(3) / 6 + 0.5 num_coeffs_max = num_stages - 5 + a_matrix = zeros(2, num_coeffs_max) - a_matrix = nothing - dt_opt = nothing - if num_stages > 5 - # Calculate coefficients of the stability polynomial in monomial form - dtmax = tspan[2] - tspan[1] - dteps = 1.0f-9 + # Calculate coefficients of the stability polynomial in monomial form + dtmax = tspan[2] - tspan[1] + dteps = 1.0f-9 - num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) - monomial_coeffs, dt_opt = bisect_stability_polynomial_PERK4(num_eig_vals, - num_stages, - dtmax, dteps, - eig_vals, cS3; - verbose) + monomial_coeffs, dt_opt = bisect_stability_polynomial_PERK4(num_eig_vals, + num_stages, + dtmax, dteps, + eig_vals, cS3; + verbose) + if num_stages > 5 a_unknown = copy(monomial_coeffs) l = 2 for _ in 5:(num_stages - 2) From 676ebace8595bd9ea99c934c789262bfd33b3fb8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 11:14:24 +0100 Subject: [PATCH 15/61] compute dt for no opt vars --- ext/TrixiConvexECOSExt.jl | 100 ++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 43 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 07647af0658..ef05b610806 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -164,27 +164,35 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, end end - # Use last optimal values for gamma in (potentially) next iteration - problem = minimize(stability_polynomials!(pnoms, consistency_order, - num_stage_evals, - normalized_powered_eigvals_scaled, - gamma)) - - solve!(problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - "feastol" => 1e-9, - "abstol" => 1e-9, - "reltol" => 1e-9, - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent = true) - - abs_p = problem.optval + # Check if there are variables to optimize + if num_stage_evals - consistency_order > 0 + # Use last optimal values for gamma in (potentially) next iteration + problem = minimize(stability_polynomials!(pnoms, consistency_order, + num_stage_evals, + normalized_powered_eigvals_scaled, + gamma)) + + solve!(problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + "feastol" => 1e-9, + "abstol" => 1e-9, + "reltol" => 1e-9, + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent = true) + + abs_p = problem.optval + else + abs_p = stability_polynomials!(pnoms, consistency_order, + num_stage_evals, + normalized_powered_eigvals_scaled, + gamma) + end if abs_p < 1 dtmin = dt @@ -197,7 +205,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, println("Concluded stability polynomial optimization \n") end - if consistency_order - num_stage_evals != 0 + if num_stage_evals - consistency_order > 0 gamma_opt = evaluate(gamma) else gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. @@ -250,27 +258,33 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, while dtmax - dtmin > dteps dt = 0.5 * (dtmax + dtmin) - # Use last optimal values for gamma in (potentially) next iteration - problem = minimize(stability_polynomials_PERK4!(pnoms, - num_stage_evals, - normalized_powered_eigvals, - gamma, dt, cS3)) - - solve!(problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - "feastol" => 1e-9, - "abstol" => 1e-9, - "reltol" => 1e-9, - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent = true) - - abs_p = problem.optval + if num_stage_evals > 5 + # Use last optimal values for gamma in (potentially) next iteration + problem = minimize(stability_polynomials_PERK4!(pnoms, + num_stage_evals, + normalized_powered_eigvals, + gamma, dt, cS3)) + + solve!(problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + "feastol" => 1e-9, + "abstol" => 1e-9, + "reltol" => 1e-9, + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent = true) + + abs_p = problem.optval + else + abs_p = stability_polynomials_PERK4!(pnoms, num_stage_evals, + normalized_powered_eigvals, + gamma, dt, cS3) + end if abs_p < 1 dtmin = dt From a1df22d19231a17e89d06a5403a090a3a7dccf11 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 11:15:08 +0100 Subject: [PATCH 16/61] float --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 6925bae0630..1fccc9c8258 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -56,7 +56,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE end initial_condition = initial_condition_isentropic_vortex -EdgeLength = 20 +EdgeLength = 20.0 N_passes = 1 T_end = EdgeLength * N_passes From 7b8b16663c581c92898451e4c1fbe47ecdde98b4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 3 Dec 2024 18:07:13 +0100 Subject: [PATCH 17/61] mention perk4 --- docs/src/time_integration.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index 3cdb7fa75ae..d493c18657c 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -212,4 +212,5 @@ Then, the stable CFL number can be computed as described above. ##### Single/Standalone methods - [`Trixi.PairedExplicitRK2`](@ref): Second-order PERK method with at least two stages. -- [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages. \ No newline at end of file +- [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages. +- [`Trixi.PairedExplicitRK4`](@ref): Fourth-order PERK method with at least _five_ stages. \ No newline at end of file From 921a3b929ac477dd6879cd151e7c72d3407a9653 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 5 Dec 2024 08:09:17 +0100 Subject: [PATCH 18/61] Apply suggestions from code review Co-authored-by: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 7b9c48b8ac6..ec5e4eb6e29 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -233,13 +233,13 @@ end @inline function last_three_stages!(integrator::PairedExplicitRK4Integrator, p, alg) for stage in 1:2 - @threaded for u_ind in eachindex(integrator.u) - integrator.u_tmp[u_ind] = integrator.u[u_ind] + + @threaded for i in eachindex(integrator.u) + integrator.u_tmp[u_ind] = integrator.u[i] + integrator.dt * (alg.a_matrix_constant[1, stage] * - integrator.k1[u_ind] + + integrator.k1[i] + alg.a_matrix_constant[2, stage] * - integrator.du[u_ind]) + integrator.du[i]) end integrator.f(integrator.du, integrator.u_tmp, p, @@ -263,11 +263,11 @@ end integrator.f(integrator.du, integrator.u_tmp, p, integrator.t + alg.c[alg.num_stages] * integrator.dt) - @threaded for u_ind in eachindex(integrator.u) + @threaded for i in eachindex(integrator.u) # Note that 'k1' carries the values of K_{S-1} # and that we construct 'K_S' "in-place" from 'integrator.du' integrator.u[u_ind] += 0.5 * integrator.dt * - (integrator.k1[u_ind] + integrator.du[u_ind]) + (integrator.k1[i] + integrator.du[i]) end end From 5c4fce264fbec87e882a3c6ea347ec39d70f2d49 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:08:53 +0100 Subject: [PATCH 19/61] prevent readdlm error for S=2 --- .../methods_PERK2.jl | 26 ++++++++++--------- .../methods_PERK3.jl | 2 +- 2 files changed, 15 insertions(+), 13 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 d75ee041e3d..edf67809abc 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -33,8 +33,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) coeffs_max = num_stages - 2 - a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] consistency_order = 2 @@ -77,22 +77,24 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) coeffs_max = num_stages - 2 - a_matrix = zeros(coeffs_max, 2) - a_matrix[:, 1] = c[3:end] - path_monomial_coeffs = joinpath(base_path_monomial_coeffs, - "gamma_" * string(num_stages) * ".txt") + if coeffs_max > 0 + a_matrix[:, 1] = c[3:end] + + path_monomial_coeffs = joinpath(base_path_monomial_coeffs, + "gamma_" * string(num_stages) * ".txt") - @assert isfile(path_monomial_coeffs) "Couldn't find file" - monomial_coeffs = readdlm(path_monomial_coeffs, Float64) - num_monomial_coeffs = size(monomial_coeffs, 1) + @assert isfile(path_monomial_coeffs) "Couldn't find file" + monomial_coeffs = readdlm(path_monomial_coeffs, Float64) + num_monomial_coeffs = size(monomial_coeffs, 1) - @assert num_monomial_coeffs == coeffs_max - A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) + @assert num_monomial_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + end return a_matrix, c 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 02bc3eeba36..2e886514184 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -79,8 +79,8 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) a_coeffs_max = num_stages - 2 - a_matrix = zeros(a_coeffs_max, 2) + a_matrix[:, 1] = c[3:end] path_a_coeffs = joinpath(base_path_a_coeffs, From 7728de06d970228c07f42e5d3e582b5449ba9ae7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:37:24 +0100 Subject: [PATCH 20/61] Warisa's comments --- ext/TrixiConvexECOSExt.jl | 78 ++++++++++--------- .../methods_PERK3.jl | 6 +- .../methods_PERK4.jl | 42 +++++----- 3 files changed, 67 insertions(+), 59 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index ef05b610806..0479237f12a 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -34,14 +34,9 @@ function Trixi.undo_normalization!(gamma_opt, num_stage_evals, end end -# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency -# order, including contributions from free coefficients for higher orders, and -# return the maximum absolute value -function stability_polynomials!(pnoms, consistency_order, num_stage_evals, - normalized_powered_eigvals_scaled, - gamma) - num_eig_vals = length(pnoms) - +@inline function stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, + normalized_powered_eigvals_scaled, + consistency_order) # Initialize with zero'th order (z^0) coefficient for i in 1:num_eig_vals pnoms[i] = 1.0 @@ -53,6 +48,19 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals, pnoms[i] += normalized_powered_eigvals_scaled[i, k] end end +end + +# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency +# order, including contributions from free coefficients for higher orders, and +# return the maximum absolute value +function stability_polynomials!(pnoms, consistency_order, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma) + stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, + normalized_powered_eigvals_scaled, + consistency_order) # Contribution from free coefficients for k in (consistency_order + 1):num_stage_evals @@ -69,26 +77,17 @@ end # Specialized form of the stability polynomials for fourth-order PERK schemes. function stability_polynomials_PERK4!(pnoms, num_stage_evals, + num_eig_vals, normalized_powered_eigvals, - gamma, dt, cS3) - num_eig_vals = length(pnoms) - + gamma, + dt, cS3) # Constants arising from the particular form of Butcher tableau chosen for the 4th order PERK methods k1 = 0.001055026310046423 / cS3 k2 = 0.03726406530405851 / cS3 # Note: `cS3` = c_{S-3} is in principle free, while the other abscissae are fixed to 1.0 - # Initialize with zero'th order (z^0) coefficient - for i in 1:num_eig_vals - pnoms[i] = 1.0 - end - - # First `consistency_order` = 4 terms of the exponential - for k in 1:4 - for i in 1:num_eig_vals - pnoms[i] += dt^k * normalized_powered_eigvals[i, k] - end - end + stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, normalized_powered_eigvals, + 4) # "Fixed" term due to choice of the PERK4 Butcher tableau # Required to un-do the normalization of the eigenvalues here @@ -109,6 +108,17 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, end end +@inline function normalized_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, 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 +end + #= The following structures and methods provide a simplified implementation to discover optimal stability polynomial for a given set of `eig_vals` @@ -136,13 +146,9 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, gamma = Variable(num_stage_evals - consistency_order) normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - - for j in 1:num_stage_evals - fac_j = factorial(j) - for i in 1:num_eig_vals - normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j - end - end + normalized_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, eig_vals, + num_stage_evals) normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) @@ -169,6 +175,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, # Use last optimal values for gamma in (potentially) next iteration problem = minimize(stability_polynomials!(pnoms, consistency_order, num_stage_evals, + num_eig_vals, normalized_powered_eigvals_scaled, gamma)) @@ -190,6 +197,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, else abs_p = stability_polynomials!(pnoms, consistency_order, num_stage_evals, + num_eig_vals, normalized_powered_eigvals_scaled, gamma) end @@ -242,13 +250,9 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, gamma = Variable(num_stage_evals - 5) 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_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, eig_vals, + num_stage_evals) if verbose println("Start optimization of stability polynomial \n") @@ -262,6 +266,7 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, # Use last optimal values for gamma in (potentially) next iteration problem = minimize(stability_polynomials_PERK4!(pnoms, num_stage_evals, + num_eig_vals, normalized_powered_eigvals, gamma, dt, cS3)) @@ -282,6 +287,7 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, abs_p = problem.optval else abs_p = stability_polynomials_PERK4!(pnoms, num_stage_evals, + num_eig_vals, normalized_powered_eigvals, gamma, dt, cS3) 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 c4fb048c6d1..d4371b61896 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -6,7 +6,7 @@ #! format: noindent # Initialize Butcher array abscissae c for PairedExplicitRK3 based on SSPRK33 base method -function compute_c_coeffs(num_stages, cS2) +function PERK3_compute_c_coeffs(num_stages, cS2) c = zeros(eltype(cS2), num_stages) # Last timesteps as for SSPRK33, see motivation in Section 3.3 of @@ -28,7 +28,7 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, cS2) # Initialize array of c - c = compute_c_coeffs(num_stages, cS2) + c = PERK3_compute_c_coeffs(num_stages, cS2) # Initialize the array of our solution a_unknown = zeros(num_stages - 2) @@ -72,7 +72,7 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, cS2) # Initialize array of c - c = compute_c_coeffs(num_stages, cS2) + c = PERK3_compute_c_coeffs(num_stages, cS2) # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) num_a_coeffs_max = num_stages - 2 diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 7b9c48b8ac6..9296eabc89d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -5,11 +5,7 @@ @muladd begin #! format: noindent -# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 4 -# using a list of eigenvalues -function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, - eig_vals::Vector{ComplexF64}; - verbose = false, cS3) +function PERK4_compute_c_coeffs(num_stages, cS3) c = ones(num_stages) # Best internal stability properties c[1] = 0.0 @@ -18,6 +14,22 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, c[num_stages - 1] = sqrt(3) / 6 + 0.5 c[num_stages] = -sqrt(3) / 6 + 0.5 + return c +end + +# Constant/non-optimized part of the Butcher matrix +function PERK4_a_matrix_constant(cS3) + return [(0.479274057836310-(0.114851811257441 / cS3)) 0.1397682537005989 0.1830127018922191 + 0.114851811257441/cS3 0.648906880894214 0.028312163512968] +end + +# Compute the Butcher tableau for a paired explicit Runge-Kutta method order 4 +# using a list of eigenvalues +function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, + eig_vals::Vector{ComplexF64}; + verbose = false, cS3) + c = PERK4_compute_c_coeffs(num_stages, cS3) + num_coeffs_max = num_stages - 5 a_matrix = zeros(2, num_coeffs_max) @@ -49,9 +61,7 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, a_matrix[2, :] = a_unknown end - # Constant/non-optimized part of the Butcher matrix - a_matrix_constant = [(0.479274057836310-(0.114851811257441 / cS3)) 0.1397682537005989 0.1830127018922191 - 0.114851811257441/cS3 0.648906880894214 0.028312163512968] + a_matrix_constant = PERK4_a_matrix_constant(cS3) return a_matrix, a_matrix_constant, c, dt_opt end @@ -61,13 +71,7 @@ end function compute_PairedExplicitRK4_butcher_tableau(num_stages, base_path_a_coeffs::AbstractString, cS3) - c = ones(num_stages) # Best internal stability properties - c[1] = 0.0 - - c[num_stages - 3] = cS3 - c[num_stages - 2] = 0.479274057836310 - c[num_stages - 1] = sqrt(3) / 6 + 0.5 - c[num_stages] = -sqrt(3) / 6 + 0.5 + c = PERK4_compute_c_coeffs(num_stages, cS3) num_coeffs_max = num_stages - 5 @@ -87,9 +91,7 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, a_matrix[2, :] = a_coeffs end - # Constant/non-optimized part of the Butcher matrix - a_matrix_constant = [(0.479274057836310-(0.114851811257441 / cS3)) 0.1397682537005989 0.1830127018922191 - 0.114851811257441/cS3 0.648906880894214 0.028312163512968] + a_matrix_constant = PERK4_a_matrix_constant(cS3) return a_matrix, a_matrix_constant, c end @@ -231,7 +233,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK4; return integrator end -@inline function last_three_stages!(integrator::PairedExplicitRK4Integrator, p, alg) +@inline function PERK4_KS2_to_KS!(integrator::PairedExplicitRK4Integrator, p, alg) for stage in 1:2 @threaded for u_ind in eachindex(integrator.u) integrator.u_tmp[u_ind] = integrator.u[u_ind] + @@ -300,7 +302,7 @@ function step!(integrator::PairedExplicitRK4Integrator) PERK_ki!(integrator, prob.p, alg, stage) end - last_three_stages!(integrator, prob.p, alg) + PERK4_KS2_to_KS!(integrator, prob.p, alg) end integrator.iter += 1 From 2a9d8c4ebc18faa6189cea94b1875d8bc8d37264 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:43:29 +0100 Subject: [PATCH 21/61] rename --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 6383f14bc4e..d82f1e5cfd1 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -233,7 +233,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK4; return integrator end -@inline function PERK4_KS2_to_KS!(integrator::PairedExplicitRK4Integrator, p, alg) +@inline function PERK4_kS2_to_kS!(integrator::PairedExplicitRK4Integrator, p, alg) for stage in 1:2 @threaded for i in eachindex(integrator.u) integrator.u_tmp[u_ind] = integrator.u[i] + @@ -302,7 +302,7 @@ function step!(integrator::PairedExplicitRK4Integrator) PERK_ki!(integrator, prob.p, alg, stage) end - PERK4_KS2_to_KS!(integrator, prob.p, alg) + PERK4_kS2_to_kS!(integrator, prob.p, alg) end integrator.iter += 1 From 6badb45aa8638fdd6f144eaf7ed40c636ddf7984 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:47:16 +0100 Subject: [PATCH 22/61] PERK2 compute c coeffs --- .../methods_PERK2.jl | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 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 01964219eae..75b443bef11 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -5,6 +5,15 @@ @muladd begin #! format: noindent +function PERK2_compute_c_coeffs(num_stages, cS2) + c = zeros(num_stages) + for k in 2:num_stages + c[k] = cS * (k - 1) / (num_stages - 1) + end + + return c +end + # Compute the coefficients of the A matrix in the Butcher tableau using # stage scaling factors and monomial coefficients function compute_a_coeffs(num_stage_evals, stage_scaling_factors, monomial_coeffs) @@ -24,11 +33,7 @@ end # using a list of eigenvalues function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, bS, cS; verbose = false) - # c Vector from Butcher Tableau (defines timestep per stage) - c = zeros(num_stages) - for k in 2:num_stages - c[k] = cS * (k - 1) / (num_stages - 1) - end + c = PERK2_compute_c_coeffs(num_stages, cS) stage_scaling_factors = bS * reverse(c[2:(end - 1)]) # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) @@ -66,11 +71,7 @@ end function compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs::AbstractString, bS, cS) - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(num_stages) - for k in 2:num_stages - c[k] = cS * (k - 1) / (num_stages - 1) - end + c = PERK2_compute_c_coeffs(num_stages, cS) stage_scaling_factors = bS * reverse(c[2:(end - 1)]) # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency) From b232eec85093f31d50e3ed83cb3e30b0278c876d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:53:10 +0100 Subject: [PATCH 23/61] comment --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index d82f1e5cfd1..1aaf64f10fd 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -233,6 +233,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK4; return integrator end +# Computes last three stages, i.e., i = S-2, S-1, S @inline function PERK4_kS2_to_kS!(integrator::PairedExplicitRK4Integrator, p, alg) for stage in 1:2 @threaded for i in eachindex(integrator.u) From b29c7006edb4e61ee46f777e7cd45e07eebea3bb Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 09:55:46 +0100 Subject: [PATCH 24/61] bugfix --- .../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 75b443bef11..081373749eb 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,7 @@ @muladd begin #! format: noindent -function PERK2_compute_c_coeffs(num_stages, cS2) +function PERK2_compute_c_coeffs(num_stages, cS) c = zeros(num_stages) for k in 2:num_stages c[k] = cS * (k - 1) / (num_stages - 1) From 283c3854853f38d04eb9b09ae0bb802157e6ed37 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 10 Dec 2024 09:57:42 +0100 Subject: [PATCH 25/61] Update docs/src/time_integration.md Co-authored-by: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> --- docs/src/time_integration.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index d493c18657c..9afbd5f65e4 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -213,4 +213,4 @@ Then, the stable CFL number can be computed as described above. - [`Trixi.PairedExplicitRK2`](@ref): Second-order PERK method with at least two stages. - [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages. -- [`Trixi.PairedExplicitRK4`](@ref): Fourth-order PERK method with at least _five_ stages. \ No newline at end of file +- [`Trixi.PairedExplicitRK4`](@ref): Fourth-order PERK method with at least five stages. \ No newline at end of file From c5b48a5429272b9f73936de0797cb3de33449b8a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 10:41:31 +0100 Subject: [PATCH 26/61] debug --- ext/TrixiConvexECOSExt.jl | 37 +++++++++---------- .../methods_PERK4.jl | 16 ++++---- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 0479237f12a..3ce690ab48b 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -34,33 +34,24 @@ function Trixi.undo_normalization!(gamma_opt, num_stage_evals, end end -@inline function stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, - normalized_powered_eigvals_scaled, - consistency_order) +# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency +# order, including contributions from free coefficients for higher orders, and +# return the maximum absolute value +function stability_polynomials!(pnoms, consistency_order, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma) # Initialize with zero'th order (z^0) coefficient for i in 1:num_eig_vals pnoms[i] = 1.0 end - # First `consistency_order` terms of the exponential for k in 1:consistency_order for i in 1:num_eig_vals pnoms[i] += normalized_powered_eigvals_scaled[i, k] end end -end - -# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency -# order, including contributions from free coefficients for higher orders, and -# return the maximum absolute value -function stability_polynomials!(pnoms, consistency_order, - num_stage_evals, - num_eig_vals, - normalized_powered_eigvals_scaled, - gamma) - stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, - normalized_powered_eigvals_scaled, - consistency_order) # Contribution from free coefficients for k in (consistency_order + 1):num_stage_evals @@ -86,8 +77,16 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, k2 = 0.03726406530405851 / cS3 # Note: `cS3` = c_{S-3} is in principle free, while the other abscissae are fixed to 1.0 - stability_polynomials_fixed_coeffs!(pnoms, num_eig_vals, normalized_powered_eigvals, - 4) + # Initialize with zero'th order (z^0) coefficient + for i in 1:num_eig_vals + pnoms[i] = 1.0 + end + # First `consistency_order` = 4 terms of the exponential + for k in 1:4 + for i in 1:num_eig_vals + pnoms[i] += dt^k * normalized_powered_eigvals[i, k] + end + end # "Fixed" term due to choice of the PERK4 Butcher tableau # Required to un-do the normalization of the eigenvalues here diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 1aaf64f10fd..3369bfc38be 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -237,12 +237,12 @@ end @inline function PERK4_kS2_to_kS!(integrator::PairedExplicitRK4Integrator, p, alg) for stage in 1:2 @threaded for i in eachindex(integrator.u) - integrator.u_tmp[u_ind] = integrator.u[i] + - integrator.dt * - (alg.a_matrix_constant[1, stage] * - integrator.k1[i] + - alg.a_matrix_constant[2, stage] * - integrator.du[i]) + integrator.u_tmp[i] = integrator.u[i] + + integrator.dt * + (alg.a_matrix_constant[1, stage] * + integrator.k1[i] + + alg.a_matrix_constant[2, stage] * + integrator.du[i]) end integrator.f(integrator.du, integrator.u_tmp, p, @@ -269,8 +269,8 @@ end @threaded for i in eachindex(integrator.u) # Note that 'k1' carries the values of K_{S-1} # and that we construct 'K_S' "in-place" from 'integrator.du' - integrator.u[u_ind] += 0.5 * integrator.dt * - (integrator.k1[i] + integrator.du[i]) + integrator.u[i] += 0.5 * integrator.dt * + (integrator.k1[i] + integrator.du[i]) end end From 44f4ebecaa0a4bb651482c4359b64fa2533a8678 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 10 Dec 2024 10:43:45 +0100 Subject: [PATCH 27/61] rename --- ext/TrixiConvexECOSExt.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 3ce690ab48b..85070d67103 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -107,9 +107,9 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, end end -@inline function normalized_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, - num_stage_evals) +@inline function normalize_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, eig_vals, + num_stage_evals) for j in 1:num_stage_evals fac_j = factorial(j) for i in 1:num_eig_vals @@ -145,9 +145,9 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, gamma = Variable(num_stage_evals - consistency_order) normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - normalized_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, - num_stage_evals) + normalize_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, eig_vals, + num_stage_evals) normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) @@ -249,9 +249,9 @@ function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, gamma = Variable(num_stage_evals - 5) normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - normalized_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, - num_stage_evals) + normalize_power_eigvals!(normalized_powered_eigvals, + num_eig_vals, eig_vals, + num_stage_evals) if verbose println("Start optimization of stability polynomial \n") From 22445722374d9f489a04a1387c79aba484119379 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 19 Dec 2024 08:44:18 +0100 Subject: [PATCH 28/61] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl Co-authored-by: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 3369bfc38be..cdf9d9a7c1a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -47,10 +47,8 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, if num_stages > 5 a_unknown = copy(monomial_coeffs) - l = 2 - for _ in 5:(num_stages - 2) - a_unknown[l] /= monomial_coeffs[l - 1] - l += 1 + for i in 5:(num_stages - 2) + a_unknown_1[i - 3] /= monomial_coeffs[i - 4] end reverse!(a_unknown) From 8771f9524072de4d238fcc9835a89e210eba803b Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 19 Dec 2024 08:46:23 +0100 Subject: [PATCH 29/61] Apply suggestions from code review Co-authored-by: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 1fccc9c8258..af0b45cbcea 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -56,14 +56,14 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE end initial_condition = initial_condition_isentropic_vortex -EdgeLength = 20.0 +edge_length = 20.0 N_passes = 1 -T_end = EdgeLength * N_passes +T_end = edge_length * N_passes tspan = (0.0, T_end) -coordinates_min = (-EdgeLength / 2, -EdgeLength / 2) -coordinates_max = (EdgeLength / 2, EdgeLength / 2) +coordinates_min = (-edge_length / 2, -edge_length / 2) +coordinates_max = (edge_length / 2, edge_length / 2) cells_per_dimension = (32, 32) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) From 1001f65fe9de709b41f4289cb59e127f6e20d196 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 19 Dec 2024 10:06:38 +0100 Subject: [PATCH 30/61] fix --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index cdf9d9a7c1a..920eb2d35d1 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -48,7 +48,7 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, if num_stages > 5 a_unknown = copy(monomial_coeffs) for i in 5:(num_stages - 2) - a_unknown_1[i - 3] /= monomial_coeffs[i - 4] + a_unknown[i - 3] /= monomial_coeffs[i - 4] end reverse!(a_unknown) From e9c988e8694c1ae0dd6fa7d52a192a7d1694554b Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 19 Dec 2024 10:09:03 +0100 Subject: [PATCH 31/61] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- ext/TrixiConvexECOSExt.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 85070d67103..9b3cc568dea 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -90,12 +90,12 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, # "Fixed" term due to choice of the PERK4 Butcher tableau # Required to un-do the normalization of the eigenvalues here - pnoms += k1 * dt^5 * normalized_powered_eigvals[:, 5] * factorial(5) + pnoms += k1 * dt^5 * view(normalized_powered_eigvals, :, 5) * factorial(5) # Contribution from free coefficients for k in 1:(num_stage_evals - 5) - pnoms += (k2 * dt^(k + 4) * normalized_powered_eigvals[:, k + 4] * gamma[k] + - k1 * dt^(k + 5) * normalized_powered_eigvals[:, k + 5] * gamma[k] * + pnoms += (k2 * dt^(k + 4) * view(normalized_powered_eigvals, :, k + 4) * gamma[k] + + k1 * dt^(k + 5) * view(normalized_powered_eigvals, :, k + 5) * gamma[k] * (k + 5)) end From 1044bac0a64c338764448f0c795b77c0b8f3176b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 19 Dec 2024 10:11:37 +0100 Subject: [PATCH 32/61] view --- ext/TrixiConvexECOSExt.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 9b3cc568dea..411cf31dfd9 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -55,7 +55,8 @@ function stability_polynomials!(pnoms, consistency_order, # Contribution from free coefficients for k in (consistency_order + 1):num_stage_evals - pnoms += gamma[k - consistency_order] * normalized_powered_eigvals_scaled[:, k] + pnoms += gamma[k - consistency_order] * + view(normalized_powered_eigvals_scaled[:, k]) end # For optimization only the maximum is relevant @@ -94,8 +95,10 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, # Contribution from free coefficients for k in 1:(num_stage_evals - 5) - pnoms += (k2 * dt^(k + 4) * view(normalized_powered_eigvals, :, k + 4) * gamma[k] + - k1 * dt^(k + 5) * view(normalized_powered_eigvals, :, k + 5) * gamma[k] * + pnoms += (k2 * dt^(k + 4) * view(normalized_powered_eigvals, :, k + 4) * + gamma[k] + + k1 * dt^(k + 5) * view(normalized_powered_eigvals, :, k + 5) * + gamma[k] * (k + 5)) end From dc76fa722db9484ea12a9a69357bbb5a7d06c2b1 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 19 Dec 2024 10:40:36 +0100 Subject: [PATCH 33/61] Update ext/TrixiConvexECOSExt.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- ext/TrixiConvexECOSExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 411cf31dfd9..2d4b780c7e4 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -56,7 +56,7 @@ function stability_polynomials!(pnoms, consistency_order, # Contribution from free coefficients for k in (consistency_order + 1):num_stage_evals pnoms += gamma[k - consistency_order] * - view(normalized_powered_eigvals_scaled[:, k]) + view(normalized_powered_eigvals_scaled, :, k) end # For optimization only the maximum is relevant From 28b6e2f1ef54332e0b9f2f02d5e27d12ac286b38 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 14:13:33 +0100 Subject: [PATCH 34/61] news --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index c187f229519..ce5ec36a4de 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,8 @@ for human readability. and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008]) - `LobattoLegendreBasis` and related datastructures made fully floating-type general, enabling calculations with higher than double (`Float64`) precision ([#2128]) +- New time integrator `PairedExplicitRK4`, implementing the fourth-order paired explicit Runge-Kutta + method with [Convex.jl](https://github.com/jump-dev/Convex.jl), [ECOS.jl](https://github.com/jump-dev/ECOS.jl) ([#2147]) ## Changes when updating to v0.9 from v0.8.x From 438d9d299022eab6fbdd777fbe5ed9faea987334 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 28 Dec 2024 11:20:25 +0100 Subject: [PATCH 35/61] comments --- .../elixir_euler_vortex_perk4.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index af0b45cbcea..95ac9b045ae 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -20,25 +20,25 @@ The classical isentropic vortex test case as presented in Section 5.1 of """ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D) # Evaluate error after full domain traversion - if t == T_end + if t == t_end t = 0 end - # initial center of the vortex + # Initial center of the vortex inicenter = SVector(0.0, 0.0) - # strength of the vortex + # Strength of the vortex S = 13.5 # Radius of vortex R = 1.5 # Free-stream Mach M = 0.4 - # base flow + # Base flow v1 = 1.0 v2 = 1.0 vel = SVector(v1, v2) - center = inicenter + vel * t # advection of center - center = x - center # distance to centerpoint + center = inicenter + vel * t # Advection of center + center = x - center # Distance to centerpoint center = SVector(center[2], -center[1]) r2 = center[1]^2 + center[2]^2 @@ -46,11 +46,12 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE rho = (1 - (S * M / pi)^2 * (gamma - 1) * exp(2 * f) / 8)^(1 / (gamma - 1)) - du = S / (2 * π * R) * exp(f) # vel. perturbation + du = S / (2 * π * R) * exp(f) # Vel. perturbation vel = vel + du * center v1, v2 = vel p = rho^gamma / (gamma * M^2) + prim = SVector(rho, v1, v2, p) return prim2cons(prim, equations) end @@ -59,8 +60,8 @@ initial_condition = initial_condition_isentropic_vortex edge_length = 20.0 N_passes = 1 -T_end = edge_length * N_passes -tspan = (0.0, T_end) +t_end = edge_length * N_passes +tspan = (0.0, t_end) coordinates_min = (-edge_length / 2, -edge_length / 2) coordinates_max = (edge_length / 2, edge_length / 2) From e106776875547033ba0c220f946c363193bc0392 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 28 Dec 2024 11:26:29 +0100 Subject: [PATCH 36/61] comment --- src/time_integration/methods_2N.jl | 2 +- src/time_integration/methods_3Sstar.jl | 2 +- src/time_integration/methods_SSP.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index fb2847e1c91..5abb8331c7d 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -89,7 +89,7 @@ mutable struct SimpleIntegrator2N{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi.jl sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg opts::SimpleIntegrator2NOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index f197ef60902..1d405cffbe0 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/src/time_integration/methods_3Sstar.jl @@ -156,7 +156,7 @@ mutable struct SimpleIntegrator3Sstar{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time step (iteration) p::Params # will be the semidiscretization from Trixi.jl sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg opts::SimpleIntegrator3SstarOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 33d1f164138..c8cd08e244e 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -89,7 +89,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg opts::SimpleIntegratorSSPOptions finalstep::Bool # added for convenience 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 081373749eb..bc7bd2a13b7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -198,7 +198,7 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience 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 d4371b61896..f52cf6ea2fc 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -191,7 +191,7 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 920eb2d35d1..5ce0cf93992 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -187,7 +187,7 @@ mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F + f::F # `rhs` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience From cdd6cb878840a4f1e8c5ccc7f35734e250a27191 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 2 Jan 2025 13:43:27 +0100 Subject: [PATCH 37/61] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 5ce0cf93992..4e7b6a44aaa 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -118,7 +118,7 @@ end The following structures and methods provide an implementation of the fourth-order paired explicit Runge-Kutta (P-ERK) method -optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). +optimized for a certain simulation setup (PDE, IC & BCs, Riemann Solver, DG Solver). The method has been proposed in - D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024). Fourth-Order Paired-Explicit Runge-Kutta Methods From b9415db7bd40cac84a95e478821232bb45e4f9c9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 2 Jan 2025 13:51:00 +0100 Subject: [PATCH 38/61] plural --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- 2 files changed, 2 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 bc7bd2a13b7..5cd820350f5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -122,7 +122,7 @@ end 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). +optimized for a certain simulation setup (PDE, IC & BCs, Riemann Solver, DG Solver). - Brian Vermeire (2019). Paired explicit Runge-Kutta schemes for stiff systems of equations 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 f52cf6ea2fc..dd41653da99 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -120,7 +120,7 @@ end The following structures and methods provide an implementation of the third-order paired explicit Runge-Kutta (PERK) method -optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). +optimized for a certain simulation setup (PDE, IC & BCs, Riemann Solver, DG Solver). The original paper is - Nasab, Vermeire (2022) Third-order Paired Explicit Runge-Kutta schemes for stiff systems of equations From 47c6ea65b934110cc26d8d1ff240e5f5ef5ebeeb Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 2 Jan 2025 13:53:09 +0100 Subject: [PATCH 39/61] explain `tdir` --- src/time_integration/methods_SSP.jl | 4 ++-- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index c8cd08e244e..e8703aceefe 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -83,14 +83,14 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, du::uType r0::uType t::RealT - tdir::RealT + tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked f::F # `rhs` of the semidiscretization - alg::Alg + alg::Alg # SimpleSSPRK33 opts::SimpleIntegratorSSPOptions finalstep::Bool # added for convenience dtchangeable::Bool 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 5cd820350f5..3634056e493 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -192,7 +192,7 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT + tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) 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 dd41653da99..00d34c5838a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -185,7 +185,7 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT + tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 4e7b6a44aaa..1b3da9123cc 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -181,7 +181,7 @@ mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT + tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) From 2ea94877df4e240758c18b9e4bd12d3f0db46a68 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Jan 2025 13:56:08 +0100 Subject: [PATCH 40/61] remove comment --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 1b3da9123cc..8ff6b48ce3f 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -135,7 +135,7 @@ struct PairedExplicitRK4 <: AbstractPairedExplicitRKSingle c::Vector{Float64} dt_opt::Union{Float64, Nothing} -end # struct PairedExplicitRK4 +end # Constructor for previously computed A Coeffs function PairedExplicitRK4(num_stages, base_path_a_coeffs::AbstractString, From dff1ffeb4df411d04d4f9397724556464a09cc07 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Jan 2025 13:58:30 +0100 Subject: [PATCH 41/61] additional if clause --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 8ff6b48ce3f..7ff3bc5c487 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -80,13 +80,15 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, "a_" * string(num_stages) * ".txt") @assert isfile(path_a_coeffs) "Couldn't find file $path_a_coeffs" - a_coeffs = readdlm(path_a_coeffs, Float64) - num_a_coeffs = size(a_coeffs, 1) + if num_stages > 5 + a_coeffs = readdlm(path_a_coeffs, Float64) + num_a_coeffs = size(a_coeffs, 1) - @assert num_a_coeffs == num_coeffs_max - if num_coeffs_max > 0 - a_matrix[1, :] -= a_coeffs - a_matrix[2, :] = a_coeffs + @assert num_a_coeffs == num_coeffs_max + if num_coeffs_max > 0 + a_matrix[1, :] -= a_coeffs + a_matrix[2, :] = a_coeffs + end end a_matrix_constant = PERK4_a_matrix_constant(cS3) From 3959cb82797d80b66cf6bd1d1ae2018b6156d1da Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 7 Jan 2025 14:23:08 +0100 Subject: [PATCH 42/61] use broadcasting for copy operation --- .../paired_explicit_runge_kutta/methods_PERK3.jl | 4 +--- .../paired_explicit_runge_kutta/methods_PERK4.jl | 4 +--- 2 files changed, 2 insertions(+), 6 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 00d34c5838a..ee17d6861d1 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -268,9 +268,7 @@ function step!(integrator::PairedExplicitRK3Integrator) end # We need to store `du` of the S-1 stage in `kS1` for the final update: - @threaded for i in eachindex(integrator.u) - integrator.kS1[i] = integrator.du[i] - end + integrator.kS1[i] .= integrator.du[i] PERK_ki!(integrator, prob.p, alg, alg.num_stages) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 7ff3bc5c487..1f67217231a 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -259,9 +259,7 @@ end end # Store K_{S-1} in `k1`: - @threaded for i in eachindex(integrator.u) - integrator.k1[i] = integrator.du[i] - end + integrator.k1[i] .= integrator.du[i] integrator.f(integrator.du, integrator.u_tmp, p, integrator.t + alg.c[alg.num_stages] * integrator.dt) From d5eea188df81499136185094ea48661e3674ec98 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 7 Jan 2025 14:28:30 +0100 Subject: [PATCH 43/61] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl --- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 ee17d6861d1..849dae87307 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -268,7 +268,7 @@ function step!(integrator::PairedExplicitRK3Integrator) end # We need to store `du` of the S-1 stage in `kS1` for the final update: - integrator.kS1[i] .= integrator.du[i] + integrator.kS1 .= integrator.du PERK_ki!(integrator, prob.p, alg, alg.num_stages) From 1ad0d61d8f66e5585f16c875715f94a0b26a6117 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 7 Jan 2025 14:29:08 +0100 Subject: [PATCH 44/61] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl --- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 1f67217231a..e45ccd9f252 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -259,7 +259,7 @@ end end # Store K_{S-1} in `k1`: - integrator.k1[i] .= integrator.du[i] + integrator.k1 .= integrator.du integrator.f(integrator.du, integrator.u_tmp, p, integrator.t + alg.c[alg.num_stages] * integrator.dt) From 985dc5f597bcdce42bd26ba88de5f0aa6550c277 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 8 Jan 2025 11:01:51 +0100 Subject: [PATCH 45/61] fix merge errors --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 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 6d6ada9f4d9..13098ba465b 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -33,7 +33,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, # - 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(2, coeffs_max) a_matrix[1, :] = c[3:end] @@ -50,14 +50,14 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, dteps, eig_vals; verbose) - if num_stages != consistency_order + if coeffs_max > 0 monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - a_matrix[1, :] -= A - a_matrix[2, :] = A + a_matrix[:, 1] -= A + a_matrix[:, 2] = A end return a_matrix, c, dt_opt @@ -82,8 +82,6 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, a_matrix[1, :] = c[3:end] if coeffs_max > 0 - a_matrix[:, 1] = c[3:end] - path_monomial_coeffs = joinpath(base_path_monomial_coeffs, "gamma_" * string(num_stages) * ".txt") From 1c787dfc1ba81db55df08b455a9e8010fc2a038d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 8 Jan 2025 11:03:08 +0100 Subject: [PATCH 46/61] fix errors --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 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 13098ba465b..01dc9428f57 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -56,8 +56,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + a_matrix[1, :] -= A + a_matrix[2, :] = A end return a_matrix, c, dt_opt @@ -92,8 +92,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + a_matrix[1, :] -= A + a_matrix[2, :] = A end return a_matrix, c From 01471d5d97541f9d4a20418eea613b0356c66d88 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 8 Jan 2025 14:52:22 +0100 Subject: [PATCH 47/61] typo --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 6 +++--- 1 file changed, 3 insertions(+), 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 070ee904be1..e63c0cfbb3d 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -55,7 +55,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, dteps, eig_vals; verbose) - if coeffs_max > 0 + if num_coeffs_max > 0 monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) num_monomial_coeffs = length(monomial_coeffs) @@ -82,7 +82,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, a_matrix = zeros(2, num_coeffs_max) a_matrix[1, :] = c[3:end] - if coeffs_max > 0 + if num_coeffs_max > 0 path_monomial_coeffs = joinpath(base_path_monomial_coeffs, "gamma_" * string(num_stages) * ".txt") @@ -90,7 +90,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, monomial_coeffs = readdlm(path_monomial_coeffs, Float64) num_monomial_coeffs = size(monomial_coeffs, 1) - @assert num_monomial_coeffs == coeffs_max + @assert num_monomial_coeffs == num_coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) a_matrix[1, :] -= A From 3bec564098a579672ccc1c5335e6be23ced69ac7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 8 Jan 2025 16:30:55 +0100 Subject: [PATCH 48/61] fix errors --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 -- 1 file changed, 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 e63c0cfbb3d..c08f8e64259 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -56,8 +56,6 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, eig_vals; verbose) if num_coeffs_max > 0 - monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, - num_stages) num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == num_coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) From 88000ead42b759ec29f07e4a89a0919fb4335480 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 9 Jan 2025 10:40:25 +0100 Subject: [PATCH 49/61] completeness --- src/time_integration/methods_2N.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index e219c6607b4..a7536de644b 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -57,6 +57,8 @@ the low-storage explicit Runge-Kutta method of Third-order 2N-storage Runge-Kutta schemes with error control URL: https://ntrs.nasa.gov/citations/19940028444 File: https://ntrs.nasa.gov/api/citations/19940028444/downloads/19940028444.pdf + +using the same interface as OrdinaryDiffEq.jl. """ struct CarpenterKennedy2N43 <: SimpleAlgorithm2N a::SVector{4, Float64} From 22b9a5021680295f01ee76ccabe730681dce1b86 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 9 Jan 2025 11:18:05 +0100 Subject: [PATCH 50/61] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- NEWS.md | 8 -------- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 3 files changed, 2 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index bbcad47ba41..aeb6898fea5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,14 +4,6 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -## Changes when updating to v0.10 from v0.9.x - -#### Changed - -- The `PairedExplicitRK2` constructor with second argument `base_path_monomial_coeffs::AbstractString` requires - now `dt_opt`, `bS`, `cS` to be given as keyword arguments ([#2184]). - Previously, those where standard function parameters, in the same order as listed above. - ## Changes in the v0.9 lifecycle #### Added diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 95ac9b045ae..358e1716abd 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -106,7 +106,7 @@ ode_algorithm = Trixi.PairedExplicitRK4(num_stages, path_coeff_file) # run the simulation sol = Trixi.solve(ode, ode_algorithm, - dt = 42.0, # Not used + dt = 42.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); summary_callback() diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index e45ccd9f252..6c006afc07e 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -113,7 +113,7 @@ end In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used. - `tspan`: Time span of the simulation. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the + - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the `ODEProblem` after the equation has been semidiscretized. - `cS3` (`Float64`, optional): Value of $c_{S-3}$ in the Butcher tableau, where $S$ is the number of stages. Default is `1.0f0`. From 164794c5d9cdb807dd3a32c0d1d81a8cf0ebfe23 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 9 Jan 2025 11:22:11 +0100 Subject: [PATCH 51/61] make t_end function --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 358e1716abd..284798da744 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -20,7 +20,7 @@ The classical isentropic vortex test case as presented in Section 5.1 of """ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D) # Evaluate error after full domain traversion - if t == t_end + if t == t_end() t = 0 end @@ -60,8 +60,8 @@ initial_condition = initial_condition_isentropic_vortex edge_length = 20.0 N_passes = 1 -t_end = edge_length * N_passes -tspan = (0.0, t_end) +t_end() = edge_length * N_passes +tspan = (0.0, t_end()) coordinates_min = (-edge_length / 2, -edge_length / 2) coordinates_max = (edge_length / 2, edge_length / 2) From 68cec101ee4d06d11edd14d4fb95e72d02551304 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 9 Jan 2025 11:25:06 +0100 Subject: [PATCH 52/61] avoid globals --- .../structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 284798da744..896f6c4f281 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -57,14 +57,14 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE end initial_condition = initial_condition_isentropic_vortex -edge_length = 20.0 +edge_length() = 20.0 -N_passes = 1 -t_end() = edge_length * N_passes +N_passes() = 1 +t_end() = edge_length() * N_passes() tspan = (0.0, t_end()) -coordinates_min = (-edge_length / 2, -edge_length / 2) -coordinates_max = (edge_length / 2, edge_length / 2) +coordinates_min = (-edge_length() / 2, -edge_length() / 2) +coordinates_max = (edge_length() / 2, edge_length() / 2) cells_per_dimension = (32, 32) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) From 926a097eb4bfa9fbc827ac828c39f39b2f12c474 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 9 Jan 2025 16:54:32 +0100 Subject: [PATCH 53/61] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 5 +++-- examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl | 3 ++- src/time_integration/methods_2N.jl | 2 +- src/time_integration/methods_3Sstar.jl | 2 +- src/time_integration/methods_SSP.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 8 files changed, 11 insertions(+), 9 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 896f6c4f281..7e0098bf318 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -1,5 +1,6 @@ -using OrdinaryDiffEq # Required for `CallbackSet` +# We use time integration methods implemented in Trixi.jl, but we need the `CallbackSet` +using OrdinaryDiffEq: CallbackSet using Trixi # Ratio of specific heats @@ -21,7 +22,7 @@ The classical isentropic vortex test case as presented in Section 5.1 of function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D) # Evaluate error after full domain traversion if t == t_end() - t = 0 + t = zero(t) end # Initial center of the vortex diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl index b43ab3c416d..ac249d23c78 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl @@ -1,5 +1,6 @@ -using OrdinaryDiffEq # Required for `CallbackSet` +# We use time integration methods implemented in Trixi.jl, but we need the `CallbackSet` +using OrdinaryDiffEq: CallbackSet using Trixi # Convex and ECOS are imported because they are used for finding the optimal time step and optimal diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index a7536de644b..65a8aa3e26d 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -103,7 +103,7 @@ mutable struct SimpleIntegrator2N{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi.jl sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg opts::SimpleIntegrator2NOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index 460b0a105e2..c2ef37d57b6 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/src/time_integration/methods_3Sstar.jl @@ -156,7 +156,7 @@ mutable struct SimpleIntegrator3Sstar{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time step (iteration) p::Params # will be the semidiscretization from Trixi.jl sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg opts::SimpleIntegrator3SstarOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index e8703aceefe..071bf3e6f50 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -89,7 +89,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg # SimpleSSPRK33 opts::SimpleIntegratorSSPOptions finalstep::Bool # added for convenience 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 c08f8e64259..6464271f1af 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -200,7 +200,7 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience 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 849dae87307..d8fac2d8e10 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -191,7 +191,7 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 6c006afc07e..2c613361756 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -189,7 +189,7 @@ mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F, iter::Int # current number of time steps (iteration) p::Params # will be the semidiscretization from Trixi sol::Sol # faked - f::F # `rhs` of the semidiscretization + f::F # `rhs!` of the semidiscretization alg::Alg # This is our own class written above; Abbreviation for ALGorithm opts::PairedExplicitRKOptions finalstep::Bool # added for convenience From 5bba4cddc1d361245d87de6649cdece4f356c6ff Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 9 Jan 2025 16:55:51 +0100 Subject: [PATCH 54/61] i.e., --- src/time_integration/methods_SSP.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- .../paired_explicit_runge_kutta/methods_PERK4.jl | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 071bf3e6f50..91bacb4a7cf 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -83,7 +83,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, du::uType r0::uType t::RealT - tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time + tdir::RealT # DIRection of time integration, i.e., if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) 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 6464271f1af..c219baeef78 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -194,7 +194,7 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time + tdir::RealT # DIRection of time integration, i.e., if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) 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 d8fac2d8e10..5dcfe63d9a7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -185,7 +185,7 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time + tdir::RealT # DIRection of time integration, i.e., if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 2c613361756..3c6f519dd16 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -183,7 +183,7 @@ mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F, du::uType u_tmp::uType t::RealT - tdir::RealT # DIRection of time integration, i.e, if one marches forward or backward in time + tdir::RealT # DIRection of time integration, i.e., if one marches forward or backward in time dt::RealT # current time step dtcache::RealT # manually set time step iter::Int # current number of time steps (iteration) From 366b84fb56f09fc720430bb32a7a2057a807ead2 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 10:30:01 +0100 Subject: [PATCH 55/61] shorten opt ext --- ext/TrixiConvexECOSExt.jl | 192 ++++++------------ .../methods_PERK2.jl | 3 +- .../methods_PERK3.jl | 3 +- .../methods_PERK4.jl | 13 +- .../polynomial_optimizer.jl | 1 - test/test_tree_1d_hypdiff.jl | 4 +- test/test_unit.jl | 18 +- 7 files changed, 80 insertions(+), 154 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 018c7a7b5cb..52d83d867de 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -10,7 +10,7 @@ using LinearAlgebra: eigvals # Use functions that are to be extended and additional symbols that are not exported using Trixi: Trixi, undo_normalization!, - bisect_stability_polynomial, bisect_stability_polynomial_PERK4, @muladd + bisect_stability_polynomial, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). # Since these FMAs can increase the performance of many numerical algorithms, @@ -31,20 +31,17 @@ end # Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency # order, including contributions from free coefficients for higher orders, and # return the maximum absolute value -function stability_polynomials!(pnoms, consistency_order, - num_stage_evals, +function stability_polynomials!(pnoms, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, - gamma) + gamma, consistency_order) # Initialize with zero'th order (z^0) coefficient for i in 1:num_eig_vals pnoms[i] = 1.0 end # First `consistency_order` terms of the exponential for k in 1:consistency_order - for i in 1:num_eig_vals - pnoms[i] += normalized_powered_eigvals_scaled[i, k] - end + pnoms += view(normalized_powered_eigvals_scaled, :, k) end # Contribution from free coefficients @@ -64,13 +61,12 @@ end # Specialized form of the stability polynomials for fourth-order PERK schemes. function stability_polynomials_PERK4!(pnoms, num_stage_evals, num_eig_vals, - normalized_powered_eigvals, - gamma, - dt, cS3) + normalized_powered_eigvals_scaled, + gamma, cS3) # Constants arising from the particular form of Butcher tableau chosen for the 4th order PERK methods + # cS3 = c_{S-3} k1 = 0.001055026310046423 / cS3 k2 = 0.03726406530405851 / cS3 - # Note: `cS3` = c_{S-3} is in principle free, while the other abscissae are fixed to 1.0 # Initialize with zero'th order (z^0) coefficient for i in 1:num_eig_vals @@ -78,22 +74,19 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, end # First `consistency_order` = 4 terms of the exponential for k in 1:4 - for i in 1:num_eig_vals - pnoms[i] += dt^k * normalized_powered_eigvals[i, k] - end + pnoms += view(normalized_powered_eigvals_scaled, :, k) end # "Fixed" term due to choice of the PERK4 Butcher tableau # Required to un-do the normalization of the eigenvalues here - pnoms += k1 * dt^5 * view(normalized_powered_eigvals, :, 5) * factorial(5) + pnoms += k1 * view(normalized_powered_eigvals_scaled, :, 5) * factorial(5) # Contribution from free coefficients for k in 1:(num_stage_evals - 5) - pnoms += (k2 * dt^(k + 4) * view(normalized_powered_eigvals, :, k + 4) * + pnoms += (k2 * view(normalized_powered_eigvals_scaled, :, k + 4) * gamma[k] + - k1 * dt^(k + 5) * view(normalized_powered_eigvals, :, k + 5) * - gamma[k] * - (k + 5)) + k1 * view(normalized_powered_eigvals_scaled, :, k + 5) * + gamma[k] * (k + 5)) # Undo normalization of the second summand end # For optimization only the maximum is relevant @@ -124,22 +117,34 @@ 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) + +For the fourth-order PERK schemes, a specialized optimization routine is required, see +- D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024). + Fourth-Order Paired-Explicit Runge-Kutta Methods + [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470) =# # Perform bisection to optimize timestep for stability of the polynomial function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals; - verbose = false) + kwargs...) dtmin = 0.0 dt = -1.0 abs_p = -1.0 + if consistency_order == 4 + # Fourth-order scheme has one additional fixed coefficient + num_reduced_unkown = 5 + else # p = 2, 3 + num_reduced_unkown = consistency_order + end + # Construct stability polynomial for each eigenvalue pnoms = ones(Complex{Float64}, num_eig_vals, 1) # Init datastructure for monomial coefficients - gamma = Variable(num_stage_evals - consistency_order) + gamma = Variable(num_stage_evals - num_reduced_unkown) normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) normalize_power_eigvals!(normalized_powered_eigvals, @@ -148,7 +153,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) - if verbose + if kwargs[:verbose] println("Start optimization of stability polynomial \n") end @@ -156,24 +161,28 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, 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 + @views normalized_powered_eigvals_scaled[:, k] = dt^k .* + normalized_powered_eigvals[:, + k] end # Check if there are variables to optimize - if num_stage_evals - consistency_order > 0 + if num_stage_evals - num_reduced_unkown > 0 # Use last optimal values for gamma in (potentially) next iteration - problem = minimize(stability_polynomials!(pnoms, consistency_order, - num_stage_evals, - num_eig_vals, - normalized_powered_eigvals_scaled, - gamma)) + if consistency_order == 4 + problem = minimize(stability_polynomials_PERK4!(pnoms, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma, kwargs[:cS3])) + else # p = 2, 3 + problem = minimize(stability_polynomials!(pnoms, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma, consistency_order)) + end solve!(problem, # Parameters taken from default values for EiCOS @@ -191,11 +200,19 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, abs_p = problem.optval else - abs_p = stability_polynomials!(pnoms, consistency_order, - num_stage_evals, - num_eig_vals, - normalized_powered_eigvals_scaled, - gamma) + if consistency_order == 4 + abs_p = stability_polynomials_PERK4!(pnoms, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma, kwargs[:cS3]) + else + abs_p = stability_polynomials!(pnoms, + num_stage_evals, + num_eig_vals, + normalized_powered_eigvals_scaled, + gamma, consistency_order) + end end if abs_p < 1 @@ -205,7 +222,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, end end - if verbose + if kwargs[:verbose] println("Concluded stability polynomial optimization \n") end @@ -215,102 +232,13 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. end - # Catch case S = 3 (only one opt. variable) + # Catch case with only one optimization variable if isa(gamma_opt, Number) gamma_opt = [gamma_opt] end undo_normalization!(gamma_opt, num_stage_evals, - consistency_order, consistency_order) - - return gamma_opt, dt -end - -# Specialized routine for PERK4. -# For details, see Section 4 in -# - D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024). -# Fourth-Order Paired-Explicit Runge-Kutta Methods -# [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470) -function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals, - num_stage_evals, - dtmax, dteps, eig_vals, cS3; - verbose = false) - dtmin = 0.0 - dt = -1.0 - abs_p = -1.0 - - # Construct stability polynomial for each eigenvalue - pnoms = ones(Complex{Float64}, num_eig_vals, 1) - - # Init datastructure for monomial coefficients - gamma = Variable(num_stage_evals - 5) - - normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - normalize_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, - num_stage_evals) - - if verbose - println("Start optimization of stability polynomial \n") - end - - # Bisection on timestep - while dtmax - dtmin > dteps - dt = 0.5 * (dtmax + dtmin) - - if num_stage_evals > 5 - # Use last optimal values for gamma in (potentially) next iteration - problem = minimize(stability_polynomials_PERK4!(pnoms, - num_stage_evals, - num_eig_vals, - normalized_powered_eigvals, - gamma, dt, cS3)) - - solve!(problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - "feastol" => 1e-9, - "abstol" => 1e-9, - "reltol" => 1e-9, - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent = true) - - abs_p = problem.optval - else - abs_p = stability_polynomials_PERK4!(pnoms, num_stage_evals, - num_eig_vals, - normalized_powered_eigvals, - gamma, dt, cS3) - end - - if abs_p < 1 - dtmin = dt - else - dtmax = dt - end - end - - if verbose - println("Concluded stability polynomial optimization \n") - end - - if num_stage_evals > 5 - gamma_opt = evaluate(gamma) - else - gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. - end - - # Catch case S = 6 (only one opt. variable) - if isa(gamma_opt, Number) - gamma_opt = [gamma_opt] - end - - undo_normalization!(gamma_opt, num_stage_evals, 5, 4) + num_reduced_unkown, consistency_order) return gamma_opt, dt end 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 c219baeef78..2263fbb0712 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -42,13 +42,12 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, a_matrix = zeros(2, num_coeffs_max) a_matrix[1, :] = c[3:end] - consistency_order = 2 - dtmax = tspan[2] - tspan[1] dteps = 1e-9 # Hyperparameter of the optimization, might be too large for systems requiring very small timesteps num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + consistency_order = 2 monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, num_eig_vals, num_stages, dtmax, 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 5dcfe63d9a7..35b0a772919 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -33,13 +33,12 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, # Initialize the array of our solution a_unknown = zeros(num_stages - 2) - # Calculate coefficients of the stability polynomial in monomial form - consistency_order = 3 dtmax = tspan[2] - tspan[1] dteps = 1.0f-9 num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + consistency_order = 3 monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, num_eig_vals, num_stages, dtmax, dteps, diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index 3c6f519dd16..6c097a0e2ee 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -33,17 +33,18 @@ function compute_PairedExplicitRK4_butcher_tableau(num_stages, tspan, num_coeffs_max = num_stages - 5 a_matrix = zeros(2, num_coeffs_max) - # Calculate coefficients of the stability polynomial in monomial form dtmax = tspan[2] - tspan[1] dteps = 1.0f-9 num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) - monomial_coeffs, dt_opt = bisect_stability_polynomial_PERK4(num_eig_vals, - num_stages, - dtmax, dteps, - eig_vals, cS3; - verbose) + consistency_order = 4 + monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, + num_eig_vals, + num_stages, + dtmax, dteps, + eig_vals; + verbose, cS3) if num_stages > 5 a_unknown = copy(monomial_coeffs) 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 6658f1ac4e6..bfd53ba2eaf 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -26,4 +26,3 @@ end # extension or by the Convex and ECOS-specific code loaded by Requires.jl function undo_normalization! end function bisect_stability_polynomial end -function bisect_stability_polynomial_PERK4 end diff --git a/test/test_tree_1d_hypdiff.jl b/test/test_tree_1d_hypdiff.jl index b2515d20a01..7e1166a1d1e 100644 --- a/test/test_tree_1d_hypdiff.jl +++ b/test/test_tree_1d_hypdiff.jl @@ -27,8 +27,8 @@ end @trixi_testset "elixir_hypdiff_nonperiodic_perk4.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic_perk4.jl"), - l2=[1.3655114989569954e-7, 1.020034502220849e-6], - linf=[7.173289721662535e-7, 4.507115265006689e-6], + l2=[1.3655114994521285e-7, 1.0200345014751413e-6], + linf=[7.173289867656862e-7, 4.507115296537023e-6], atol=2.5e-13) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index 48f53572516..13654a19534 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1778,15 +1778,15 @@ end ode_algorithm = Trixi.PairedExplicitRK4(14, tspan, vec(eig_vals)) @test isapprox(transpose(ode_algorithm.a_matrix), - [0.9935760056654522 0.006423994334547779 - 0.984991598524171 0.01500840147582901 - 0.9731962964227893 0.026803703577210732 - 0.9564649429772978 0.04353505702270225 - 0.9319644395990909 0.0680355604009091 - 0.8955295433261026 0.10447045667389748 - 0.8444449101874965 0.15555508981250354 - 0.7923618582881918 0.20763814171180817 - 0.7723161199637355 0.2276838800362645], atol = 1e-13) + [0.9935765040401348 0.0064234959598652 + 0.9849926812139576 0.0150073187860425 + 0.9731978940975923 0.0268021059024077 + 0.9564664284695985 0.0435335715304015 + 0.9319632992510594 0.0680367007489407 + 0.8955171743167522 0.1044828256832478 + 0.8443975130657495 0.1556024869342504 + 0.7922561745278265 0.2077438254721735 + 0.7722324105428290 0.2277675894571710], atol = 1e-13) end @testset "Sutherlands Law" begin From 40c4f097261d892409e32b92ed04faa0a01d6ea1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 10:32:16 +0100 Subject: [PATCH 56/61] typos --- ext/TrixiConvexECOSExt.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 52d83d867de..f3bfff73449 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -86,7 +86,7 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, pnoms += (k2 * view(normalized_powered_eigvals_scaled, :, k + 4) * gamma[k] + k1 * view(normalized_powered_eigvals_scaled, :, k + 5) * - gamma[k] * (k + 5)) # Undo normalization of the second summand + gamma[k] * (k + 5)) # Ensure same normalization of both summands end # For optimization only the maximum is relevant @@ -135,16 +135,16 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, if consistency_order == 4 # Fourth-order scheme has one additional fixed coefficient - num_reduced_unkown = 5 + num_reduced_unknown = 5 else # p = 2, 3 - num_reduced_unkown = consistency_order + num_reduced_unknown = consistency_order end # Construct stability polynomial for each eigenvalue pnoms = ones(Complex{Float64}, num_eig_vals, 1) # Init datastructure for monomial coefficients - gamma = Variable(num_stage_evals - num_reduced_unkown) + gamma = Variable(num_stage_evals - num_reduced_unknown) normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) normalize_power_eigvals!(normalized_powered_eigvals, @@ -168,7 +168,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, end # Check if there are variables to optimize - if num_stage_evals - num_reduced_unkown > 0 + if num_stage_evals - num_reduced_unknown > 0 # Use last optimal values for gamma in (potentially) next iteration if consistency_order == 4 problem = minimize(stability_polynomials_PERK4!(pnoms, @@ -238,7 +238,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, end undo_normalization!(gamma_opt, num_stage_evals, - num_reduced_unkown, consistency_order) + num_reduced_unknown, consistency_order) return gamma_opt, dt end From 4195f72d2b1961545c51560d1589742461c58635 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 10:35:06 +0100 Subject: [PATCH 57/61] revert --- ext/TrixiConvexECOSExt.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index f3bfff73449..1d45e0bdd0f 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -9,8 +9,7 @@ using ECOS: Optimizer using LinearAlgebra: eigvals # Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, undo_normalization!, - bisect_stability_polynomial, @muladd +using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). # Since these FMAs can increase the performance of many numerical algorithms, From 99f754b7c226e1dad59c00a0485c78c9ca9b170b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 10:42:53 +0100 Subject: [PATCH 58/61] shorten --- ext/TrixiConvexECOSExt.jl | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 1d45e0bdd0f..e9fbd53fec8 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -31,13 +31,11 @@ end # order, including contributions from free coefficients for higher orders, and # return the maximum absolute value function stability_polynomials!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, consistency_order) # Initialize with zero'th order (z^0) coefficient - for i in 1:num_eig_vals - pnoms[i] = 1.0 - end + pnoms .= 1 + # First `consistency_order` terms of the exponential for k in 1:consistency_order pnoms += view(normalized_powered_eigvals_scaled, :, k) @@ -59,7 +57,6 @@ end # Specialized form of the stability polynomials for fourth-order PERK schemes. function stability_polynomials_PERK4!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, cS3) # Constants arising from the particular form of Butcher tableau chosen for the 4th order PERK methods @@ -68,9 +65,8 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, k2 = 0.03726406530405851 / cS3 # Initialize with zero'th order (z^0) coefficient - for i in 1:num_eig_vals - pnoms[i] = 1.0 - end + pnoms .= 1 + # First `consistency_order` = 4 terms of the exponential for k in 1:4 pnoms += view(normalized_powered_eigvals_scaled, :, k) @@ -97,13 +93,10 @@ function stability_polynomials_PERK4!(pnoms, num_stage_evals, end @inline function normalize_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, + 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 + @views normalized_powered_eigvals[:, j] = eig_vals[:] .^ j ./ factorial(j) end end @@ -147,7 +140,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) normalize_power_eigvals!(normalized_powered_eigvals, - num_eig_vals, eig_vals, + eig_vals, num_stage_evals) normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) @@ -172,13 +165,11 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, if consistency_order == 4 problem = minimize(stability_polynomials_PERK4!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, kwargs[:cS3])) else # p = 2, 3 problem = minimize(stability_polynomials!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, consistency_order)) end @@ -202,13 +193,11 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, if consistency_order == 4 abs_p = stability_polynomials_PERK4!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, kwargs[:cS3]) else abs_p = stability_polynomials!(pnoms, num_stage_evals, - num_eig_vals, normalized_powered_eigvals_scaled, gamma, consistency_order) end From c416a8112b342bda62fcf9c5f53fec93574632c4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 10:46:41 +0100 Subject: [PATCH 59/61] do not export --- ext/TrixiConvexECOSExt.jl | 6 +++--- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index e9fbd53fec8..d9e4aeedec5 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -9,7 +9,7 @@ using ECOS: Optimizer using LinearAlgebra: eigvals # Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd +using Trixi: Trixi, bisect_stability_polynomial, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). # Since these FMAs can increase the performance of many numerical algorithms, @@ -20,8 +20,8 @@ using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd # Undo normalization of stability polynomial coefficients by index factorial # relative to consistency order. -function Trixi.undo_normalization!(gamma_opt, num_stage_evals, - num_reduced_coeffs, fac_offset) +function undo_normalization!(gamma_opt, num_stage_evals, + num_reduced_coeffs, fac_offset) for k in 1:(num_stage_evals - num_reduced_coeffs) gamma_opt[k] /= factorial(k + fac_offset) end 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 bfd53ba2eaf..0598427ad05 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -24,5 +24,4 @@ end # Add definitions of functions related to polynomial optimization by Convex and ECOS here # such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package # extension or by the Convex and ECOS-specific code loaded by Requires.jl -function undo_normalization! end function bisect_stability_polynomial end From 0ba9120748d60783196bbb7251ff07e1ca951c5c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 11:30:37 +0100 Subject: [PATCH 60/61] convert float --- .../elixir_euler_vortex_perk4.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 7e0098bf318..9ab40c9ad6b 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -25,17 +25,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE t = zero(t) end + RealT = eltype(x) # Initial center of the vortex - inicenter = SVector(0.0, 0.0) + inicenter = SVector(0, 0) # Strength of the vortex - S = 13.5 + S = convert(RealT, 13.5) # Radius of vortex - R = 1.5 + R = convert(RealT, 1.5) # Free-stream Mach - M = 0.4 + M = convert(RealT, 0.4) # Base flow - v1 = 1.0 - v2 = 1.0 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) center = inicenter + vel * t # Advection of center @@ -45,9 +46,10 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE f = (1 - r2) / (2 * R^2) - rho = (1 - (S * M / pi)^2 * (gamma - 1) * exp(2 * f) / 8)^(1 / (gamma - 1)) + rho = (1 - (S * M / convert(RealT, pi))^2 * (gamma - 1) * exp(2 * f) / 8)^(1 / + (gamma - 1)) - du = S / (2 * π * R) * exp(f) # Vel. perturbation + du = S / (2 * convert(RealT, pi) * R) * exp(f) # Vel. perturbation vel = vel + du * center v1, v2 = vel From cb23fe6a6f94fdf692222c748998b4d57d402bcc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 10 Jan 2025 17:36:07 +0100 Subject: [PATCH 61/61] remove first blank line --- examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl | 1 - examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl index 9ab40c9ad6b..e35dec5854a 100644 --- a/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl +++ b/examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl @@ -1,4 +1,3 @@ - # We use time integration methods implemented in Trixi.jl, but we need the `CallbackSet` using OrdinaryDiffEq: CallbackSet using Trixi diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl index ac249d23c78..55c1d38e3a6 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl @@ -1,4 +1,3 @@ - # We use time integration methods implemented in Trixi.jl, but we need the `CallbackSet` using OrdinaryDiffEq: CallbackSet using Trixi