Implement simple CFL condition

examples/fluid/dam_break_2d.jl

@@ -68,24 +68,16 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
-info_callback = InfoCallback(interval=250)
+info_callback = InfoCallback(interval=100)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
 
 use_reinit = false
 density_reinit_cb = use_reinit ? DensityReinitializationCallback(semi.systems[1], dt=0.01) :
                     nothing
+stepsize_callback = StepsizeCallback(cfl=1.1)
 
-callbacks = CallbackSet(info_callback, saving_callback, density_reinit_cb)
-
-# 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, RDPK3SpFSAL35(),
-            abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
-            reltol=1e-5, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
-            dtmax=1e-3, # 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

@@ -13,7 +13,7 @@ using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @series
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
-                 get_tmp_cache, ODESolution, ODEProblem
+                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
@@ -40,7 +40,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,116 @@
+struct StepsizeCallback{ISCONSTANT, ELTYPE}
+    cfl_number::ELTYPE
+end
+
+@inline is_constant(::StepsizeCallback{ISCONSTANT}) where {ISCONSTANT} = ISCONSTANT
+
+@doc raw"""
+    StepsizeCallback(; cfl::Real)
+
+Set the time step size according to a CFL condition if the time integration method isn't
+adaptive itself.
+
+The current implementation is using the simplest form of CFL condition, which chooses a
+time step size that is constant during the simulation.
+The step size is therefore only applied once at the beginning of the simulation.
+
+The step size ``\Delta t`` is chosen as the minimum
+```math
+    \Delta t = \min(\Delta t_\eta, \Delta t_a, \Delta t_c),
+```
+where
+```math
+    \Delta t_\eta = 0.125 \, h^2 / \eta, \quad \Delta t_a = 0.25 \sqrt{h / \lVert g \rVert},
+    \quad \Delta t_c = \text{CFL} \, h / c,
+```
+with ``\nu = \alpha h c / (2n + 4)``, where ``\alpha`` is the parameter of the viscosity
+and ``n`` is the number of dimensions.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+## References
+- 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.
+  [doi: 10.1016/j.cpc.2012.07.006](https://doi.org/10.1016/j.cpc.2012.07.006)
+- S. Adami, X. Y. Hu, N. A. Adams.
+  "A generalized wall boundary condition for smoothed particle hydrodynamics".
+  In: Journal of Computational Physics 231, 21 (2012), pages 7057--7075.
+  [doi: 10.1016/J.JCP.2012.05.005](https://doi.org/10.1016/J.JCP.2012.05.005)
+- 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.
+  [doi: 10.1016/j.cma.2016.10.028](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.
+  [doi: 10.1016/j.cma.2015.02.004](https://doi.org/10.1016/j.cma.2015.02.004)
+"""
+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" => string(is_constant(stepsize_callback)),
+            "CFL number" => stepsize_callback.cfl_number,
+        ]
+        summary_box(io, "StepsizeCallback", setup)
+    end
+end

src/general/semidiscretization.jl

@@ -296,6 +296,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

@@ -259,3 +259,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,29 @@ 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
+
+    kinematic_viscosity = kinematic_viscosity(system, viscosity)
+    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).
+    #
+    # See docstring of the callback for the references.
+    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")

src/schemes/fluid/viscosity.jl

@@ -100,6 +100,19 @@ end
     return 0.0
 end
 
+# See, e.g.,
+# 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
+function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan)
+    (; smoothing_length, state_equation) = system
+    (; sound_speed) = state_equation
+    (; alpha) = viscosity
+
+    return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4)
+end
+
 @doc raw"""
     ViscosityAdami(; nu)
 
@@ -180,4 +193,8 @@ end
     return visc .* v_diff
 end
 
+function kinematic_viscosity(system, viscosity::ViscosityAdami)
+    return viscosity.nu
+end
+
 @inline viscous_velocity(v, system, particle) = current_velocity(v, system, particle)

test/callbacks/callbacks.jl

@@ -1,3 +1,4 @@
 @testset verbose=true "Callbacks" begin
     include("info.jl")
+    include("stepsize.jl")
 end

test/callbacks/stepsize.jl

@@ -0,0 +1,17 @@
+@testset verbose=true "StepsizeCallback" begin
+    @testset verbose=true "show" begin
+        callback = StepsizeCallback(cfl=1.2)
+
+        show_compact = "StepsizeCallback(is_constant=true, cfl_number=1.2)"
+        @test repr(callback) == show_compact
+
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ StepsizeCallback                                                                                 │
+        │ ════════════════                                                                                 │
+        │ is constant: ……………………………………………… true                                                             │
+        │ CFL number: ………………………………………………… 1.2                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+        @test repr("text/plain", callback) == show_box
+    end
+end

test/examples/dam_break_2d_corrections.jl

@@ -60,14 +60,16 @@
     )
 
     @testset "continuity_reinit" begin
-        trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-                      fluid_particle_spacing=particle_spacing,
-                      smoothing_length=3.0 * particle_spacing,
-                      boundary_density_calculator=ContinuityDensity(),
-                      fluid_density_calculator=ContinuityDensity(),
-                      correction=nothing, use_reinit=true,
-                      prefix="continuity_reinit", tspan=tspan,
-                      fluid_density=fluid_density, density_diffusion=nothing)
+        @test_nowarn_mod trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+                                       fluid_particle_spacing=particle_spacing,
+                                       smoothing_length=3.0 * particle_spacing,
+                                       boundary_density_calculator=ContinuityDensity(),
+                                       fluid_density_calculator=ContinuityDensity(),
+                                       correction=nothing, use_reinit=true,
+                                       prefix="continuity_reinit", tspan=tspan,
+                                       fluid_density=fluid_density,
+                                       density_diffusion=nothing)
 
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
@@ -82,17 +84,27 @@
         println("="^100)
         println("fluid/dam_break_2d.jl with ", correction_name)
 
-        trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-                      fluid_particle_spacing=particle_spacing,
-                      smoothing_length=smoothing_length,
-                      boundary_density_calculator=SummationDensity(),
-                      fluid_density_calculator=fluid_density_calculator,
-                      correction=correction, use_reinit=false,
-                      clip_negative_pressure=(fluid_density_calculator isa SummationDensity),
-                      smoothing_kernel=smoothing_kernel,
-                      prefix="$(correction_name)", tspan=tspan,
-                      fluid_density=fluid_density, density_diffusion=nothing,
-                      boundary_layers=5)
+        if correction_name == "blended_gradient_continuity_correction"
+            # This simulation requires smaller time steps for some reason
+            cfl_ = 0.7
+        else
+            cfl_ = 1.1
+        end
+
+        @test_nowarn_mod trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+                                       fluid_particle_spacing=particle_spacing,
+                                       smoothing_length=smoothing_length,
+                                       boundary_density_calculator=SummationDensity(),
+                                       fluid_density_calculator=fluid_density_calculator,
+                                       correction=correction, use_reinit=false,
+                                       clip_negative_pressure=(fluid_density_calculator isa
+                                                               SummationDensity),
+                                       smoothing_kernel=smoothing_kernel,
+                                       prefix="$(correction_name)", tspan=tspan,
+                                       fluid_density=fluid_density,
+                                       density_diffusion=nothing,
+                                       boundary_layers=5, cfl=cfl_)
 
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0