constructor changed

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,

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_

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