something wrong with b values... will have to reconstruct that part a…

…gain from the stability polynomial in TrixiNLsolveExt

ext/TrixiNLsolveExt.jl

@@ -86,19 +86,23 @@ function EmbeddedPairedExplicitRK3_butcher_tableau_objective_function!(b_eq, x,
         0
     ]
 
+    println("Length of b_coeff")
+    println(length(b_coeff)) # Debugging
+
     b_eq_count = 0
 
     #TODO: check the logic of this loop cos b is getting absurdly large
     for i in 3:num_stage_evals-1
         sum = 0.0
+        fac_i = factorial(i)
         for j in (i + num_stages - num_stage_evals):num_stages-1
             prod = 1.0
             for k in (3 + j - i):j
                 prod *= a[k]
             end
-            sum += prod * b_coeff[j] * c[j - i + 2]
+            sum += prod * b_coeff[j] * c[j - i + 2] 
         end
-        b_eq[i-2] = embedded_monomial_coeffs[i-2] - sum
+        b_eq[i-2] = embedded_monomial_coeffs[i-2] - sum * fac_i
         b_eq_count += 1
     end
 
@@ -174,6 +178,8 @@ end
 function Trixi.solve_b_butcher_coeffs_unknown!(b_unknown, num_stages, num_stage_evals, embedded_monomial_coeffs, c, a_unknown;
     verbose, max_iter = 100000)
 
+    verbose = true
+
     # Construct a full a coefficient vector
     a = zeros(num_stages)
     num_a_unknown = length(a_unknown)
@@ -201,17 +207,22 @@ function Trixi.solve_b_butcher_coeffs_unknown!(b_unknown, num_stages, num_stage_
     # Due to the nature of the nonlinear solver, different initial guesses can lead to 
     # small numerical differences in the solution.
 
-    # There is e-2 free variables of b of the embedded scheme
-    x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stage_evals - 2)
+    for _ in 1:max_iter
+
+        # There is e-2 free variables of b of the embedded scheme
+        x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stage_evals - 2)
+
+        sol = nlsolve(embedded_scheme_objective_function!, x0, method = :trust_region,
+                    ftol = 4.0e-16, # Enforce objective up to machine precision
+                    iterations = 10^4, xtol = 1.0e-13, autodiff = :forward)
 
-    sol = nlsolve(embedded_scheme_objective_function!, x0, method = :newton,
-                  ftol = 4.0e-16, # Enforce objective up to machine precision
-                  iterations = 10^4, xtol = 1.0e-13, autodiff = :forward)
+        b_unknown = sol.zero # Retrieve solution (root = zero)
 
-    b_unknown = sol.zero # Retrieve solution (root = zero)
+        is_sol_valid = all(x -> !isnan(x) && x >= 0, b_unknown) && (sum(b_unknown) <= 1.0)
 
-
-    return b_unknown
+        return b_unknown
+
+    end
 end
 end # @muladd
 

src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl

@@ -17,7 +17,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
 
     # Initialize the array of our solution
     a_unknown = zeros(num_stage_evals - 2)
-    b_unknown = zeros(num_stage_evals - 2)
+    b = zeros(num_stage_evals - 2)
 
     # Special case of e = 3
     if num_stage_evals == 3
@@ -54,14 +54,13 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
 
         println("dt_opt_b = ", dt_opt_b)
 
-        b_unknown = Trixi.solve_b_butcher_coeffs_unknown!(b_unknown, num_stages, num_stage_evals, embedded_monomial_coeffs, c, a_unknown; verbose)
+        b = Trixi.solve_b_butcher_coeffs_unknown!(b, num_stages, num_stage_evals, embedded_monomial_coeffs, c, a_unknown; verbose)
 
         println("Embedded monomial after normalization: ", embedded_monomial_coeffs)
-        println("This is b_unknown = ",b_unknown)
         println("dt_opt_a = ", dt_opt_a)
 
 
-        error("Stability polynomial optimization process is complete.")
+        #error("Stability polynomial optimization process is complete.")
     end
     # Fill A-matrix in P-ERK style
     a_matrix = zeros(num_stage_evals - 2, 2)