Merge branch 'embedded_PERK3' of https://github.com/warisa-r/Trixi.jl …

…into embedded_PERK3

examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl

@@ -56,13 +56,12 @@ 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, 11, tspan, semi)
+ode_algorithm = Trixi.EmbeddedPairedRK3(16, 5, 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).
 #cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
-stepsize_callback = StepsizeCallback() # Warisa: This number is quite small in contrast the other one from optimizing A
-# I've tried using cfl of 1.5 and the error is very similar.
+stepsize_callback = StepsizeCallback()
 
 callbacks = CallbackSet(summary_callback,
                         alive_callback,
@@ -78,24 +77,3 @@ sol = Trixi.solve(ode, ode_algorithm,
 
 # Print the timer summary
 summary_callback()
-
-# Some function defined so that I can check if the second order condition is met. This will be removed later.
-function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded)
-    # Construct the b vector
-    b = [
-        b_unknown[1],
-        zeros(Float64, num_stages_embedded - num_stage_evals_embedded)...,
-        b_unknown[2:end]...,
-        0
-    ]
-    return b
-end
-
-b = construct_b_vector(ode_algorithm.b, ode_algorithm.num_stages - 1,
-                       ode_algorithm.num_stage_evals - 1)
-println("dot(b, c) = ", dot(b, ode_algorithm.c))
-println("sum(b) = ", sum(b))
-
-println("dt_opt_a = ", ode_algorithm.dt_opt_a)
-println("dt_opt_b = ", ode_algorithm.dt_opt_b)
-println("ratio = ", ode_algorithm.dt_opt_b / ode_algorithm.dt_opt_a * 100)

ext/TrixiConvexECOSExt.jl

@@ -3,11 +3,11 @@ module TrixiConvexECOSExt
 
 # Required for coefficient optimization in P-ERK scheme integrators
 if isdefined(Base, :get_extension)
-    using Convex: MOI, solve!, Variable, minimize, evaluate
+    using Convex: MOI, solve!, Variable, minimize, evaluate, sumsquares
     using ECOS: Optimizer
 else
     # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
-    using ..Convex: MOI, solve!, Variable, minimize, evaluate
+    using ..Convex: MOI, solve!, Variable, minimize, evaluate, sumsquares
     using ..ECOS: Optimizer
 end
 
@@ -16,6 +16,7 @@ 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, solve_b_embedded, @muladd
+using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, compute_b_embedded_coeffs, @muladd
 
 # By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
 # Since these FMAs can increase the performance of many numerical algorithms,
@@ -294,6 +295,59 @@ function Trixi.solve_b_embedded(consistency_order, num_eig_vals, num_stage_evals
 
     return b_embedded, dt
 end
+
+function Trixi.compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+
+    A = zeros(num_stage_evals-1, num_stage_evals)
+    rhs = [1, 1/2, embedded_monomial_coeffs...]
+
+    # Define the variables
+    b_embedded = Variable(num_stage_evals)  # the unknown coefficients we want to solve for
+
+    # Initialize A matrix for the constraints
+    A[1, :] .= 1  # sum(b) = 1 constraint row
+    for i in 1:num_stage_evals-1
+        A[2, i+1] = c[num_stages - num_stage_evals + i]  # second order constraint: dot(b, c) = 1/2
+    end
+
+    # Fill the A matrix with other constraints based on monomial coefficients
+    for i in 3:(num_stage_evals-1)
+        for j in i+1:num_stage_evals
+            A[i, j] = c[num_stages - num_stage_evals + j - 2]
+            for k in 1:i-2
+                A[i, j] *= a_unknown[j - 1 - k]  # recursive dependence on `a_unknown`
+            end
+        end
+    end
+
+    println("A matrix of finding b")
+    display(A)
+
+    # Set up weights to prioritize the first two constraints
+    weights = [2, 2, ones(num_stage_evals-3)...]  # Heavier weight for the first two rows
+    weighted_residual = weights .* (A * b_embedded - rhs)  # Element-wise multiplication for weighting
+
+    # Set up the objective to minimize the weighted norm of the difference
+    problem = minimize(sumsquares(weighted_residual), [b_embedded >= 0])
+
+    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-7,
+                                "abstol_inacc" => 5e-6,
+                                "reltol_inacc" => 5e-7,
+                                "nitref" => 9,
+                                "maxit" => 1000000,
+                                "verbose" => 3); silent_solver = true)
+
+
+    ot = problem.optval
+    println("Optimal value of the objective function: ", ot)
+    return evaluate(b_embedded)
+end
 end # @muladd
 
 end # module TrixiConvexECOSExt

ext/TrixiNLsolveExt.jl

@@ -174,56 +174,6 @@ function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, num_stage_
 
     error("Maximum number of iterations ($max_iter) reached. Cannot find valid sets of coefficients.")
 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)
-
-    for i in 1:num_a_unknown
-        a[num_stages - i + 1] = a_unknown[num_a_unknown - i + 1]
-    end
-
-    # Define the objective_function
-    function embedded_scheme_objective_function!(b_eq, x)
-        return EmbeddedPairedExplicitRK3_butcher_tableau_objective_function!(b_eq, x,
-                                                                     num_stages,
-                                                                     num_stage_evals,
-                                                                     embedded_monomial_coeffs,
-                                                                     c, a)
-    end
-
-    # RealT is determined as the type of the first element in monomial_coeffs to ensure type consistency
-    RealT = typeof(embedded_monomial_coeffs[1])
-
-    # To ensure consistency and reproducibility of results across runs, we use 
-    # a seeded random initial guess.
-    rng = StableRNG(55555)
-
-    # Due to the nature of the nonlinear solver, different initial guesses can lead to 
-    # small numerical differences in the solution.
-
-    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)
-
-        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
-
-    end
-end
 end # @muladd
 
 end # module TrixiNLsolveExt