From 3609f39b979de2e894bb78ecb9c67e8fe75ddf9f Mon Sep 17 00:00:00 2001
From: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com>
Date: Mon, 25 Nov 2024 10:50:13 +0100
Subject: [PATCH] Modification on PERK time integrator to handle cases where s
 = e (#2177)

* handle s=consistency_order

* minor change
---
 ext/TrixiConvexECOSExt.jl                     | 12 +++++++--
 .../methods_PERK2.jl                          | 17 +++++++------
 .../methods_PERK3.jl                          | 25 ++++++++++---------
 3 files changed, 32 insertions(+), 22 deletions(-)

diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl
index 8251fe3eed9..9b897436dbb 100644
--- a/ext/TrixiConvexECOSExt.jl
+++ b/ext/TrixiConvexECOSExt.jl
@@ -60,7 +60,11 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals,
     end
 
     # For optimization only the maximum is relevant
-    return maximum(abs(pnoms))
+    if consistency_order - num_stage_evals == 0
+        return maximum(abs.(pnoms)) # If there is no variable to optimize, we need to use the broadcast operator.
+    else
+        return maximum(abs(pnoms))
+    end
 end
 
 #=
@@ -151,7 +155,11 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
         println("Concluded stability polynomial optimization \n")
     end
 
-    gamma_opt = evaluate(gamma)
+    if consistency_order - num_stage_evals != 0
+        gamma_opt = evaluate(gamma)
+    else
+        gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing.
+    end
 
     # Catch case S = 3 (only one opt. variable)
     if isa(gamma_opt, Number)
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 0ff648556c4..2451680a505 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
@@ -49,15 +49,16 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
                                                           dtmax,
                                                           dteps,
                                                           eig_vals; verbose)
-    monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
-                                          num_stages)
 
-    num_monomial_coeffs = length(monomial_coeffs)
-    @assert num_monomial_coeffs == coeffs_max
-    A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs)
-
-    a_matrix[:, 1] -= A
-    a_matrix[:, 2] = A
+    if num_stages != consistency_order
+        monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
+                                              num_stages)
+        num_monomial_coeffs = length(monomial_coeffs)
+        @assert num_monomial_coeffs == coeffs_max
+        A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs)
+        a_matrix[:, 1] -= A
+        a_matrix[:, 2] = A
+    end
 
     return a_matrix, c, dt_opt
 end
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 477748b28c5..02bc3eeba36 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
@@ -33,21 +33,22 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan,
     # Initialize the array of our solution
     a_unknown = zeros(num_stages - 2)
 
+    # Calculate coefficients of the stability polynomial in monomial form
+    consistency_order = 3
+    dtmax = tspan[2] - tspan[1]
+    dteps = 1.0f-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, dteps,
+                                                          eig_vals; verbose)
+
     # Special case of e = 3
-    if num_stages == 3
+    if num_stages == consistency_order
         a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau
     else
-        # Calculate coefficients of the stability polynomial in monomial form
-        consistency_order = 3
-        dtmax = tspan[2] - tspan[1]
-        dteps = 1.0f-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, dteps,
-                                                              eig_vals; verbose)
         monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
                                               num_stages)