From b2358e5a3cddd20fbbece88e68a7ca05fae3a670 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 25 Oct 2024 20:52:37 +0200 Subject: [PATCH 1/3] Finally some positive b --- .../elixir_advection_embedded_perk3.jl | 7 ++- ext/TrixiConvexECOSExt.jl | 59 ++++++++++++++++++- ext/TrixiNLsolveExt.jl | 50 ---------------- .../methods_embedded_PERK3.jl | 21 ++++++- .../paired_explicit_runge_kutta.jl | 2 +- 5 files changed, 81 insertions(+), 58 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl index 0221817e1cf..61e202f6984 100644 --- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl +++ b/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(10, 7, 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). @@ -79,6 +79,9 @@ sol = Trixi.solve(ode, ode_algorithm, # Print the timer summary summary_callback() +ode_algorithm.b + + # 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 @@ -98,4 +101,4 @@ 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) \ No newline at end of file +println("ratio = ", ode_algorithm.dt_opt_b / ode_algorithm.dt_opt_a * 100) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 0a5b4c44c5d..518eaa6c14f 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/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 @@ -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, 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, @@ -160,6 +160,59 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, return gamma_opt, 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 diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl index 319134e3700..ded4e7e281a 100644 --- a/ext/TrixiNLsolveExt.jl +++ b/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 diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl index d859014d712..f29c1cffee3 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl @@ -7,8 +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) @@ -45,6 +45,19 @@ function compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomia b_embedded = A \ rhs return b_embedded end +=# + +# 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-1)..., + b_unknown[2:end]..., + 0 + ] + return b +end # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 # using a list of eigenvalues @@ -94,13 +107,17 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, println("dt_opt_b = ", dt_opt_b) b = compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c) + + b_full = construct_b_vector(b, num_stages - 1, num_stage_evals - 1) + println("dot(b, c) = ", dot(b_full, c)) + println("sum(b) = ", sum(b_full)) println("b: ", b) println("Embedded monomial after normalization: ", embedded_monomial_coeffs) println("dt_opt_a = ", dt_opt_a) - #error("b found.") + error("b found.") end # Fill A-matrix in P-ERK style a_matrix = zeros(num_stage_evals - 2, 2) diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl index 52d2d307db7..af396cb3ae7 100644 --- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -16,5 +16,5 @@ include("polynomial_optimizer.jl") # such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package # extension or by the NLsolve-specific code loaded by Requires.jl function solve_a_butcher_coeffs_unknown! end -function solve_b_butcher_coeffs_unknown! end +function compute_b_embedded_coeffs end end # @muladd From 2af55dd94c40dad4f5654cc64ee14fc166de167c Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 25 Oct 2024 21:25:12 +0200 Subject: [PATCH 2/3] change with the integrator --- .../elixir_advection_embedded_perk3.jl | 4 +- .../methods_embedded_PERK3.jl | 39 ++++++++++--------- 2 files changed, 23 insertions(+), 20 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl index 61e202f6984..0f90ea22199 100644 --- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl +++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl @@ -56,12 +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(10, 7, tspan, semi) +ode_algorithm = Trixi.EmbeddedPairedRK3(10, 9, 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 +stepsize_callback = StepsizeCallback() # I've tried using cfl of 1.5 and the error is very similar. callbacks = CallbackSet(summary_callback, diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl index f29c1cffee3..1f145fb76cb 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl @@ -117,7 +117,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, println("dt_opt_a = ", dt_opt_a) - error("b found.") + #error("b found.") end # Fill A-matrix in P-ERK style a_matrix = zeros(num_stage_evals - 2, 2) @@ -370,7 +370,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator) end # Higher stages where the weight of b in the butcher tableau is zero - for stage in 2:(alg.num_stages - alg.num_stage_evals + 1) + for stage in 2:(alg.num_stages - alg.num_stage_evals) # Construct current state @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + alg.c[stage] * integrator.k1[i] @@ -384,25 +384,28 @@ function step!(integrator::EmbeddedPairedRK3Integrator) end end - # k_(s-e+2) + # #k_(s-e+1) and k_(s-e+2) # Construct current state - @threaded for i in eachindex(integrator.du) - integrator.u_tmp[i] = integrator.u[i] + - alg.c[alg.num_stages - alg.num_stage_evals + 2] * - integrator.k1[i] - end - integrator.f(integrator.du, integrator.u_tmp, prob.p, - integrator.t + - alg.c[alg.num_stages - alg.num_stage_evals + 2] * integrator.dt) + for j in 1:2 + @threaded for i in eachindex(integrator.du) + integrator.u_tmp[i] = integrator.u[i] + + alg.c[alg.num_stages - alg.num_stage_evals + 2] * + integrator.k1[i] + end - @threaded for i in eachindex(integrator.du) - integrator.k_higher[i] = integrator.du[i] * integrator.dt - end + integrator.f(integrator.du, integrator.u_tmp, prob.p, + integrator.t + + alg.c[alg.num_stages - alg.num_stage_evals + j] * integrator.dt) - @threaded for i in eachindex(integrator.u) - integrator.u[i] += integrator.k_higher[i] * - alg.b[2] + @threaded for i in eachindex(integrator.du) + integrator.k_higher[i] = integrator.du[i] * integrator.dt + end + + @threaded for i in eachindex(integrator.u) + integrator.u[i] += integrator.k_higher[i] * + alg.b[j+1] + end end # Higher stages after num_stage_evals where b is non-zero @@ -427,7 +430,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator) @threaded for i in eachindex(integrator.u) integrator.u[i] += integrator.k_higher[i] * - alg.b[stage - alg.num_stages + alg.num_stage_evals] + alg.b[stage - alg.num_stages + alg.num_stage_evals+1] end end From 6360b8d66e9163d55534eef1e3bc14e27e7f8792 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 25 Oct 2024 21:25:48 +0200 Subject: [PATCH 3/3] remove some comments and unneeded file --- .../elixir_advection_embedded_perk3.jl | 25 --- ext/TrixiConvexClarabelExt.jl | 173 ------------------ 2 files changed, 198 deletions(-) delete mode 100644 ext/TrixiConvexClarabelExt.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl index 0f90ea22199..dc1ab1c3516 100644 --- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl +++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl @@ -62,7 +62,6 @@ ode_algorithm = Trixi.EmbeddedPairedRK3(10, 9, tspan, semi) # b values in the Butcher tableau of the ODE algorithm). #cfl_number = Trixi.calculate_cfl(ode_algorithm, ode) stepsize_callback = StepsizeCallback() -# I've tried using cfl of 1.5 and the error is very similar. callbacks = CallbackSet(summary_callback, alive_callback, @@ -78,27 +77,3 @@ sol = Trixi.solve(ode, ode_algorithm, # Print the timer summary summary_callback() - -ode_algorithm.b - - -# 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) diff --git a/ext/TrixiConvexClarabelExt.jl b/ext/TrixiConvexClarabelExt.jl deleted file mode 100644 index bedc05031f7..00000000000 --- a/ext/TrixiConvexClarabelExt.jl +++ /dev/null @@ -1,173 +0,0 @@ -# Package extension for adding Convex-based features to Trixi.jl -module TrixiConvexClarabelExt - -# Required for coefficient optimization in P-ERK scheme integrators -if isdefined(Base, :get_extension) - using Convex: MOI, solve!, Variable, minimize, evaluate, dot - using Clarabel: 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 ..Clarabel: Optimizer -end - -# Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, @muladd - -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# Compute new version of stability polynomial of the embedded scheme for paired explicit Runge-Kutta -# up to specified consistency order(p = 2), including contributions from free coefficients for higher -# orders, and return the maximum absolute value -function embedded_scheme_stability_polynomials!(pnoms, - num_stages_embedded, - num_stage_evals_embedded, - normalized_powered_eigvals_scaled, - a, b, c) - - # Construct a full b coefficient vector #TODO: is there a way to not do this and just use b directly? - b_coeff = [ - 1 - sum(b), - zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., - b..., - 0 - ] - num_eig_vals = length(pnoms) - - # Initialize with 1 + z - for i in 1:num_eig_vals - pnoms[i] = 1.0 + normalized_powered_eigvals_scaled[i, 1] - end - # z^2: b^T * c - #pnoms += dot(b, c) * normalized_powered_eigvals_scaled[:, 2] - pnoms += 0.5 * normalized_powered_eigvals_scaled[:, 2] - - # Contribution from free coefficients - for i in 3:num_stage_evals_embedded - sum = 0.0 - for j in (i + num_stages_embedded - num_stage_evals_embedded):num_stages_embedded - prod = 1.0 - for k in (3 + j - i):j - prod *= a[k] - end - sum += prod * b_coeff[j] * c[j - i + 2] - end - pnoms += sum * normalized_powered_eigvals_scaled[:, i] - end - - # For optimization only the maximum is relevant - return maximum(abs(pnoms)) -end - -function solve_b_butcher_coeffs_unknown_new(num_eig_vals, eig_vals, - num_stages, num_stage_evals, - num_stages_embedded, - num_stage_evals_embedded, - a_unknown, c, dtmax, dteps) - dtmin = 0.0 - dt = -1.0 - abs_p = -1.0 - - # 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 - - # Construct stability polynomial for each eigenvalue - pnoms = ones(Complex{Float64}, num_eig_vals, 1) - - # There are e - 2 free variables for the stability polynomial of the embedded scheme - b = Variable(num_stage_evals - 2) - - normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals-1) - - for j in 1:num_stage_evals-1 - #fac_j = factorial(j) - for i in 1:num_eig_vals - #normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j - # Try first without factorial normalization - normalized_powered_eigvals[i, j] = eig_vals[i]^j - end - end - - normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals) - - # 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-1 - 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 - - # Second-order constraint - # Since c[1] is always 0 we can ignore the contribution of b[1] and only account for the ones from other non-zero entries of b - constraints = [b >= 0, - 2 * dot(b, c[(num_stages - num_stage_evals + 2):(num_stages - 1)]) == 1.0] - - # Use last optimal values for b in (potentially) next iteration - problem = minimize(embedded_scheme_stability_polynomials!(pnoms, - num_stages_embedded, - num_stage_evals_embedded, - normalized_powered_eigvals_scaled, - a, b, c), constraints) - - #= - solve!(problem, - # Parameters taken from default values for EiCOS - MOI.OptimizerWithAttributes(Optimizer, "b" => 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) - =# - - - solve!(problem, - MOI.OptimizerWithAttributes(Optimizer, - "tol_feas" => 1e-4); - silent_solver = true) - - print("dt: ", dt, "\n") - print("Problem value: ", problem.optval, "\n") - abs_p = problem.optval - - if abs_p < 1 - dtmin = dt - else - dtmax = dt - end - end - - b_opt = evaluate(b) - - # Catch case S = 3 (only one opt. variable) - if isa(b_opt, Number) - b_opt = [b_opt] - end - - return b_opt, dt -end -end # @muladd - -end # module TrixiConvexClarabelExt