From 1672c2a7af12a866126bc2894e6f43c5afa9e9ae Mon Sep 17 00:00:00 2001 From: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> Date: Sat, 28 Sep 2024 11:59:48 +0700 Subject: [PATCH] Implement a function to calculate the optimal CFL number for the PERK2 integrator and apply the necessary related updates (#2083) * Update stepsize_callback CFL to 3.0 and add calculate_cfl * delete unnecessary line * update calculation of cfl number * update tests * update cfl number calculation for PERK (put this in the constructor) * revert an unnecessary change on TrixiConvexECOS * revert another unnecessary change i made in stepsize.jl * update test values but need to be changed again according to CI workflow * revert changes made in test_tree_1d_advaction.jl since we shouldn't alter the existing example * revert changes made in stepsize.jl * bring back the old example * add new example and create a function that simply return cfl number calculated from dt_opt +fmt * revert unnecessary change from formatting I made in stepsize.jl * revert unnecessary fmt changes * revert another change in test_unit.jl * correct a constructor + add and delete some comments. * add a test for optimal cfl number of PERK2 * correct the values of test to math thosein CI workflow * Update test/test_unit.jl Co-authored-by: Daniel Doehring * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring * use amr with the current example * fix test values + fmt * Update test/test_tree_1d_advection.jl --------- Co-authored-by: Daniel Doehring Co-authored-by: Michael Schlottke-Lakemper --- .../elixir_advection_perk2_optimal_cfl.jl | 84 +++++++++++++++++++ .../methods_PERK2.jl | 41 +++++++-- test/test_tree_1d_advection.jl | 19 +++++ test/test_unit.jl | 2 +- 4 files changed, 136 insertions(+), 10 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl new file mode 100644 index 0000000000..2b9602ab4a --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl @@ -0,0 +1,84 @@ + +using Convex, ECOS +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = 1.0 +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = -1.0 # minimum coordinate +coordinates_max = 1.0 # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 20.0 +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(alive_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.1, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 4, + med_level = 5, med_threshold = 0.1, + max_level = 6, max_threshold = 0.6) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +# Construct second order paired explicit Runge-Kutta method with 6 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.PairedExplicitRK2(6, tspan, semi) + +# For Paired Explicit Runge-Kutta methods, we receive an optimized timestep for a given reference semidiscretization. +# To allow for e.g. adaptivity, we reverse-engineer the corresponding CFL number to make it available during the simulation. +cfl_number = Trixi.calculate_cfl(ode_algorithm, ode) +stepsize_callback = StepsizeCallback(cfl = cfl_number) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, + alive_callback, + save_solution, + analysis_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation +sol = Trixi.solve(ode, ode_algorithm, + dt = 1.0, # Manual time step value, will be overwritten by the stepsize_callback when it is specified. + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 23a3ceba76..ad385e6df2 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -68,7 +68,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, a_matrix[:, 1] -= A a_matrix[:, 2] = A - return a_matrix, c + return a_matrix, c, dt_opt end # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2 @@ -76,7 +76,6 @@ end function compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs::AbstractString, bS, cS) - # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) for k in 2:num_stages @@ -107,7 +106,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -118,6 +117,7 @@ end - `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing monomial coefficients of the stability polynomial of PERK method. The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(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. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the @@ -144,16 +144,19 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle b1::Float64 bS::Float64 cS::Float64 + dt_opt::Float64 end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + dt_opt, bS = 1.0, cS = 0.5) + # If the user has the monomial coefficients, they also must have the optimal time step a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs, bS, cS) - return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) end # Constructor that calculates the coefficients with polynomial optimizer from a @@ -171,12 +174,12 @@ end function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, bS = 1.0, cS = 0.5) - a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, - eig_vals, tspan, - bS, cS; - verbose) + a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages, + eig_vals, tspan, + bS, cS; + verbose) - return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 @@ -232,6 +235,26 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, k_higher::uType end +""" + calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode) + +This function computes the CFL number once using the initial condition of the problem and the optimal timestep (`dt_opt`) from the ODE algorithm. +""" +function calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode) + t0 = first(ode.tspan) + u_ode = ode.u0 + semi = ode.p + dt_opt = ode_algorithm.dt_opt + + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + cfl_number = dt_opt / max_dt(u, t0, mesh, + have_constant_speed(equations), equations, + solver, cache) + return cfl_number +end + """ add_tstop!(integrator::PairedExplicitRK2Integrator, t) Add a time stop during the time integration process. diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 115c5f3c69..f061e2e1c3 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -123,6 +123,25 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 end end + +# Testing the second-order paired explicit Runge-Kutta (PERK) method with the optimal CFL number +@trixi_testset "elixir_advection_perk2_optimal_cfl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2_optimal_cfl.jl"), + l2=[0.0009700887119146429], + linf=[0.00137209242077041]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 + end +end end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 70e2e2ed10..5831122ffe 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1671,7 +1671,7 @@ end Trixi.download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", joinpath(path_coeff_file, "gamma_6.txt")) - ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file) + ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file, 42) # dummy optimal time step (dt_opt plays no role in determining `a_matrix`) @test isapprox(ode_algorithm.a_matrix, [0.12405417889682908 0.07594582110317093