From 44ea471bc5d25bb1365d22f9e47b87d14b8079ab Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 14 Oct 2024 22:46:03 +0200 Subject: [PATCH] add some comments --- ...lixir_euler_source_terms_embedded_perk3.jl | 16 ++++++- .../elixir_advection_embedded_perk3.jl | 3 +- .../methods_embedded_PERK3.jl | 44 +++++++++++-------- 3 files changed, 42 insertions(+), 21 deletions(-) diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_embedded_perk3.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_embedded_perk3.jl index c07296df888..c3eb657eaad 100644 --- a/examples/structured_1d_dgsem/elixir_euler_source_terms_embedded_perk3.jl +++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_embedded_perk3.jl @@ -44,7 +44,9 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) - +# Construct embedded order paired explicit Runge-Kutta method with 10 stages and 6 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, 6, tspan, semi) cfl_number = Trixi.calculate_cfl(ode_algorithm, ode) @@ -64,4 +66,16 @@ sol = Trixi.solve(ode, ode_algorithm, dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); summary_callback() # print the timer summary + +# 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 = [1 - sum(b_unknown), zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., b_unknown..., 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("cfl_number = ", cfl_number) \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl index 178534d58b5..9b8d02c1198 100644 --- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl +++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl @@ -62,7 +62,7 @@ ode_algorithm = Trixi.EmbeddedPairedRK3(10, 6, 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(cfl = cfl_number) # Warisa: This number is very small in contrast the other one from optimizing A +stepsize_callback = StepsizeCallback(cfl = cfl_number) # 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. callbacks = CallbackSet(summary_callback, @@ -80,6 +80,7 @@ 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 = [1 - sum(b_unknown), zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., b_unknown..., 0] 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 f850720ef84..e0731477152 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 @@ -51,7 +51,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, a_matrix[:, 2] = a_unknown - #TODO: something wrong when s = e + # Find the optimal b coeffficient from the Butcher tableau and its time step b, dt_opt_b = solve_b_butcher_coeffs_unknown(num_eig_vals, eig_vals, num_stages, num_stage_evals, num_stages - 1, # num_stages_embedded = num_stages - 1 @@ -59,44 +59,51 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, a_unknown, c, dtmax, dteps) - return a_matrix, b, c, dt_opt_b + return a_matrix, b, c, dt_opt_b # Return the optimal time step from the b coefficients for testing purposes end # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 # using provided values of coefficients a in A-matrix of Butcher tableau -function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, - base_path_a_coeffs::AbstractString; +function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, + base_path_coeffs::AbstractString; cS2) # Initialize array of c c = compute_c_coeffs(num_stages, cS2) # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) - a_coeffs_max = num_stages - 2 + coeffs_max = num_stage_evals - 2 - #TODO: Read the values from the file - b = zeros(num_stages) - - a_matrix = zeros(a_coeffs_max, 2) + a_matrix = zeros(coeffs_max, 2) a_matrix[:, 1] = c[3:end] - path_a_coeffs = joinpath(base_path_a_coeffs, - "a_" * string(num_stages) * ".txt") + b = zeros(coeffs_max) + + path_a_coeffs = joinpath(base_path_coeffs, + "a_" * string(num_stages) * "_" * string(num_stage_evals) * ".txt") + + path_b_coeffs = joinpath(base_path_coeffs, + "b_" * string(num_stages) * "_" * string(num_stage_evals) * ".txt") @assert isfile(path_a_coeffs) "Couldn't find file $path_a_coeffs" a_coeffs = readdlm(path_a_coeffs, Float64) num_a_coeffs = size(a_coeffs, 1) - @assert num_a_coeffs == a_coeffs_max + @assert num_a_coeffs == coeffs_max # Fill A-matrix in P-ERK style a_matrix[:, 1] -= a_coeffs a_matrix[:, 2] = a_coeffs + @assert isfile(path_b_coeffs) "Couldn't find file $path_b_coeffs" + b = readdlm(path_b_coeffs, Float64) + num_b_coeffs = size(b, 1) + @assert num_b_coeffs == coeffs_max + return a_matrix, b, c end @doc raw""" - EmbeddedPairedRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt; + EmbeddedPairedRK3(num_stages, base_path_coeffs::AbstractString, dt_opt; cS2 = 1.0f0) EmbeddedPairedRK3(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, cS2 = 1.0f0) @@ -105,8 +112,8 @@ end Parameters: - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method. - - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in - the Butcher tableau of the Runge Kutta method. + - `base_path_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix and a file constaining + in some coefficients in the b vector of the Butcher tableau of the Runge Kutta method. The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks. - `dt_opt` (`Float64`): Optimal time step size for the simulation setup. - `tspan`: Time span of the simulation. @@ -141,11 +148,10 @@ end # struct EmbeddedPairedRK3 # Constructor for previously computed A Coeffs function EmbeddedPairedRK3(num_stages, num_stage_evals, - base_path_a_coeffs::AbstractString, dt_opt; + base_path_coeffs::AbstractString, dt_opt; cS2 = 1.0f0) - #TODO: Update this constructor to have num_stage_evals as well and also return correct things - a_matrix, b, c = compute_EmbeddedPairedRK3_butcher_tableau(num_stages, - base_path_a_coeffs; + a_matrix, b, c = compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, + base_path_coeffs; cS2) return EmbeddedPairedRK3(num_stages, num_stage_evals, a_matrix, b, c, dt_opt)