From 87399dc00081ce39ceae78b7fddd1d2176e54581 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 14 Feb 2024 09:57:40 +0100 Subject: [PATCH 1/4] add constructor of PERK3 --- .../methods_PERK2.jl | 1 - .../methods_PERK3.jl | 121 ++++++++++++++++++ .../paired_explicit_runge_kutta.jl | 1 + .../polynomial_optimizer.jl | 22 ++-- 4 files changed, 133 insertions(+), 12 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 11cbeba7154..f8e061c0872 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,7 +2,6 @@ # Since these FMAs can increase the performance of many numerical algorithms, # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -include("polynomial_optimizer.jl") @muladd begin abstract type PERK end diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 61cefbe8ac0..fd3d1705cd6 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -2,8 +2,117 @@ # Since these FMAs can increase the performance of many numerical algorithms, # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +using NLsolve @muladd begin +function F_own(aUnknown, NumStages, NumStageEvals, MonCoeffs) + cS2 = 1 # cS2 = c_ts(S-2) + c_ts = c_Own(cS2, NumStages) # ts = timestep + + c_eq = zeros(NumStageEvals - 2) # Add equality constraint that cS2 is equal to 1 + # Both terms should be present + for i in 1:(NumStageEvals - 4) + term1 = aUnknown[NumStageEvals - 1] + term2 = aUnknown[NumStageEvals] + for j in 1:i + term1 *= aUnknown[NumStageEvals - 1 - j] + term2 *= aUnknown[NumStageEvals - j] + end + term1 *= c_ts[NumStages - 2 - i] * 1/6 + term2 *= c_ts[NumStages - 1 - i] * 4/6 + + c_eq[i] = MonCoeffs[i] - (term1 + term2) + end + + # Highest coefficient: Only one term present + i = NumStageEvals - 3 + term2 = aUnknown[NumStageEvals] + for j in 1:i + term2 *= aUnknown[NumStageEvals - j] + end + term2 *= c_ts[NumStages - 1 - i] * 4/6 + + c_eq[i] = MonCoeffs[i] - term2 + + c_eq[NumStageEvals - 2] = 1.0 - 4 * aUnknown[NumStageEvals] - aUnknown[NumStageEvals - 1] + + return c_eq +end + +function c_Own(cS2, NumStages) + c_ts = zeros(NumStages); + + # Last timesteps as for SSPRK33 + c_ts[NumStages] = 0.5; + c_ts[NumStages - 1] = 1; + + # Linear increasing timestep for remainder + for i = 2:NumStages-2 + c_ts[i] = cS2 * (i-1)/(NumStages - 3); + end + return c_ts +end + +function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSemidiscretization) + + #Initialize array of c + c_ts = c_Own(1.0, NumStages) + + #Initialize the array of our solution + aUnknown = zeros(NumStageEvals) + + #Case of e = 3 + if NumStageEvals == 3 + aUnknown = [0, c_ts[2], 0.25] + + else + #Calculate MonCoeffs from polynomial Optimizer + ConsOrder = 3 + dtMax = 1.0 + dtEps = 1e-9 + filter_thres = 1e-12 + + #Get EigVals + J = jacobian_ad_forward(semi) + EigVals = eigvals(J) + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + + MonCoeffs, _, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + + #Define the objective_function + function objective_function(x) + return F_own(x, NumStages, NumStageEvals, MonCoeffs) + end + + #call nlsolver to solve until the answer is not NaN or negative values + is_sol_valid = false + while is_sol_valid == false + #Initialize initial guess + x0 = 0.1 .* rand(NumStageEvals) + x0[1] = Base.big(0) + x0[2] = c_ts[2] + + sol = nlsolve(objective_function, method = :trust_region, x0, + ftol = 1e-15, iterations = 10^4, xtol = 1e-13) + + aUnknown = sol.zero + + is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) + end + end + + println("aUnknown") + println(aUnknown[3:end]) #To debug + + AMatrix = zeros(NumStages - 2, 2) + AMatrix[:, 1] = c_ts[3:end] + AMatrix[:, 1] -= aUnknown[3:end] + AMatrix[:, 2] = aUnknown[3:end] + + return AMatrix, c_ts +end + function ComputePERK3_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, cS2::Float64) # c Vector form Butcher Tableau (defines timestep per stage) @@ -67,6 +176,18 @@ mutable struct PERK3 <: PERKSingle ComputePERK3_ButcherTableau(NumStages_, BasePathMonCoeffs_, cS2_) return newPERK3 end + + #Constructor that compute A coefficients from semidiscretization + function PERK3(NumStages_::Int, NumStageEvals_ ::Int, semi_::AbstractSemidiscretization) + + newPERK3 = new(NumStages_) + + newPERK3.AMatrix, newPERK3.c = + ComputePERK3_ButcherTableau(NumStages_, NumStageEvals_, semi_) + + return newPERK3 + end + end # struct PERK3 diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl index 03689b75321..9bb3d0b817c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -7,5 +7,6 @@ include("methods_PERK2.jl") include("methods_PERK3.jl") +include("polynomial_optimizer.jl") end # @muladd \ No newline at end of file diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 7c22114dea7..1ad29cd2f0c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -46,7 +46,7 @@ function ReadInEigVals(Path_toEvalFile::AbstractString) return NumEigVals, EigVals end -function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) +function Polynoms(ConsOrder::Int, NumStageEvals::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) for i in 1:NumEigVals pnoms[i] = 1.0 @@ -58,7 +58,7 @@ function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_po end end - for k in ConsOrder + 1:NumStages + for k in ConsOrder + 1:NumStageEvals pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] end @@ -66,7 +66,7 @@ function Polynoms(ConsOrder::Int, NumStages::Int, NumEigVals::Int, normalized_po end -function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) +function Bisection(ConsOrder::Int, NumEigVals::Int, NumStageEvals::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) dtMin = 0.0 @@ -76,23 +76,23 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float pnoms = ones(Complex{Float64}, NumEigVals, 1) - gamma = Variable(NumStages - ConsOrder) # Init datastructure for results + gamma = Variable(NumStageEvals - ConsOrder) # Init datastructure for results - normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStages) + normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStageEvals) - for j in 1:NumStages + for j in 1:NumStageEvals fac_j = factorial(j) for i in 1:NumEigVals normalized_powered_Eigvals[i, j] = EigVals[i]^j / fac_j end end - normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStages) + normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStageEvals) while dtMax - dtMin > dtEps dt = 0.5 * (dtMax + dtMin) - for k in 1:NumStages + for k in 1:NumStageEvals dt_k = dt^k for i in 1:NumEigVals normalized_powered_EigvalsScaled[i,k] = dt_k * normalized_powered_Eigvals[i,k] @@ -100,7 +100,7 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float end # Use last optimal values for gamm0 in (potentially) next iteration - problem = minimize(Polynoms(ConsOrder, NumStages, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) + problem = minimize(Polynoms(ConsOrder, NumStageEvals, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) Convex.solve!( problem, @@ -133,8 +133,8 @@ function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float return evaluate(gamma), AbsP, dt end -function undo_normalization(ConsOrder::Int, NumStages::Int, gammaOpt) - for k in ConsOrder + 1:NumStages +function undo_normalization(ConsOrder::Int, NumStageEvals::Int, gammaOpt) + for k in ConsOrder + 1:NumStageEvals gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) end return gammaOpt From e518386946efe48d9c7d306e19eca4adcec47566 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 14 Feb 2024 16:42:20 +0100 Subject: [PATCH 2/4] minor fixes --- .../paired_explicit_runge_kutta/methods_PERK3.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index fd3d1705cd6..79f090711a8 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -70,14 +70,14 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem ConsOrder = 3 dtMax = 1.0 dtEps = 1e-9 - filter_thres = 1e-12 + filter_threshold = 1e-12 #Get EigVals J = jacobian_ad_forward(semi) EigVals = eigvals(J) - NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_threshold) - MonCoeffs, _, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs, dt, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) #Define the objective_function @@ -90,7 +90,7 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem while is_sol_valid == false #Initialize initial guess x0 = 0.1 .* rand(NumStageEvals) - x0[1] = Base.big(0) + x0[1] = 0.0 x0[2] = c_ts[2] sol = nlsolve(objective_function, method = :trust_region, x0, @@ -98,6 +98,7 @@ function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSem aUnknown = sol.zero + #First check the values of a[i, i-1] which are acquired by nonlinear solver and c[i] - a[i, i-1]. is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) end end From a1222d7d37bb65e4114eee323e999902807b3df8 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 15 Feb 2024 20:43:22 +0100 Subject: [PATCH 3/4] function and variable name adjustments --- .../methods_PERK3.jl | 205 +++++++++--------- .../polynomial_optimizer.jl | 188 ++++++++-------- 2 files changed, 198 insertions(+), 195 deletions(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 79f090711a8..5901d3b5e65 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -5,150 +5,149 @@ using NLsolve @muladd begin -function F_own(aUnknown, NumStages, NumStageEvals, MonCoeffs) - cS2 = 1 # cS2 = c_ts(S-2) - c_ts = c_Own(cS2, NumStages) # ts = timestep +function f_own(a_unknown, num_stages, num_stage_evals, mon_coeffs) + c_s2 = 1 # c_s2 = c_ts(S-2) + c_ts = c_own(c_s2, num_stages) # ts = timestep - c_eq = zeros(NumStageEvals - 2) # Add equality constraint that cS2 is equal to 1 + c_eq = zeros(num_stage_evals - 2) # Add equality constraint that c_s2 is equal to 1 # Both terms should be present - for i in 1:(NumStageEvals - 4) - term1 = aUnknown[NumStageEvals - 1] - term2 = aUnknown[NumStageEvals] + for i in 1:(num_stage_evals - 4) + term1 = a_unknown[num_stage_evals - 1] + term2 = a_unknown[num_stage_evals] for j in 1:i - term1 *= aUnknown[NumStageEvals - 1 - j] - term2 *= aUnknown[NumStageEvals - j] + term1 *= a_unknown[num_stage_evals - 1 - j] + term2 *= a_unknown[num_stage_evals - j] end - term1 *= c_ts[NumStages - 2 - i] * 1/6 - term2 *= c_ts[NumStages - 1 - i] * 4/6 + term1 *= c_ts[num_stages - 2 - i] * 1/6 + term2 *= c_ts[num_stages - 1 - i] * 4/6 - c_eq[i] = MonCoeffs[i] - (term1 + term2) + c_eq[i] = mon_coeffs[i] - (term1 + term2) end # Highest coefficient: Only one term present - i = NumStageEvals - 3 - term2 = aUnknown[NumStageEvals] + i = num_stage_evals - 3 + term2 = a_unknown[num_stage_evals] for j in 1:i - term2 *= aUnknown[NumStageEvals - j] + term2 *= a_unknown[num_stage_evals - j] end - term2 *= c_ts[NumStages - 1 - i] * 4/6 + term2 *= c_ts[num_stages - 1 - i] * 4/6 - c_eq[i] = MonCoeffs[i] - term2 + c_eq[i] = mon_coeffs[i] - term2 - c_eq[NumStageEvals - 2] = 1.0 - 4 * aUnknown[NumStageEvals] - aUnknown[NumStageEvals - 1] + c_eq[num_stage_evals - 2] = 1.0 - 4 * a_unknown[num_stage_evals] - a_unknown[num_stage_evals - 1] return c_eq end -function c_Own(cS2, NumStages) - c_ts = zeros(NumStages); +function c_own(c_s2, num_stages) + c_ts = zeros(num_stages) # Last timesteps as for SSPRK33 - c_ts[NumStages] = 0.5; - c_ts[NumStages - 1] = 1; + c_ts[num_stages] = 0.5 + c_ts[num_stages - 1] = 1 # Linear increasing timestep for remainder - for i = 2:NumStages-2 - c_ts[i] = cS2 * (i-1)/(NumStages - 3); + for i in 2:num_stages-2 + c_ts[i] = c_s2 * (i-1)/(num_stages - 3) end return c_ts end -function ComputePERK3_ButcherTableau(NumStages, NumStageEvals, semi::AbstractSemidiscretization) +function compute_PERK3_butcher_tableau(num_stages, num_stage_evals, semi::AbstractSemidiscretization) - #Initialize array of c - c_ts = c_Own(1.0, NumStages) + # Initialize array of c + c_ts = c_own(1.0, num_stages) - #Initialize the array of our solution - aUnknown = zeros(NumStageEvals) + # Initialize the array of our solution + a_unknown = zeros(num_stage_evals) - #Case of e = 3 - if NumStageEvals == 3 - aUnknown = [0, c_ts[2], 0.25] + # Special case of e = 3 + if num_stage_evals == 3 + a_unknown = [0, c_ts[2], 0.25] else - #Calculate MonCoeffs from polynomial Optimizer - ConsOrder = 3 - dtMax = 1.0 - dtEps = 1e-9 + # Calculate coefficients of the stability polynomial in monomial form + cons_order = 3 + dt_max = 1.0 + dt_eps = 1e-9 filter_threshold = 1e-12 - #Get EigVals + # Compute spectrum J = jacobian_ad_forward(semi) - EigVals = eigvals(J) - NumEigVals, EigVals = filter_Eigvals(EigVals, filter_threshold) + eig_vals = eigvals(J) + num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_threshold) - MonCoeffs, dt, _ = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) - MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + mon_coeffs, dt, _ = bisection(cons_order, num_eig_vals, num_stages, dt_max, dt_eps, eig_vals) + mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs) - #Define the objective_function + # Define the objective_function function objective_function(x) - return F_own(x, NumStages, NumStageEvals, MonCoeffs) + return f_own(x, num_stages, num_stage_evals, mon_coeffs) end - #call nlsolver to solve until the answer is not NaN or negative values + # Call nlsolver to solve repeatedly until the result is not NaN or negative values is_sol_valid = false - while is_sol_valid == false - #Initialize initial guess - x0 = 0.1 .* rand(NumStageEvals) + while !is_sol_valid + # Initialize initial guess + x0 = 0.1 .* rand(num_stage_evals) x0[1] = 0.0 x0[2] = c_ts[2] - sol = nlsolve(objective_function, method = :trust_region, x0, - ftol = 1e-15, iterations = 10^4, xtol = 1e-13) + sol = nlsolve(objective_function, x0, method = :trust_region, ftol = 1e-15, iterations = 10^4, xtol = 1e-13) - aUnknown = sol.zero + a_unknown = sol.zero - #First check the values of a[i, i-1] which are acquired by nonlinear solver and c[i] - a[i, i-1]. - is_sol_valid = all(x -> !isnan(x) && x >= 0, aUnknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- aUnknown[3:end]) + # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 + is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && all(x -> !isnan(x) && x >= 0 , c_ts[3:end] .- a_unknown[3:end]) end end - println("aUnknown") - println(aUnknown[3:end]) #To debug + println("a_unknown") + println(a_unknown[3:end]) # To debug - AMatrix = zeros(NumStages - 2, 2) - AMatrix[:, 1] = c_ts[3:end] - AMatrix[:, 1] -= aUnknown[3:end] - AMatrix[:, 2] = aUnknown[3:end] + a_matrix = zeros(num_stages - 2, 2) + a_matrix[:, 1] = c_ts[3:end] + a_matrix[:, 1] -= a_unknown[3:end] + a_matrix[:, 2] = a_unknown[3:end] - return AMatrix, c_ts + return a_matrix, c_ts end -function ComputePERK3_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, cS2::Float64) +function compute_PERK3_butcher_tableau(num_stages::Int, base_path_mon_coeffs::AbstractString, c_s2::Float64) - # c Vector form Butcher Tableau (defines timestep per stage) - c = zeros(NumStages) - for k in 2:NumStages-2 - c[k] = cS2 * (k - 1)/(NumStages - 3) # Equidistant timestep distribution (similar to PERK2) - end + # c Vector form Butcher Tableau (defines timestep per stage) + c = zeros(num_stages) + for k in 2:num_stages-2 + c[k] = c_s2 * (k - 1)/(num_stages - 3) # Equidistant timestep distribution (similar to PERK2) + end - # Original third order proposed PERK - #= - c[NumStages - 1] = 1.0/3.0 - c[NumStages] = 1.0 - =# - # Own third order PERK based on SSPRK33 - c[NumStages - 1] = 1.0 - c[NumStages] = 0.5 - - println("Timestep-split: "); display(c); println("\n") + # Original third order proposed PERK + #= + c[num_stages - 1] = 1.0/3.0 + c[num_stages] = 1.0 + =# + # Own third order PERK based on SSPRK33 + c[num_stages - 1] = 1.0 + c[num_stages] = 0.5 + + println("Timestep-split: "); display(c); println("\n") - # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) - CoeffsMax = NumStages - 2 + # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) + coeffs_max = num_stages - 2 - AMatrix = zeros(CoeffsMax, 2) - AMatrix[:, 1] = c[3:end] + a_matrix = zeros(coeffs_max, 2) + a_matrix[:, 1] = c[3:end] - PathMonCoeffs = BasePathMonCoeffs * "a_" * string(NumStages) * "_" * string(NumStages) * ".txt" - NumMonCoeffs, A = read_file(PathMonCoeffs, Float64) - @assert NumMonCoeffs == CoeffsMax + path_mon_coeffs = base_path_mon_coeffs * "a_" * string(num_stages) * "_" * string(num_stages) * ".txt" + num_mon_coeffs, A = read_file(path_mon_coeffs, Float64) + @assert num_mon_coeffs == coeffs_max - AMatrix[:, 1] -= A - AMatrix[:, 2] = A + a_matrix[:, 1] -= A + a_matrix[:, 2] = A - println("A matrix: "); display(AMatrix); println() + println("A matrix: "); display(AMatrix); println() - return AMatrix, c + return a_matrix, c end """ @@ -163,31 +162,31 @@ CarpenterKennedy2N{54, 43} methods. """ mutable struct PERK3 <: PERKSingle - const NumStages::Int + const num_stages::Int - AMatrix::Matrix{Float64} - c::Vector{Float64} + a_matrix::Matrix{Float64} + c::Vector{Float64} - # Constructor for previously computed A Coeffs - function PERK3(NumStages_::Int, BasePathMonCoeffs_::AbstractString, cS2_::Float64=1.0) + # Constructor for previously computed A Coeffs + function PERK3(num_stages_::Int, base_path_mon_coeffs_::AbstractString, c_s2_::Float64=1.0) + newPERK3 = new(num_stages_) - newPERK3 = new(NumStages_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, base_path_mon_coeffs_, c_s2_) + + return newPERK3 + end - newPERK3.AMatrix, newPERK3.c = - ComputePERK3_ButcherTableau(NumStages_, BasePathMonCoeffs_, cS2_) - return newPERK3 - end + # Constructor that computes Butcher matrix A coefficients from a semidiscretization + function PERK3(num_stages_::Int, num_stage_evals_ ::Int, semi_::AbstractSemidiscretization) - #Constructor that compute A coefficients from semidiscretization - function PERK3(NumStages_::Int, NumStageEvals_ ::Int, semi_::AbstractSemidiscretization) + newPERK3 = new(num_stages_) - newPERK3 = new(NumStages_) + newPERK3.a_matrix, newPERK3.c = compute_PERK3_butcher_tableau(num_stages_, num_stage_evals_, semi_) - newPERK3.AMatrix, newPERK3.c = - ComputePERK3_ButcherTableau(NumStages_, NumStageEvals_, semi_) + display(newPERK3.a_matrix) - return newPERK3 - end + return newPERK3 + end end # struct PERK3 diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index 1ad29cd2f0c..acfb9128f71 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -3,139 +3,143 @@ using Convex using ECOS const MOI = Convex.MOI -function filter_Eigvals(EigVals::Array{Complex{Float64}}, threshold:: Float64) +function filter_eigvals(eig_vals::Array{Complex{Float64}}, threshold:: Float64) + filtered_eigvals_counter = 0 - #FilteredEigVals is not an Array but the code seems to work fine regardless - FilteredEigVals = Complex{Float64}[] - for EigVal in EigVals - if abs(EigVal) < threshold + filtered_eig_vals = Complex{Float64}[] + + for eig_val in eig_vals + if abs(eig_val) < threshold filtered_eigvals_counter += 1 else - push!(FilteredEigVals, EigVal) + push!(filtered_eig_vals, eig_val) end end println("Total number of $filtered_eigvals_counter eigenvalues has not been recorded because they were smaller than $threshold") - return length(FilteredEigVals), FilteredEigVals + return length(filtered_eig_vals), filtered_eig_vals end -function ReadInEigVals(Path_toEvalFile::AbstractString) +function read_in_eig_vals(path_to_eval_file::AbstractString) + + # Declare and set to some value + num_eig_vals = -1 + + open(path_to_eval_file, "r") do eval_file + num_eig_vals = countlines(eval_file) + end - NumEigVals = -1 # Declare and set to some value - open(Path_toEvalFile, "r") do EvalFile - NumEigVals = countlines(EvalFile) - end + eig_vals = Array{Complex{Float64}}(undef, num_eig_vals) - EigVals = Array{Complex{Float64}}(undef, NumEigVals) + line_index = 0 - LineIndex = 0 - open(Path_toEvalFile, "r") do EvalFile - # read till end of file - while !eof(EvalFile) + open(path_to_eval_file, "r") do eval_file + # Read till end of file + while !eof(eval_file) - # read a new / next line for every iteration - LineContent = readline(EvalFile) + # Read a new / next line for every iteration + line_content = readline(eval_file) - EigVals[LineIndex + 1] = parse(Complex{Float64}, LineContent) + eig_vals[line_index + 1] = parse(Complex{Float64}, line_content) - LineIndex += 1 + line_index += 1 + end end - end - return NumEigVals, EigVals + return num_eig_vals, eig_vals end -function Polynoms(ConsOrder::Int, NumStageEvals::Int, NumEigVals::Int, normalized_powered_EigvalsScaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) - - for i in 1:NumEigVals - pnoms[i] = 1.0 - end +function polynoms(cons_order::Int, num_stage_evals::Int, num_eig_vals::Int, normalized_powered_eigvals_scaled::Array{Complex{Float64}}, pnoms::Array{Complex{Float64}}, gamma::Variable) + for i in 1:num_eig_vals + pnoms[i] = 1.0 + end - for k in 1:ConsOrder - for i in 1:NumEigVals - pnoms[i] += normalized_powered_EigvalsScaled[i,k] + for k in 1:cons_order + for i in 1:num_eig_vals + pnoms[i] += normalized_powered_eigvals_scaled[i,k] + end end - end - for k in ConsOrder + 1:NumStageEvals - pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] - end + for k in cons_order + 1:num_stage_evals + pnoms += gamma[k - cons_order] * normalized_powered_eigvals_scaled[:,k] + end - return maximum(abs(pnoms)) + return maximum(abs(pnoms)) end -function Bisection(ConsOrder::Int, NumEigVals::Int, NumStageEvals::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) +function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_max::Float64, dt_eps::Float64, eig_vals::Array{Complex{Float64}}) - dtMin = 0.0 + dt_min = 0.0 - dt = -1.0 + dt = -1.0 - AbsP = -1.0 + abs_p = -1.0 - pnoms = ones(Complex{Float64}, NumEigVals, 1) + pnoms = ones(Complex{Float64}, num_eig_vals, 1) - gamma = Variable(NumStageEvals - ConsOrder) # Init datastructure for results + # Init datastructure for results + gamma = Variable(num_stage_evals - cons_order) - normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStageEvals) + normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - for j in 1:NumStageEvals - fac_j = factorial(j) - for i in 1:NumEigVals - normalized_powered_Eigvals[i, j] = EigVals[i]^j / fac_j + for j in 1:num_stage_evals + fac_j = factorial(j) + for i in 1:num_eig_vals + normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j + end end - end - normalized_powered_EigvalsScaled = zeros(Complex{Float64}, NumEigVals, NumStageEvals) + normalized_powered_eigvals_scaled = zeros(Complex{Float64}, num_eig_vals, num_stage_evals) - while dtMax - dtMin > dtEps - dt = 0.5 * (dtMax + dtMin) + while dt_max - dt_min > dt_eps + dt = 0.5 * (dt_max + dt_min) - for k in 1:NumStageEvals - dt_k = dt^k - for i in 1:NumEigVals - normalized_powered_EigvalsScaled[i,k] = dt_k * normalized_powered_Eigvals[i,k] - end - end + for k in 1:num_stage_evals + dt_k = dt^k + for i in 1:num_eig_vals + normalized_powered_eigvals_scaled[i,k] = dt_k * normalized_powered_eigvals[i,k] + end + end - # Use last optimal values for gamm0 in (potentially) next iteration - problem = minimize(Polynoms(ConsOrder, NumStageEvals, NumEigVals, normalized_powered_EigvalsScaled, pnoms, gamma)) - - Convex.solve!( - problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, - "delta" => 2e-7, - #"eps" => 1e9, # 1e-13 - "feastol" => 1e-9, # 1e-9 - "abstol" => 1e-9, # 1e-9 - "reltol" => 1e-9, # 1e-9 - "feastol_inacc" => 1e-4, - "abstol_inacc" => 5e-5, - "reltol_inacc" => 5e-5, - "nitref" => 9, - "maxit" => 100, - "verbose" => 3); silent_solver = false - ) - - AbsP = problem.optval - - println("Current MaxAbsP: ", AbsP, "\nCurrent dt: ", dt, "\n") - - if AbsP < 1.0 - dtMin = dt - else - dtMax = dt + # Use last optimal values for gamm0 in (potentially) next iteration + problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma)) + + Convex.solve!( + problem, + # Parameters taken from default values for EiCOS + MOI.OptimizerWithAttributes(ECOS.Optimizer, "gamma" => 0.99, + "delta" => 2e-7, + #"eps" => 1e9, # 1e-13 + "feastol" => 1e-9, # 1e-9 + "abstol" => 1e-9, # 1e-9 + "reltol" => 1e-9, # 1e-9 + "feastol_inacc" => 1e-4, + "abstol_inacc" => 5e-5, + "reltol_inacc" => 5e-5, + "nitref" => 9, + "maxit" => 100, + "verbose" => 3); silent_solver = false + ) + + abs_p = problem.optval + + println("Current MaxAbsP: ", abs_p, "\nCurrent dt: ", dt, "\n") + + if abs_p < 1.0 + dt_min = dt + else + dt_max = dt + end end - end - return evaluate(gamma), AbsP, dt + return evaluate(gamma), abs_p, dt end -function undo_normalization(ConsOrder::Int, NumStageEvals::Int, gammaOpt) - for k in ConsOrder + 1:NumStageEvals - gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) - end - return gammaOpt +function undo_normalization!(cons_order::Int, num_stage_evals::Int, gamma_opt) + for k in cons_order + 1:num_stage_evals + gamma_opt[k - cons_order] = gamma_opt[k - cons_order] / factorial(k) + end + return gamma_opt end \ No newline at end of file From 483015bdafdca28ff5b9e445a496f5eb5239e488 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 15 Feb 2024 20:52:25 +0100 Subject: [PATCH 4/4] spelling fix --- .../paired_explicit_runge_kutta/polynomial_optimizer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl index acfb9128f71..534ee883909 100644 --- a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -103,7 +103,7 @@ function bisection(cons_order::Int, num_eig_vals::Int, num_stage_evals::Int, dt_ end end - # Use last optimal values for gamm0 in (potentially) next iteration + # Use last optimal values for gamma in (potentially) next iteration problem = minimize(polynoms(cons_order, num_stage_evals, num_eig_vals, normalized_powered_eigvals_scaled, pnoms, gamma)) Convex.solve!(