Updated cfl number calculation

examples/tree_1d_dgsem/elixir_advection_perk2.jl

@@ -39,16 +39,22 @@ summary_callback = SummaryCallback()
 analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
-stepsize_callback = StepsizeCallback(cfl = 2.5)
-
 alive_callback = AliveCallback(alive_interval = analysis_interval)
 
 save_solution = SaveSolutionCallback(dt = 0.1,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
+# 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)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+# For PERK schemes, the CFL number is calculated from the optimal time step of the scheme.
+stepsize_callback = StepsizeCallback(ode, ode_algorithm)
+
 # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
 callbacks = CallbackSet(summary_callback,
                         alive_callback,
@@ -58,12 +64,6 @@ callbacks = CallbackSet(summary_callback,
 
 ###############################################################################
 # run the simulation
-
-# 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)
-
 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);

src/callbacks_step/stepsize.jl

@@ -47,6 +47,22 @@ function StepsizeCallback(; cfl::Real = 1.0)
                      initialize = initialize!)
 end
 
+# For Paired-Explicit Runge-Kutta methods, the CFL number is calculated based on the optimal timestep
+function StepsizeCallback(ode, ode_algorithm::AbstractPairedExplicitRKSingle)
+    # TODO: Loop over all the time step and choose the minimum CFL number? from the loop???
+    t = first(ode.tspan)
+    u_ode = ode.u0
+    semi = ode.p
+    dt_opt = ode_algorithm.dt_opt
+    cfl = calculate_cfl(u_ode, t, dt_opt, semi)
+
+    stepsize_callback = StepsizeCallback(cfl)
+
+    DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect!
+                     save_positions = (false, false),
+                     initialize = initialize!)
+end
+
 function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t,
                      integrator) where {Condition, Affect! <: StepsizeCallback}
     cb.affect!(integrator)
@@ -90,6 +106,19 @@ function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization)
                 solver, cache)
 end
 
+# For Paired Explicit Runge-Kutta methods, use the CFL number calculated from the optimal timestep of the
+# scheme.
+function calculate_cfl(u_ode, t, dt_opt, semi::AbstractSemidiscretization)
+    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, t, mesh,
+                        have_constant_speed(equations), equations,
+                        solver, cache)
+
+    return cfl_number
+end
+
 # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own
 # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have
 # an integrator at this stage but only the ODE, this method will be used there. It's called in

src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl

@@ -60,15 +60,15 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
                                                           eig_vals; verbose)
     monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
                                           num_stages)
-
+    
     num_monomial_coeffs = length(monomial_coeffs)
     @assert num_monomial_coeffs == coeffs_max
     A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs)
 
     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
@@ -77,6 +77,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages,
                                                    base_path_monomial_coeffs::AbstractString,
                                                    bS, cS)
 
+    #TODO: If the user has a specific set of monomial coefficients, they must also have already obtained dt_opt
+    #TODO: where did I get this monomial in unit test....
     # c Vector form Butcher Tableau (defines timestep per stage)
     c = zeros(num_stages)
     for k in 2:num_stages
@@ -107,7 +109,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 +120,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.
     - `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,10 +147,11 @@ 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,
+function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt,
                            bS = 1.0, cS = 0.5)
     a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages,
                                                             base_path_monomial_coeffs,
@@ -171,12 +175,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,
+    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

test/test_tree_1d_advection.jl

@@ -82,10 +82,11 @@ end
     end
 end
 
+#TODO: replace l2 and linf errors with the values in CI
 @trixi_testset "elixir_advection_perk2.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"),
-                        l2=[0.011288030389423475],
-                        linf=[0.01596735472556976])
+                        l2=[1.45247275e-02],
+                        linf=[2.05428436e-02])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -96,6 +97,7 @@ end
     end
 end
 
+# TODO: edit these values as well
 # Testing the second-order paired explicit Runge-Kutta (PERK) method without stepsize callback
 @trixi_testset "elixir_advection_perk2.jl(fixed time step)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"),

test/test_unit.jl

@@ -1670,8 +1670,8 @@ end
     path_coeff_file = mktempdir()
     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