Add constructor of PERK3

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
@@ -166,7 +165,7 @@
   callback::Callback # callbacks; used in Trixi
   adaptive::Bool # whether the algorithm is adaptive; ignored
   dtmax::Float64 # ignored
   maxiters::Int # maximal numer of time steps
   tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored
 end
 

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
 
 

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

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,15 +58,15 @@ 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
 
   return maximum(abs(pnoms))
 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,31 +76,31 @@ 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]
       end
     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