From cde00a83fa1fef3e6d490927ecb29e1d88019c4a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 17 Sep 2024 19:20:10 +0200 Subject: [PATCH] =?UTF-8?q?Revert=20"Add=20function=20to=20calculate=20opt?= =?UTF-8?q?imal=20CFL=20number=20for=20PERK2=20integrator=20and=E2=80=A6"?= =?UTF-8?q?=20(#2082)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This reverts commit 8cfa2d3c1becc3dab75cb4d3569746ca7fa7b674. --- .../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, 10 insertions(+), 136 deletions(-) delete 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 deleted file mode 100644 index 2b9602ab4a4..00000000000 --- a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl +++ /dev/null @@ -1,84 +0,0 @@ - -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 ad385e6df24..23a3ceba76c 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, dt_opt + return a_matrix, c end # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2 @@ -76,6 +76,7 @@ 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 @@ -106,7 +107,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -117,7 +118,6 @@ 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,19 +144,16 @@ 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, dt_opt) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) end # Constructor that calculates the coefficients with polynomial optimizer from a @@ -174,12 +171,12 @@ end function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, bS = 1.0, cS = 0.5) - a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages, - eig_vals, tspan, - bS, cS; - verbose) + a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, + eig_vals, tspan, + bS, cS; + verbose) - return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 @@ -235,26 +232,6 @@ 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 f061e2e1c30..115c5f3c69c 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -123,25 +123,6 @@ 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 5831122ffe2..70e2e2ed107 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, 42) # dummy optimal time step (dt_opt plays no role in determining `a_matrix`) + ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file) @test isapprox(ode_algorithm.a_matrix, [0.12405417889682908 0.07594582110317093