From 15cb528f836ede520a1ac3c2e995c95861c377ba Mon Sep 17 00:00:00 2001 From: Warisa Date: Sun, 24 Nov 2024 12:41:04 +0100 Subject: [PATCH 1/2] handle s=consistency_order --- ext/TrixiConvexECOSExt.jl | 12 ++++++++-- .../methods_PERK2.jl | 17 +++++++------- .../methods_PERK3.jl | 23 ++++++++++--------- 3 files changed, 31 insertions(+), 21 deletions(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 8251fe3eed9..9b897436dbb 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -60,7 +60,11 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals, end # For optimization only the maximum is relevant - return maximum(abs(pnoms)) + if consistency_order - num_stage_evals == 0 + 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 #= @@ -151,7 +155,11 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, println("Concluded stability polynomial optimization \n") end - gamma_opt = evaluate(gamma) + if consistency_order - num_stage_evals != 0 + gamma_opt = evaluate(gamma) + else + gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. + end # Catch case S = 3 (only one opt. variable) if isa(gamma_opt, Number) 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..6bdaaa44b3c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -49,15 +49,16 @@ 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 - A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + if num_stages != 2 + 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 + end return a_matrix, c, dt_opt 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 477748b28c5..5759fa5db98 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -33,21 +33,22 @@ 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) + + monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, + num_eig_vals, num_stages, + dtmax, dteps, + eig_vals; verbose) + # Special case of e = 3 if num_stages == 3 a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau else - # 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) - - monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, - num_eig_vals, num_stages, - dtmax, dteps, - eig_vals; verbose) monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) From f2431579335aa520f076b15b7304d0f18c25c2b3 Mon Sep 17 00:00:00 2001 From: Warisa Date: Sun, 24 Nov 2024 12:46:12 +0100 Subject: [PATCH 2/2] minor change --- .../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 6bdaaa44b3c..2451680a505 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -50,7 +50,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, dteps, eig_vals; verbose) - if num_stages != 2 + if num_stages != consistency_order monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) num_monomial_coeffs = length(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 5759fa5db98..02bc3eeba36 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -46,7 +46,7 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, eig_vals; verbose) # Special case of e = 3 - if num_stages == 3 + if num_stages == consistency_order a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau else monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,