Implement simple CFL condition

examples/fluid/rectangular_tank_2d.jl

@@ -68,18 +68,10 @@ ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+stepsize_callback = StepsizeCallback(cfl=0.8)
 
-callbacks = CallbackSet(info_callback, saving_callback)
-
-# Use a Runge-Kutta method with automatic (error based) time step size control.
-# Limiting of the maximum stepsize is necessary to prevent crashing.
-# When particles are approaching a wall in a uniform way, they can be advanced
-# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
-# Sometimes, the method fails to do so because forces become extremely large when
-# fluid particles are very close to boundary particles, and the time integration method
-# interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
-            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
-            reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
-            dtmax=1e-2, # Limit stepsize to prevent crashing
+callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # This is overwritten by the stepsize callback
             save_everystep=false, callback=callbacks);

src/TrixiParticles.jl

@@ -11,7 +11,7 @@ using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
-                 get_tmp_cache
+                 get_tmp_cache, set_proposed_dt!
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
@@ -37,7 +37,8 @@ export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
        BoundarySPHSystem
-export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback
+export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
+       StepsizeCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,

src/callbacks/callbacks.jl

@@ -8,3 +8,4 @@ end
 include("info.jl")
 include("solution_saving.jl")
 include("density_reinit.jl")
+include("stepsize.jl")

src/callbacks/stepsize.jl

@@ -0,0 +1,73 @@
+struct StepsizeCallback{ISCONSTANT, ELTYPE}
+    cfl_number::ELTYPE
+end
+
+@inline is_constant(::StepsizeCallback{ISCONSTANT}) where {ISCONSTANT} = ISCONSTANT
+
+function StepsizeCallback(; cfl::Real)
+    # TODO adapt for non-constant CFL conditions
+    is_constant = true
+    stepsize_callback = StepsizeCallback{is_constant, typeof(cfl)}(cfl)
+
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(stepsize_callback, stepsize_callback,
+                            save_positions=(false, false),
+                            initialize=initialize_stepsize_callback)
+end
+
+function initialize_stepsize_callback(discrete_callback, u, t, integrator)
+    stepsize_callback = discrete_callback.affect!
+
+    stepsize_callback(integrator)
+end
+
+# `condition`
+function (stepsize_callback::StepsizeCallback)(u, t, integrator)
+    # Only apply the callback when the stepsize is not constant and the time integrator
+    # is not adaptive.
+    return !is_constant(stepsize_callback) && !integrator.opts.adaptive
+end
+
+# `affect!`
+function (stepsize_callback::StepsizeCallback)(integrator)
+    (; cfl_number) = stepsize_callback
+
+    v_ode, u_ode = integrator.u.x
+    semi = integrator.p
+
+    dt = @trixi_timeit timer() "calculate dt" calculate_dt(v_ode, u_ode, cfl_number, semi)
+
+    set_proposed_dt!(integrator, dt)
+    integrator.opts.dtmax = dt
+    integrator.dtcache = dt
+
+    # Tell OrdinaryDiffEq that `u` has not been modified
+    u_modified!(integrator, false)
+
+    return stepsize_callback
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:StepsizeCallback})
+    @nospecialize cb # reduce precompilation time
+
+    stepsize_callback = cb.affect!
+    print(io, "StepsizeCallback(is_constant=", is_constant(stepsize_callback),
+          "cfl_number=", stepsize_callback.cfl_number, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, <:StepsizeCallback})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        stepsize_callback = cb.affect!
+
+        setup = [
+            "is constant" => is_constant(stepsize_callback),
+            "CFL number" => stepsize_callback.cfl_number,
+        ]
+        summary_box(io, "StepsizeCallback", setup)
+    end
+end

src/general/semidiscretization.jl

@@ -285,6 +285,12 @@ end
                     (StaticInt(v_nvariables(system)), n_moving_particles(system)))
 end
 
+function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization)
+    (; systems) = semi
+
+    return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems)
+end
+
 function drift!(du_ode, v_ode, u_ode, semi, t)
     @trixi_timeit timer() "drift!" begin
         @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode)

src/schemes/boundary/system.jl

@@ -264,3 +264,7 @@ end
 function viscosity_model(system::BoundarySPHSystem)
     return system.boundary_model.viscosity
 end
+
+function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem)
+    return Inf
+end

src/schemes/fluid/fluid.jl

@@ -15,6 +15,46 @@ end
 
 @inline viscosity_model(system::FluidSystem) = system.viscosity
 
+function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem)
+    (; smoothing_length, state_equation, viscosity, acceleration) = system
+    (; sound_speed) = state_equation
+    # TODO Specific to artificial viscosity
+    (; alpha) = viscosity
+
+    # TODO This is based on:
+    # M. Antuono, A. Colagrossi, S. Marrone.
+    # "Numerical Diffusive Terms in Weakly-Compressible SPH Schemes."
+    # In: Computer Physics Communications 183, no. 12 (2012), pages 2570-80.
+    # https://doi.org/10.1016/j.cpc.2012.07.006
+    # But in the docs for the artificial viscosity, it says that we have to divide by rho.
+    kinematic_viscosity = alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4)
+
+    dt_viscosity = 0.125 * smoothing_length^2 / kinematic_viscosity
+
+    # TODO Adami et al. (2012) just use the gravity here, but Antuono et al. (2012)
+    # are using a per-particle acceleration. Is that supposed to be the previous RHS?
+    dt_acceleration = 0.25 * sqrt(smoothing_length / norm(acceleration))
+
+    # TODO Everyone seems to be doing this differently.
+    # Sun et al. (2017) only use h / c (because c depends on v_max as c >= 10 v_max).
+    # Adami et al. (2012) use h / (c + v_max) with a fixed CFL of 0.25.
+    # Antuono et al. (2012) use h / (c + v_max + h * pi_max), where pi is the viscosity coefficient.
+    # Antuono et al. (2015) use h / (c + h * pi_max).
+    #
+    # P. N. Sun, A. Colagrossi, S. Marrone, A. M. Zhang.
+    # "The δplus-SPH Model: Simple Procedures for a Further Improvement of the SPH Scheme."
+    # In: Computer Methods in Applied Mechanics and Engineering 315 (2017), pages 25–49.
+    # https://doi.org/10.1016/j.cma.2016.10.028.
+    #
+    # M. Antuono, S. Marrone, A. Colagrossi, B. Bouscasse.
+    # "Energy Balance in the δ-SPH Scheme."
+    # In: Computer Methods in Applied Mechanics and Engineering 289 (2015), pages 209–26.
+    # https://doi.org/10.1016/j.cma.2015.02.004.
+    dt_sound_speed = cfl_number * smoothing_length / sound_speed
+
+    return min(dt_viscosity, dt_acceleration, dt_sound_speed)
+end
+
 include("pressure_acceleration.jl")
 include("viscosity.jl")
 include("weakly_compressible_sph/weakly_compressible_sph.jl")