Modification on PERK time integrator to handle cases where s = e

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)

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

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)