From f4274ee7d098331d5639758ecc891afc0be479cd Mon Sep 17 00:00:00 2001
From: Warisa <compund555@gmail.com>
Date: Sat, 6 Jul 2024 09:59:17 +0200
Subject: [PATCH] add optimum time step for special case and add some stuff in
 bisection to make this work as well

---
 ext/TrixiConvexECOSExt.jl                     | 63 ++++++++++++-------
 .../methods_PERK3.jl                          | 19 +++---
 2 files changed, 53 insertions(+), 29 deletions(-)

diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl
index 948dbf103cd..43a88413bb3 100644
--- a/ext/TrixiConvexECOSExt.jl
+++ b/ext/TrixiConvexECOSExt.jl
@@ -54,6 +54,10 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals,
         end
     end
 
+    if consistency_order == num_stage_evals
+        return maximum(abs.(pnoms))
+    end
+
     # Contribution from free coefficients
     for k in (consistency_order + 1):num_stage_evals
         pnoms += gamma[k - consistency_order] * normalized_powered_eigvals_scaled[:, k]
@@ -118,27 +122,36 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
             end
         end
 
-        # Use last optimal values for gamma in (potentially) next iteration
-        problem = minimize(stability_polynomials!(pnoms, consistency_order,
-                                                  num_stage_evals,
-                                                  normalized_powered_eigvals_scaled,
-                                                  gamma))
-
-        solve!(problem,
-               # Parameters taken from default values for EiCOS
-               MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99,
-                                           "delta" => 2e-7,
-                                           "feastol" => 1e-9,
-                                           "abstol" => 1e-9,
-                                           "reltol" => 1e-9,
-                                           "feastol_inacc" => 1e-4,
-                                           "abstol_inacc" => 5e-5,
-                                           "reltol_inacc" => 5e-5,
-                                           "nitref" => 9,
-                                           "maxit" => 100,
-                                           "verbose" => 3); silent = true)
-
-        abs_p = problem.optval
+        stability_polynomials = stability_polynomials!(pnoms, consistency_order,
+                                                       num_stage_evals,
+                                                       normalized_powered_eigvals_scaled,
+                                                       gamma)
+
+        if num_stage_evals != consistency_order
+
+            # Use last optimal values for gamma in (potentially) next iteration
+            problem = minimize(stability_polynomials)
+
+            solve!(problem,
+                # Parameters taken from default values for EiCOS
+                MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99,
+                                            "delta" => 2e-7,
+                                            "feastol" => 1e-9,
+                                            "abstol" => 1e-9,
+                                            "reltol" => 1e-9,
+                                            "feastol_inacc" => 1e-4,
+                                            "abstol_inacc" => 5e-5,
+                                            "reltol_inacc" => 5e-5,
+                                            "nitref" => 9,
+                                            "maxit" => 100,
+                                            "verbose" => 3); silent = true)
+
+            abs_p = problem.optval
+            gamma_opt = evaluate(gamma)
+        else
+            abs_p = stability_polynomials
+            gamma_opt = []
+        end
 
         if abs_p < 1
             dtmin = dt
@@ -151,7 +164,13 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
         println("Concluded stability polynomial optimization \n")
     end
 
-    return evaluate(gamma), dt
+    
+
+    if isa(gamma_opt, Number)
+        gamma_opt = [gamma_opt]
+    end
+
+    return gamma_opt, dt
 end
 end # @muladd
 
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 e66570c79f1..5fe5e5415fc 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
@@ -80,18 +80,23 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan,
     # Initialize the array of our solution
     a_unknown = zeros(num_stages - 2)
 
+    consistency_order = 3
+    dtmax = tspan[2] - tspan[1]
+    dteps = 1e-9
+
+    num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose)
+
     # Special case of e = 3
     if num_stages == 3
+        # Use bisection to get the optimal time step
+        _, dt_opt = bisect_stability_polynomial(consistency_order,
+                                                num_eig_vals, num_stages,
+                                                dtmax,
+                                                dteps,
+                                                eig_vals; verbose)
         a_unknown = [0.25]
-        dt_opt = 42.0 # TODO! This is a placeholder value
     else
         # Calculate coefficients of the stability polynomial in monomial form
-        consistency_order = 3
-        dtmax = tspan[2] - tspan[1]
-        dteps = 1e-9
-
-        num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose)
-
         monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order,
                                                               num_eig_vals, num_stages,
                                                               dtmax,