some changes. Set constraint to b positive

examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl

@@ -56,7 +56,7 @@ save_solution = SaveSolutionCallback(dt = 0.1,
 # Construct embedded order paired explicit Runge-Kutta method with 10 stages and 7 evaluation stages for given simulation setup.
 # Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used 
 # in calculating the polynomial coefficients in the ODE algorithm.
-ode_algorithm = Trixi.EmbeddedPairedRK3(16, 5, tspan, semi)
+ode_algorithm = Trixi.EmbeddedPairedRK3(16, 11, tspan, semi)
 
 # Calculate the CFL number for the given ODE algorithm and ODE problem (cfl_number calculate from dt_opt of the optimization of
 # b values in the Butcher tableau of the ODE algorithm).

ext/TrixiConvexECOSExt.jl

@@ -15,7 +15,7 @@ end
 using LinearAlgebra: eigvals
 
 # Use functions that are to be extended and additional symbols that are not exported
-using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd
+using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, solve_b_embedded, @muladd
 
 # By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
 # Since these FMAs can increase the performance of many numerical algorithms,
@@ -160,6 +160,140 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
 
     return gamma_opt, dt
 end
+
+function compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+
+    A = zeros(num_stage_evals - 1, num_stage_evals - 1)
+    b_embedded = zeros(num_stage_evals - 1)
+    rhs = [1, 1/2, embedded_monomial_coeffs...]  # Cast embedded_monomial_coeffs to Float64 if needed
+
+    # sum(b) = 1
+    A[1, :] .= 1
+
+    # The second order constraint: dot(b,c) = 1/2
+    for i in 2:num_stage_evals - 1
+        A[2, i] = c[num_stages - num_stage_evals + i]
+    end
+
+    # Fill the A matrix
+    for i in 3:(num_stage_evals - 1)
+        # z^i
+        for j in i: (num_stage_evals - 1)
+            println("i = ", i, ", j = ", j)
+            println("[num_stages - num_stage_evals + j - 1] = ", num_stages - num_stage_evals + j - 1)
+            A[i,j] = c[num_stages - num_stage_evals + j - 1]
+            # number of times a_unknown should be multiplied in each power of z
+            for k in 1: i-2
+                # so we want to get from a[k] * ... i-2 times (1 time is already accounted for by c-coeff)
+                # j-1 - k + 1 = j - k
+                println("a_unknown at index: ", j - k)
+                A[i, j] *= a_unknown[j - k] # a[k] is in the same stage as b[k-1] -> since we also store b_1
+            end
+        end
+        rhs[i] /= factorial(i)
+    end
+
+    A_inv = inv(A)
+
+    b_embedded = A_inv * rhs
+
+    return b_embedded
+end
+
+#TODO: Add an optimization of the embedded scheme that subject the b solved from the set of gamme to be +
+#      the same as the b from the original scheme. This will be done by using the same optimization as the normal scheme but with cosntraint and the function
+#      that will be used to construct the b vector from the gamma vector.
+function Trixi.solve_b_embedded(consistency_order, num_eig_vals, num_stage_evals, num_stages,
+                                           dtmax, dteps, eig_vals, a_unknown, c;
+                                           verbose = false)
+    dtmin = 0.0
+    dt = -1.0
+    abs_p = -1.0
+    consistency_order_embedded = consistency_order - 1
+
+    num_stage_evals_embedded = num_stage_evals -1
+
+    # Construct stability polynomial for each eigenvalue
+    pnoms = ones(Complex{Float64}, num_eig_vals, 1)
+
+    # Init datastructure for monomial coefficients
+    gamma = Variable(num_stage_evals_embedded - consistency_order_embedded)
+
+    normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals_embedded)
+
+    for j in 1:num_stage_evals_embedded
+        fac_j = factorial(j)
+        for i in 1:num_eig_vals
+            normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j
+        end
+    end
+
+    normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals)
+
+    if verbose
+        println("Start optimization of stability polynomial \n")
+    end
+
+    # Bisection on timestep
+    while dtmax - dtmin > dteps
+        dt = 0.5 * (dtmax + dtmin)
+
+        # Compute stability polynomial for current timestep
+        for k in 1:num_stage_evals_embedded
+            dt_k = dt^k
+            for i in 1:num_eig_vals
+                normalized_powered_eigvals_scaled[i, k] = dt_k *
+                                                          normalized_powered_eigvals[i,
+                                                                                     k]
+            end
+        end
+
+        constraint = compute_b_embedded_coeffs(num_stage_evals, num_stages, gamma, a_unknown, c).>= -1e-5
+
+        # Use last optimal values for gamma in (potentially) next iteration
+        problem = minimize(stability_polynomials!(pnoms, consistency_order_embedded,
+                                                  num_stage_evals_embedded,
+                                                  normalized_powered_eigvals_scaled,
+                                                  gamma), constraint)
+
+        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_solver = true)
+
+        abs_p = problem.optval
+
+        if abs_p < 1
+            dtmin = dt
+        else
+            dtmax = dt
+        end
+    end
+
+    if verbose
+        println("Concluded stability polynomial optimization \n")
+    end
+
+    gamma_opt = evaluate(gamma)
+
+    # Catch case S = 3 (only one opt. variable)
+    if isa(gamma_opt, Number)
+        gamma_opt = [gamma_opt]
+    end
+
+    b_embedded = compute_b_embedded_coeffs(num_stage_evals, num_stages, gamma_opt, a_unknown, c)
+
+    return b_embedded, dt
+end
 end # @muladd
 
 end # module TrixiConvexECOSExt

src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl

@@ -7,44 +7,8 @@ using DelimitedFiles: readdlm
 @muladd begin
 #! format: noindent
 
-function compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
-    # Different scheme for c?
-
-    A = zeros(num_stage_evals - 1, num_stage_evals - 1)
-    b_embedded = zeros(num_stage_evals - 1)
-    rhs = [1, 1/2, embedded_monomial_coeffs...]
-
-    # sum(b) = 1
-    A[1, :] .= 1
-
-    # The second order constraint: dot(b,c) = 1/2
-    for i in 2:num_stage_evals - 1
-        A[2, i] = c[num_stages - num_stage_evals + i]
-    end
-
-    # Fill the A matrix
-    for i in 3:(num_stage_evals - 1)
-        # z^i
-        for j in i: (num_stage_evals - 1)
-            println("i = ", i, ", j = ", j)
-            println("[num_stages - num_stage_evals + j - 1] = ", num_stages - num_stage_evals + j - 1)
-            A[i,j] = c[num_stages - num_stage_evals + j - 1]
-            # number of times a_unknown should be multiplied in each power of z
-            for k in 1: i-2
-                # so we want to get from a[k] * ... i-2 times (1 time is already accounted for by c-coeff)
-                # j-1 - k + 1 = j - k
-                println("a_unknown at index: ", j - k)
-                A[i, j] *= a_unknown[j - k] # a[k] is in the same stage as b[k-1] -> since we also store b_1
-            end
-        end
-        #rhs[i] /= factorial(i)
-    end
-
-    display(A)
-
-    b_embedded = A \ rhs
-    return b_embedded
-end
+# To be overwritten in TrixiConvexECOSExt.jl
+function solve_b_embedded end
 
 # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3
 # using a list of eigenvalues
@@ -76,13 +40,6 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
         monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
                                               num_stage_evals)
 
-        embedded_monomial_coeffs, dt_opt_b = bisect_stability_polynomial(consistency_order-1,
-                                                                        num_eig_vals, num_stage_evals-1,
-                                                                        dtmax, dteps,
-                                                                        eig_vals; verbose)
-
-        embedded_monomial_coeffs = undo_normalization!(embedded_monomial_coeffs, consistency_order-1, num_stage_evals-1)
-
         # Solve the nonlinear system of equations from monomial coefficient and
         # Butcher array abscissae c to find Butcher matrix A
         # This function is extended in TrixiNLsolveExt.jl
@@ -91,12 +48,13 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
                                                     monomial_coeffs, cS2, c;
                                                     verbose)
 
-        println("dt_opt_b = ", dt_opt_b)
+        b, dt_opt_b = solve_b_embedded(consistency_order, num_eig_vals, num_stage_evals, num_stages,
+        dtmax, dteps, eig_vals, a_unknown, c;
+        verbose = true)
 
-        b = compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+        println("dt_opt_b = ", dt_opt_b)
 
         println("b: ", b)
-        println("Embedded monomial after normalization: ", embedded_monomial_coeffs)
         println("dt_opt_a = ", dt_opt_a)