Skip to content

Commit

Permalink
shorten opt ext
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed Jan 10, 2025
1 parent 5bba4cd commit 366b84f
Show file tree
Hide file tree
Showing 7 changed files with 80 additions and 154 deletions.
192 changes: 60 additions & 132 deletions ext/TrixiConvexECOSExt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -64,36 +61,32 @@ 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
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
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
Expand Down Expand Up @@ -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

Check warning on line 138 in ext/TrixiConvexECOSExt.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"unkown" should be "unknown".
else # p = 2, 3
num_reduced_unkown = consistency_order

Check warning on line 140 in ext/TrixiConvexECOSExt.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"unkown" should be "unknown".
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)

Check warning on line 147 in ext/TrixiConvexECOSExt.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"unkown" should be "unknown".

normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals)
normalize_power_eigvals!(normalized_powered_eigvals,
Expand All @@ -148,32 +153,36 @@ 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

# Bisection on timestep
while dtmax - dtmin > dteps
dt = 0.5 * (dtmax + dtmin)

# Compute stability polynomial for current timestep
for k in 1:num_stage_evals
dt_k = dt^k
for i in 1:num_eig_vals
normalized_powered_eigvals_scaled[i, k] = dt_k *
normalized_powered_eigvals[i,
k]
end
@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

Check warning on line 171 in ext/TrixiConvexECOSExt.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"unkown" should be "unknown".
# 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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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)

Check warning on line 241 in ext/TrixiConvexECOSExt.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"unkown" should be "unknown".

return gamma_opt, dt
end
Expand Down
3 changes: 1 addition & 2 deletions src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
3 changes: 1 addition & 2 deletions src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
13 changes: 7 additions & 6 deletions src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
1 change: 0 additions & 1 deletion src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions test/test_tree_1d_hypdiff.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 366b84f

Please sign in to comment.