diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl index 07aa047d751..7c22fed20c0 100644 --- a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl @@ -47,12 +47,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # run the simulation -MonCoeffsFile = "gamma_6.txt" - -download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", -MonCoeffsFile) - -ode_algorithm = PERK2(6, "./") +ode_algorithm = PERK2(6, semi) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = Trixi.solve(ode, ode_algorithm, 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 3d247fd65e4..b13b8204f6b 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -2,6 +2,7 @@ # 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 @@ -22,7 +23,7 @@ function ComputeACoeffs(NumStageEvals::Int, return reverse(ACoeffs) end -function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::AbstractString, bS::Float64, cEnd::Float64) +function ComputePERK2_ButcherTableau(NumStages::Int, semi::AbstractSemidiscretization, bS::Float64, cEnd::Float64) # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(NumStages) @@ -39,8 +40,20 @@ function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::Abstract AMatrix[:, 1] = c[3:end] - PathMonCoeffs = BasePathMonCoeffs * "gamma_" * string(NumStages) * ".txt" - NumMonCoeffs, MonCoeffs = read_file(PathMonCoeffs, Float64) + ConsOrder = 2 + filter_thres = 1e-12 + dtMax = 1.0 + dtEps = 1e-9 + + J = jacobian_ad_forward(semi) + EigVals = eigvals(J) + + NumEigVals, EigVals = filter_Eigvals(EigVals, filter_thres) + + MonCoeffs, PWorstCase, dtOpt = Bisection(ConsOrder, NumEigVals, NumStages, dtMax, dtEps, EigVals) + MonCoeffs = undo_normalization(ConsOrder, NumStages, MonCoeffs) + + NumMonCoeffs = length(MonCoeffs) @assert NumMonCoeffs == CoeffsMax A = ComputeACoeffs(NumStages, SE_Factors, MonCoeffs) @@ -81,12 +94,12 @@ mutable struct PERK2 <: PERKSingle cEnd::Float64 # Constructor for previously computed A Coeffs - function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5) + function PERK2(NumStages_::Int, semi_::AbstractSemidiscretization, bS_::Float64=1.0, cEnd_::Float64=0.5) newPERK2 = new(NumStages_) newPERK2.AMatrix, newPERK2.c = - ComputePERK2_ButcherTableau(NumStages_, BasePathMonCoeffs_, bS_, cEnd_) + ComputePERK2_ButcherTableau(NumStages_, semi_, bS_, cEnd_) newPERK2.b1 = one(bS_) - bS_ newPERK2.bS = bS_ diff --git a/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl new file mode 100644 index 00000000000..9bf795d2b3e --- /dev/null +++ b/src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl @@ -0,0 +1,147 @@ +using LinearAlgebra +using Convex +using ECOS +const MOI = Convex.MOI + +function filter_Eigvals(EigVals::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_eigvals_counter += 1 + else + push!(FilteredEigVals, EigVal) + end + end + + println("Total number of $filtered_eigvals_counter eigenvalues has not been recorded because they were smaller than $threshold") + + return length(FilteredEigVals), FilteredEigVals +end + +function ReadInEigVals(Path_toEvalFile::AbstractString) + + NumEigVals = -1 # Declare and set to some value + open(Path_toEvalFile, "r") do EvalFile + NumEigVals = countlines(EvalFile) + end + + EigVals = Array{Complex{Float64}}(undef, NumEigVals) + + LineIndex = 0 + open(Path_toEvalFile, "r") do EvalFile + # read till end of file + while !eof(EvalFile) + + # read a new / next line for every iteration + LineContent = readline(EvalFile) + + EigVals[LineIndex + 1] = parse(Complex{Float64}, LineContent) + + LineIndex += 1 + end + end + + return NumEigVals, EigVals +end + +function Polynoms(ConsOrder::Int, NumStages::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 + + for k in 1:ConsOrder + for i in 1:NumEigVals + pnoms[i] += normalized_powered_EigvalsScaled[i,k] + end + #pnoms += 1/factorial(k) * EigValsScaled.^k + end + + for k in ConsOrder + 1:NumStages + #Attemp to use nested loop instead of [:,k] has been made. Convex.jl doesn't accept it. + pnoms += gamma[k - ConsOrder] * normalized_powered_EigvalsScaled[:,k] + end + + return maximum(abs(pnoms)) +end + + +function Bisection(ConsOrder::Int, NumEigVals::Int, NumStages::Int, dtMax::Float64, dtEps::Float64, EigVals::Array{Complex{Float64}}) + #dtMin = Base.big(0) + dtMin = 0.0 + + #dt = Base.big(-1.0) + dt = -1.0 + + #AbsP = Base.big(-1.0) + AbsP = -1.0 + + pnoms = ones(Complex{Float64}, NumEigVals, 1) + + gamma = Variable(NumStages - ConsOrder) # Init datastructure for results + + normalized_powered_Eigvals = zeros(Complex{Float64}, NumEigVals, NumStages) + + for j in 1:NumStages + 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) + + while dtMax - dtMin > dtEps + #dt = Base.big(0.5) * (dtMax + dtMin) + dt = 0.5 * (dtMax + dtMin) + + for k in 1:NumStages + dt_k = dt^k + for i in 1:NumEigVals + normalized_powered_EigvalsScaled[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, NumStages, 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 < Base.big(1.0) + if AbsP < 1.0 + dtMin = dt + else + dtMax = dt + end + end + + return evaluate(gamma), AbsP, dt +end + +function undo_normalization(ConsOrder::Int, NumStages::Int, gammaOpt) + for k in ConsOrder + 1:NumStages + gammaOpt[k - ConsOrder] = gammaOpt[k - ConsOrder] / factorial(k) + end + return gammaOpt +end \ No newline at end of file