Time-dependent CFL (#2248)

* Variable CFL

* Example

* fix text

* fix

* CFL for coupled semi

* nnews

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/callbacks_step/stepsize.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* shorten code

* coupled glm

* fmt

* Update test/test_structured_2d.jl

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: DanielDoehring

NEWS.md

@@ -14,6 +14,11 @@ for human readability.
 - `LobattoLegendreBasis` and related datastructures made fully floating-type general,
   enabling calculations with higher than double (`Float64`) precision ([#2128])
 - In 2D, quadratic elements, i.e., 8-node (quadratic) quadrilaterals are now supported in standard Abaqus `inp` format ([#2217])
+- The `cfl` value supplied in the `StepsizeCallback` and `GlmStepsizeCallback` can now be a function of simulation 
+  time `t` to enable e.g. a ramp-up of the CFL value.
+  This is useful for simulations that are initialized with an "unphysical" initial condition, but do not permit the usage of
+  adaptive, error-based timestepping.
+  Examples for this are simulations involving the MHD equations which require in general the `GlmStepsizeCallback` ([#2248])
 
 #### Changed
 

examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl

@@ -0,0 +1,141 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+equations = IdealGlmMhdEquations2D(1.4)
+
+"""
+    initial_condition_rotor(x, t, equations::IdealGlmMhdEquations2D)
+
+The classical MHD rotor test case. Here, the setup is taken from
+- Dominik Derigs, Gregor J. Gassner, Stefanie Walch & Andrew R. Winters (2018)
+  Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics
+  [doi: 10.1365/s13291-018-0178-9](https://doi.org/10.1365/s13291-018-0178-9)
+"""
+function initial_condition_rotor(x, t, equations::IdealGlmMhdEquations2D)
+    # setup taken from Derigs et al. DMV article (2018)
+    # domain must be [0, 1] x [0, 1], γ = 1.4
+    dx = x[1] - 0.5
+    dy = x[2] - 0.5
+    r = sqrt(dx^2 + dy^2)
+    f = (0.115 - r) / 0.015
+    if r <= 0.1
+        rho = 10.0
+        v1 = -20.0 * dy
+        v2 = 20.0 * dx
+    elseif r >= 0.115
+        rho = 1.0
+        v1 = 0.0
+        v2 = 0.0
+    else
+        rho = 1.0 + 9.0 * f
+        v1 = -20.0 * f * dy
+        v2 = 20.0 * f * dx
+    end
+    v3 = 0.0
+    p = 1.0
+    B1 = 5.0 / sqrt(4.0 * pi)
+    B2 = 0.0
+    B3 = 0.0
+    psi = 0.0
+    return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+initial_condition = initial_condition_rotor
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Affine type mapping to take the [-1,1]^2 domain from the mesh file
+# and put it onto the rotor domain [0,1]^2 and then warp it with a mapping
+# as described in https://arxiv.org/abs/2012.12040
+function mapping_twist(xi, eta)
+    y = 0.5 * (eta + 1.0) +
+        0.05 * cos(1.5 * pi * (2.0 * xi - 1.0)) * cos(0.5 * pi * (2.0 * eta - 1.0))
+    x = 0.5 * (xi + 1.0) + 0.05 * cos(0.5 * pi * (2.0 * xi - 1.0)) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+                           joinpath(@__DIR__, "square_unstructured_2.inp"))
+
+mesh = P4estMesh{2}(mesh_file,
+                    polydeg = 4,
+                    mapping = mapping_twist,
+                    initial_refinement_level = 1)
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:all => boundary_condition)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.15)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorLöhner(semi,
+                                variable = density_pressure)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 1,
+                                      med_level = 3, med_threshold = 0.05,
+                                      max_level = 5, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 3,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+# increase the CFL number linearly from cfl_0() at time 0 
+# to cfl_t_ramp() at time t = t_ramp(), keep it constant afterward
+cfl_0() = 0.5
+cfl_t_ramp() = 1.2
+t_ramp() = 0.1
+cfl(t) = min(cfl_0() + (cfl_t_ramp() - cfl_0()) / t_ramp() * t, cfl_t_ramp())
+
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        amr_callback,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode,
+            CarpenterKennedy2N54(thread = OrdinaryDiffEq.True(),
+                                 williamson_condition = false),
+            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

src/callbacks_step/glm_speed.jl

@@ -11,8 +11,11 @@
 Update the divergence cleaning wave speed `c_h` according to the time step
 computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations, the multi-component
 GLM-MHD equations, and the multi-ion GLM-MHD equations.
-The `cfl` number should be set to the same value as for the time step size calculation. The
-`glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
+The `cfl` number should be set to the same value as for the time step size calculation. 
+As for standard [`StepsizeCallback`](@ref) `cfl` can be either set to a `Real` number or
+a function of time `t` returning a `Real` number.
+
+The `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
 
@@ -22,9 +25,9 @@ semidiscretization, the divergence cleaning should be applied. See also
 Note: `SemidiscretizationCoupled` and all related features are considered experimental and
 may change at any time.
 """
-struct GlmSpeedCallback{RealT <: Real}
+struct GlmSpeedCallback{RealT <: Real, CflType}
     glm_scale::RealT
-    cfl::RealT
+    cfl::CflType
     semi_indices::Vector{Int}
 end
 
@@ -58,7 +61,9 @@ end
 function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = Int[])
     @assert 0<=glm_scale<=1 "glm_scale must be between 0 and 1"
 
-    glm_speed_callback = GlmSpeedCallback(glm_scale, cfl, semi_indices)
+    glm_speed_callback = GlmSpeedCallback{typeof(glm_scale), typeof(cfl)}(glm_scale,
+                                                                          cfl,
+                                                                          semi_indices)
 
     DiscreteCallback(glm_speed_callback, glm_speed_callback, # the first one is the condition, the second the affect!
                      save_positions = (false, false),
@@ -75,13 +80,17 @@ function (glm_speed_callback::GlmSpeedCallback)(u, t, integrator)
     return true
 end
 
-function update_cleaning_speed!(semi, glm_speed_callback, dt)
+function update_cleaning_speed!(semi, glm_speed_callback, dt, t)
     @unpack glm_scale, cfl = glm_speed_callback
 
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
     # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
-    c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+    if cfl isa Real # Case for constant CFL
+        c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+    else # Variable CFL
+        c_h_deltat = calc_dt_for_cleaning_speed(cfl(t), mesh, equations, solver, cache)
+    end
 
     # c_h is proportional to its own time step divided by the complete MHD time step
     # We use @reset here since the equations are immutable (to work on GPUs etc.).
@@ -99,7 +108,7 @@ end
 
     # Call the appropriate update function (this indirection allows to specialize on,
     # e.g., the semidiscretization type)
-    update_cleaning_speed!(semi, glm_speed_callback, dt)
+    update_cleaning_speed!(semi, glm_speed_callback, dt, integrator.t)
 
     # avoid re-evaluating possible FSAL stages
     u_modified!(integrator, false)

src/callbacks_step/stepsize.jl

@@ -10,9 +10,12 @@
 
 Set the time step size according to a CFL condition with CFL number `cfl`
 if the time integration method isn't adaptive itself.
+
+The supplied keyword argument `cfl` must be either a `Real` number or
+a function of time `t` returning a `Real` number.
 """
-mutable struct StepsizeCallback{RealT}
-    cfl_number::RealT
+mutable struct StepsizeCallback{CflType}
+    cfl_number::CflType
 end
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:StepsizeCallback})
@@ -39,8 +42,8 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
-function StepsizeCallback(; cfl::Real = 1.0)
-    stepsize_callback = StepsizeCallback(cfl)
+function StepsizeCallback(; cfl = 1.0)
+    stepsize_callback = StepsizeCallback{typeof(cfl)}(cfl)
 
     DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect!
                      save_positions = (false, false),
@@ -80,16 +83,6 @@ end
     return nothing
 end
 
-# General case for a single semidiscretization
-function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization)
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
-    u = wrap_array(u_ode, mesh, equations, solver, cache)
-
-    dt = cfl_number * max_dt(u, t, mesh,
-                have_constant_speed(equations), equations,
-                solver, cache)
-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
@@ -103,11 +96,28 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond
     u_ode = ode.u0
     t = first(ode.tspan)
     semi = ode.p
+
+    dt = calculate_dt(u_ode, t, cfl_number, semi)
+end
+
+# General case for a single (i.e., non-coupled) semidiscretization
+# Case for constant `cfl_number`.
+function calculate_dt(u_ode, t, cfl_number::Real, semi::AbstractSemidiscretization)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-    return cfl_number *
-           max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache)
+    dt = cfl_number * max_dt(u, t, mesh,
+                have_constant_speed(equations), equations,
+                solver, cache)
+end
+# Case for `cfl_number` as a function of time `t`.
+function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    dt = cfl_number(t) * max_dt(u, t, mesh,
+                have_constant_speed(equations), equations,
+                solver, cache)
 end
 
 include("stepsize_dg1d.jl")

src/semidiscretization/semidiscretization_coupled.jl

@@ -351,17 +351,26 @@ end
 ################################################################################
 
 # In case of coupled system, use minimum timestep over all systems
-function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled)
+# Case for constant `cfl_number`.
+function calculate_dt(u_ode, t, cfl_number::Real, semi::SemidiscretizationCoupled)
     dt = minimum(eachsystem(semi)) do i
         u_ode_slice = get_system_u_ode(u_ode, i, semi)
         calculate_dt(u_ode_slice, t, cfl_number, semi.semis[i])
     end
 
     return dt
 end
+# Case for `cfl_number` as a function of time `t`.
+function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled)
+    cfl_number_ = cfl_number(t)
+    dt = minimum(eachsystem(semi)) do i
+        u_ode_slice = get_system_u_ode(u_ode, i, semi)
+        calculate_dt(u_ode_slice, t, cfl_number_, semi.semis[i])
+    end
+end
 
 function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
-                                glm_speed_callback, dt)
+                                glm_speed_callback, dt, t)
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
     if length(semi_indices) == 0
@@ -376,12 +385,19 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
         end
     end
 
+    if cfl isa Real # Case for constant CFL
+        cfl_number = cfl
+    else # Variable CFL
+        cfl_number = cfl(t)
+    end
+
     for semi_index in semi_indices
         semi = semi_coupled.semis[semi_index]
         mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
-        c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+        c_h_deltat = calc_dt_for_cleaning_speed(cfl_number,
+                                                mesh, equations, solver, cache)
 
         # c_h is proportional to its own time step divided by the complete MHD time step
         # We use @reset here since the equations are immutable (to work on GPUs etc.).

test/test_p4est_2d.jl

@@ -666,6 +666,41 @@ end
     end
 end
 
+@trixi_testset "elixir_mhd_rotor_cfl_ramp.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor_cfl_ramp.jl"),
+                        l2=[
+                            0.45519051169507474,
+                            0.8917985468745363,
+                            0.8324681609772325,
+                            0.0,
+                            0.9801426190285389,
+                            0.10476233464125001,
+                            0.15551270692826116,
+                            0.0,
+                            2.0201603821472296e-5
+                        ],
+                        linf=[
+                            10.196786739705292,
+                            18.267539012179128,
+                            10.046104290498878,
+                            0.0,
+                            19.668302849210974,
+                            1.395022093528294,
+                            1.8717844606331189,
+                            0.0,
+                            0.001651262488701531
+                        ],
+                        tspan=(0.0, 0.02))
+    # 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)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_linearizedeuler_gaussian_source.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_linearizedeuler_gaussian_source.jl"),