add tspan as a parameter and adjust the elixir accordingly

DanielDoehring · Mar 8, 2024 · 88e8168 · 88e8168
1 parent e7060e5
commit 88e8168
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 39 deletions.
diff --git a/examples/tree_1d_dgsem/elixir_advection_PERK2.jl b/examples/tree_1d_dgsem/elixir_advection_PERK2.jl
@@ -46,9 +46,12 @@ callbacks = CallbackSet(summary_callback, analysis_callback,
 ###############################################################################
 # run the simulation
 
+# Define tspan to calculate maximum time step allowed for the bisection algorithm used in calculate polynomial coefficients in ODE algorithm
+tspan = [0.0, 1.0]
+
 # Construct second order P-ERK method with 6 stages for 
 # given simulation setup
-ode_algorithm = PERK2(6, semi)
+ode_algorithm = PERK2(6, tspan, semi)
 
 sol = Trixi.solve(ode, ode_algorithm,
                   dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
@@ -58,4 +61,4 @@ sol = Trixi.solve(ode, ode_algorithm,
 summary_callback()
 
 # Show the error of the elixir
-analysis_callback(sol)
+analysis_callback(sol)
diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
@@ -23,7 +23,7 @@ function compute_a_coeffs(num_stage_evals, se_factors, mon_coeffs)
     return reverse(a_coeffs)
 end
 
-function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization,
+function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretization, tspan,
                                        bS, c_end)
 
     # c Vector form Butcher Tableau (defines timestep per stage)
@@ -41,14 +41,14 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat
 
     cons_order = 2
     filter_thres = 1e-12
-    dtmax = 1.0
-    dt_eps = 1e-9
+    dtmax = tspan[2] - tspan[1]
+    dteps = 1e-9
 
     eig_vals = eigvals(jacobian_ad_forward(semi))
 
     num_eig_vals, eig_vals = filter_eigvals(eig_vals, filter_thres)
 
-    mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dt_eps,
+    mon_coeffs, dt_opt = bisection(cons_order, num_eig_vals, num_stages, dtmax, dteps,
                                    eig_vals)
     mon_coeffs = undo_normalization!(cons_order, num_stages, mon_coeffs)
 
@@ -59,7 +59,6 @@ function compute_PERK2_Butcher_tableau(num_stages, semi::AbstractSemidiscretizat
     a_matrix[:, 1] -= A
     a_matrix[:, 2] = A
 
-
     return a_matrix, c
 end
 
@@ -87,7 +86,6 @@ function compute_PERK2_Butcher_tableau(num_stages, base_path_mon_coeffs::Abstrac
     a_matrix[:, 1] -= A
     a_matrix[:, 2] = A
 
-
     return a_matrix, c
 end
 
@@ -128,11 +126,12 @@ mutable struct PERK2 <: PERKSingle
     end
 
     #Constructor that calculate the coefficients with polynomial optimizer
-    function PERK2(num_stages, semi::AbstractSemidiscretization, bS = 1.0, c_end = 0.5)
+    function PERK2(num_stages, tspan, semi::AbstractSemidiscretization, bS = 1.0,
+                   c_end = 0.5)
         newPERK2 = new(num_stages)
 
         newPERK2.a_matrix, newPERK2.c = compute_PERK2_Butcher_tableau(num_stages,
-                                                                      semi,
+                                                                      semi, tspan,
                                                                       bS, c_end)
 
         newPERK2.b1 = one(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
@@ -3,27 +3,6 @@ using ECOS
 using Convex
 const MOI = Convex.MOI
 
-function read_in_eig_vals(path_to_eval_file)
-    # 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
-    eig_vals = Array{Complex{Float64}}(undef, num_eig_vals)
-    line_index = 0
-    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          
-            line_content = readline(eval_file)
-            eig_vals[line_index + 1] = parse(Complex{Float64}, line_content)
-            line_index += 1
-        end
-    end
-
-    return num_eig_vals, eig_vals
-end
-
 function filter_eigvals(eig_vals, threshold)
     filtered_eigvals_counter = 0
     filtered_eig_vals = Complex{Float64}[]
@@ -60,8 +39,8 @@ function polynoms(cons_order, num_stage_evals, num_eig_vals,
     return maximum(abs(pnoms))
 end
 
-function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, eig_vals)
-    dt_min = 0.0
+function bisection(cons_order, num_eig_vals, num_stage_evals, dtmax, dteps, eig_vals)
+    dtmin = 0.0
     dt = -1.0
     abs_p = -1.0
 
@@ -84,8 +63,8 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei
 
     println("Start optimization of stability polynomial \n")
 
-    while dt_max - dt_min > dt_eps
-        dt = 0.5 * (dt_max + dt_min)
+    while dtmax - dtmin > dteps
+        dt = 0.5 * (dtmax + dtmin)
 
         for k in 1:num_stage_evals
             dt_k = dt^k
@@ -116,12 +95,10 @@ function bisection(cons_order, num_eig_vals, num_stage_evals, dt_max, dt_eps, ei
 
         abs_p = problem.optval
 
-        println("MaxAbsP: ", abs_p, "\ndt: ", dt, "\n")
-
         if abs_p < 1.0
-            dt_min = dt
+            dtmin = dt
         else
-            dt_max = dt
+            dtmax = dt
         end
     end