diff --git a/NEWS.md b/NEWS.md index 40eb3db390..dc49855738 100644 --- a/NEWS.md +++ b/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 diff --git a/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl b/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl new file mode 100644 index 0000000000..d33cbacc9b --- /dev/null +++ b/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 diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 075e084949..285506186e 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/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) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 0a95cb35ad..f6d7a8fc2a 100644 --- a/src/callbacks_step/stepsize.jl +++ b/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") diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 745a8d3f6f..f4130b3fa5 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -351,7 +351,8 @@ 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]) @@ -359,9 +360,17 @@ function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled) 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.). diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 5d9f92a7c0..51e9bc015b 100644 --- a/test/test_p4est_2d.jl +++ b/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"), diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index f78f127f43..ea3584db92 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -111,6 +111,31 @@ end end end +@trixi_testset "elixir_advection_meshview.jl with time-dependent CFL" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_meshview.jl"), + l2=[ + 8.311947673083206e-6, + 8.311947673068427e-6 + ], + linf=[ + 6.627000273318195e-5, + 6.62700027264096e-5 + ], + coverage_override=(maxiters = 10^5,), + stepsize_callback=StepsizeCallback(cfl = x -> 1.6)) + + @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin + # 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 +end + @trixi_testset "elixir_advection_extended.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), l2=[4.220397559713772e-6],