Implement the time integration method used in DualSPHysics and add docs for time integration

Project.toml

@@ -31,6 +31,13 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
+[weakdeps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+
+[extensions]
+TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"]
+
 [compat]
 Adapt = "4"
 CSV = "0.10"
@@ -44,6 +51,7 @@ GPUArraysCore = "0.2"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
+OrdinaryDiffEq = "6.91"
 PointNeighbors = "0.4.7"
 Polyester = "0.7.10"
 RecipesBase = "1"

docs/src/refs.bib

@@ -171,6 +171,13 @@ @Article{Bonet1999
   publisher = {Elsevier BV},
 }
 
+@misc{CarpenterKennedy,
+  author    = {Carpenter, Mark H. and Kennedy, Christopher A.},
+  title     = {Fourth Order 2N-storage Runge-Kutta Schemes},
+  year      = {1994},
+  url       = {https://ntrs.nasa.gov/citations/19940028444},
+}
+
 
 @Article{Clausen2013,
   author    = {Clausen, Jonathan R.},
@@ -566,6 +573,17 @@ @Article{Ramachandran2019
   publisher = {Elsevier BV},
 }
 
+@Article{Ranocha2022,
+  author    = {Ranocha, Hendrik and Dalcin, Lisandro and Parsani, Matteo and Ketcheson, David I.},
+  journal   = {Communications on Applied Mathematics and Computation},
+  title     = {Optimized Runge-Kutta Methods with Automatic Step Size Control for Compressible Computational Fluid Dynamics},
+  year      = {2022},
+  month     = dec,
+  pages     = {1191--1228},
+  volume    = {4},
+  doi       = {10.1007/s42967-021-00159-w},
+}
+
 @Article{Schoenberg1946,
   author    = {Schoenberg, I. J.},
   journal   = {Quarterly of Applied Mathematics},

docs/src/time_integration.md

@@ -0,0 +1,144 @@
+# [Time integration](@id time_integration)
+
+TrixiParticles.jl uses a modular approach where time integration is just another module
+that can be customized and exchanged.
+The function [`semidiscretize`](@ref) returns an `ODEProblem`
+(see [the OrdinaryDiffEq.jl docs](https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/)),
+which can be integrated with [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
+
+In particular, a [`DynamicalODEProblem`](https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/)
+is returned, where the right-hand side is split into two functions, the `kick!`, which
+computes the derivative of the particle velocities and the `drift!`, which computes
+the derivative of the particle positions.
+This approach allows us to use specialized time integration methods that do not work with
+general `ODEProblem`s.
+Note that this is not a true `DynamicalODEProblem` where the kick does not depend
+on the velocity. Therefore, not all integrators designed for `DynamicalODEProblem`s
+will work (properly) (see [below](@ref kick_drift_kick)).
+However, all integrators designed for general `ODEProblem`s can be used.
+
+## Usage
+
+After obtaining an `ODEProblem` from [`semidiscretize`](@ref), let us call it `ode`,
+we can pass it to the function `solve` of OrdinaryDiffEq.jl.
+For most schemes, we do the following:
+```julia
+using OrdinaryDiffEq
+sol = solve(ode, Euler(),
+            dt=1.0,
+            save_everystep=false, callback=callbacks);
+```
+Here, `Euler()` should in practice be replaced by a more useful scheme.
+`callbacks` should be a `CallbackSet` containing callbacks like the [`InfoCallback`](@ref).
+For callbacks, please refer to [the docs](@ref Callbacks) and the example files.
+In this case, we need to either set a reasonable, problem- and resolution-dependent
+step size `dt`, or use the [`StepsizeCallback`](@ref), which overwrites the step size
+dynamically during the simulation based on a CFL-number.
+
+Some schemes, e.g. the two schemes `RDPK3SpFSAL35` and `RDPK3SpFSAL49` mentioned below,
+support automatic time stepping, where the step size is determined automatically based on
+error estimates during the simulation.
+These schemes do not use the keyword argument `dt` and will ignore the step size set by
+the [`StepsizeCallback`](@ref).
+Instead, they will work with the tolerances `abstol` and `reltol`, which can be passed as
+keyword arguments to `solve`. The default tolerances are `abstol=1e-6` and `reltol=1e-3`.
+
+## Recommended schemes
+
+A list of schemes for general `ODEProblem`s can be found
+[here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).
+We commonly use the following three schemes:
+- `CarpenterKennedy2N54(williamson_condition=false)`: A five-stage, fourth order
+  low-storage Runge-Kutta method designed by [Carpenter and Kennedy](@cite CarpenterKennedy)
+  for hyperbolic problems.
+- `RDPK3SpFSAL35()`: A 5-stage, third order low-storage Runge-Kutta scheme with embedded
+  error estimator, optimized for compressible fluid mechanics [Ranocha2022](@cite).
+- `RDPK3SpFSAL49()`: A 9-stage, fourth order low-storage Runge-Kutta scheme with embedded
+  error estimator, optimized for compressible fluid mechanics [Ranocha2022](@cite).
+
+## Symplectic schemes
+
+Symplectic schemes like the leapfrog method are often used for SPH.
+
+### [Leapfrog kick-drift-kick](@id kick_drift_kick)
+
+The kick-drift-kick scheme of the leapfrog method, updating positions ``u``
+and velocities ``v`` with the functions ``\operatorname{kick}`` and ``\operatorname{drift}``,
+reads:
+```math
+\begin{align*}
+v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(u^0, t^0), \\
+u^1 &= u^0 + \Delta t\, \operatorname{drift} \left( v^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\
+v^1 &= v^{1/2} + \frac{1}{2} \Delta t\, \operatorname{kick}(u^{1}, t^0 + \Delta t).
+\end{align*}
+```
+In this form, it is also identical to the velocity Verlet scheme.
+Note that this only works as long as ``\operatorname{kick}`` does not depend on ``v``, i.e.,
+in the inviscid case.
+Once we add viscosity, ``\operatorname{kick}`` depends on both ``u`` and ``v``.
+Then, the calculation of ``v^1`` requires ``v^1`` and becomes implicit.
+
+The way this scheme is implemented in OrdinaryDiffEq.jl as `VerletLeapfrog`,
+the intermediate velocity ``v^{1/2}`` is passed to ``\operatorname{kick}`` in the last stage,
+resulting in first-order convergence when the scheme is used in the viscid case.
+Note that the method `VelocityVerlet` does not work with TrixiParticles.jl for the same and
+other additional reasons.
+
+!!! warning
+    Please do not use `VelocityVerlet` and `VerletLeapfrog` with TrixiParticles.jl.
+    They will require very small time steps due to first-order convergence.
+
+### Leapfrog drift-kick-drift
+
+The drift-kick-drift scheme of the leapfrog method reads:
+```math
+\begin{align*}
+u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, \operatorname{drift}(v^0, t^0), \\
+v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\
+u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, t^0 + \Delta t).
+\end{align*}
+```
+In the viscid case where ``\operatorname{kick}`` depends on ``v``, we can add another
+half step for ``v``, yielding
+```math
+\begin{align*}
+u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, \operatorname{drift}(v^0, u^0, t^0), \\
+v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\
+v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\
+u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t).
+\end{align*}
+```
+This scheme is implemented in OrdinaryDiffEq.jl as `LeapfrogDriftKickDrift` and yields
+the desired second order as long as ``\operatorname{drift}`` does not depend on ``u``,
+which is always the case.
+
+### Symplectic position Verlet
+
+When the density is integrated (with [`ContinuityDensity`](@ref)), the density is appended
+to ``v`` as additional dimension, so all previously mentioned schemes treat the density
+similar to the velocity.
+The SPH code [DualSPHysics](https://github.com/DualSPHysics/DualSPHysics) implements
+a variation of the drift-kick-drift scheme where the density is updated separately.
+In the following, we will call the derivative of the density ``R(v, u, t)``,
+even though it is actually included in the ``\operatorname{kick}`` as additional dimension.
+
+This scheme reads
+```math
+\begin{align*}
+u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, \operatorname{drift}(v^0, u^0, t^0), \\
+v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\
+\rho^{1/2} &= \rho^0 + \frac{1}{2} \Delta t\, R(v^0, u^0, t^0), \\
+v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\
+\rho^1 &= \rho^0 \frac{2 - \varepsilon^{1/2}}{2 + \varepsilon^{1/2}}, \\
+u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t),
+\end{align*}
+```
+where
+```math
+\varepsilon^{1/2} = - \frac{R(v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t)}{\rho^{1/2}} \Delta t.
+```
+This scheme is implemented in TrixiParticles.jl as [`SymplecticPositionVerlet`](@ref).
+
+```@docs
+    SymplecticPositionVerlet
+```

ext/TrixiParticlesOrdinaryDiffEqExt.jl

@@ -0,0 +1,170 @@
+module TrixiParticlesOrdinaryDiffEqExt
+
+# This package extension defines the `SymplecticPositionVerlet` scheme from DualSPHysics.
+# The scheme is similar to the `LeapfrogDriftKickDrift` scheme, but with a different
+# update for the density.
+# See https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme
+# and the TrixiParticles.jl docs on time integration for more details.
+
+# We need to load the name `PointNeighbors` because `@threaded` translates
+# to `PointNeighbors.parallel_foreach`, so `PointNeighbors` must be available.
+using TrixiParticles: TrixiParticles, @threaded, each_moving_particle,
+                      WeaklyCompressibleSPHSystem, ContinuityDensity,
+                      PointNeighbors
+
+using OrdinaryDiffEq.OrdinaryDiffEqSymplecticRK: alloc_symp_state, load_symp_state,
+                                                 store_symp_state!
+
+using OrdinaryDiffEqCore: OrdinaryDiffEqCore, @.., @muladd, @cache,
+                          OrdinaryDiffEqPartitionedAlgorithm,
+                          OrdinaryDiffEqMutableCache
+
+# Define a new struct for the SymplecticPositionVerlet scheme
+struct SymplecticPositionVerlet <: OrdinaryDiffEqPartitionedAlgorithm end
+
+# Overwrite the function in TrixiParticles to use the new scheme
+TrixiParticles.SymplecticPositionVerlet() = SymplecticPositionVerlet()
+
+# The following is similar to the definition of the `VerletLeapfrog` scheme
+# and the corresponding cache in OrdinaryDiffEq.jl.
+OrdinaryDiffEqCore.default_linear_interpolation(alg::SymplecticPositionVerlet, prob) = true
+
+@cache struct SymplecticPositionVerletCache{uType, rateType, uEltypeNoUnits} <:
+              OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    tmp::uType
+    k::rateType
+    fsalfirst::rateType
+    half::uEltypeNoUnits
+end
+
+function OrdinaryDiffEqCore.get_fsalfirstlast(cache::SymplecticPositionVerletCache, u)
+    return (cache.fsalfirst, cache.k)
+end
+
+function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype,
+                                      ::Type{uEltypeNoUnits},
+                                      ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
+                                      uprev, uprev2, f, t,
+                                      dt, reltol, p, calck,
+                                      ::Val{true}) where {uEltypeNoUnits,
+                                                          uBottomEltypeNoUnits,
+                                                          tTypeNoUnits}
+    tmp = zero(u)
+    k = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    half = uEltypeNoUnits(1 // 2)
+    SymplecticPositionVerletCache(u, uprev, k, tmp, fsalfirst, half)
+end
+
+function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype,
+                                      ::Type{uEltypeNoUnits},
+                                      ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
+                                      uprev, uprev2, f, t,
+                                      dt, reltol, p, calck,
+                                      ::Val{false}) where {uEltypeNoUnits,
+                                                           uBottomEltypeNoUnits,
+                                                           tTypeNoUnits}
+    # We only use inplace functions in TrixiParticles, so there is no point
+    # in implementing the non-inplace version.
+    error("`SymplecticPositionVerlet` only supports inplace functions")
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator,
+                                        cache::C) where {C <: SymplecticPositionVerletCache}
+    integrator.kshortsize = 2
+    resize!(integrator.k, integrator.kshortsize)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+
+    duprev, uprev = integrator.uprev.x
+    integrator.f.f1(integrator.k[2].x[1], duprev, uprev, integrator.p, integrator.t)
+    # verify_f2(integrator.f.f2, integrator.k[2].x[2], duprev, uprev, integrator.p,
+    #     integrator.t, integrator, cache)
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+    integrator.stats.nf2 += 1
+end
+
+@muladd function OrdinaryDiffEqCore.perform_step!(integrator,
+                                                  cache::SymplecticPositionVerletCache,
+                                                  repeat_step=false)
+    (; t, dt, f, p) = integrator
+    duprev, uprev, _, _ = load_symp_state(integrator)
+    du, u, kdu, ku = alloc_symp_state(integrator)
+
+    # update position half step
+    half = cache.half
+    f.f2(ku, duprev, uprev, p, t)
+    @.. broadcast=false u=uprev + dt * half * ku
+
+    # update velocity half step
+    f.f1(kdu, duprev, uprev, p, t)
+    @.. broadcast=false du=duprev + dt * half * kdu
+
+    # update velocity (add to previous full step velocity)
+    f.f1(kdu, du, u, p, t + half * dt)
+
+    # The following is equivalent to `du = duprev + dt * kdu` for the velocity, but when
+    # the density is integrated, a different update is used for the density.
+    semi = p
+    TrixiParticles.foreach_system(semi) do system
+        kdu_system = TrixiParticles.wrap_v(kdu, system, semi)
+        du_system = TrixiParticles.wrap_v(du, system, semi)
+        duprev_system = TrixiParticles.wrap_v(duprev, system, semi)
+
+        update_velocity!(du_system, kdu_system, duprev_system, system, dt)
+        update_density!(du_system, kdu_system, duprev_system, system, dt)
+    end
+
+    # update position (add to half step position)
+    f.f2(ku, du, u, p, t + dt)
+    @.. broadcast=false u=u + dt * half * ku
+
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2)
+    integrator.stats.nf2 += 2
+    store_symp_state!(integrator, cache, kdu, ku)
+end
+
+@muladd function update_velocity!(du_system, kdu_system, duprev_system, system, dt)
+    @.. broadcast=false du_system=duprev_system + dt * kdu_system
+end
+
+@inline function update_density!(du_system, kdu_system, duprev_system, system, dt)
+    return du_system
+end
+
+@muladd function update_velocity!(du_system, kdu_system, duprev_system,
+                                  system::WeaklyCompressibleSPHSystem, dt)
+    @threaded system for particle in each_moving_particle(system)
+        for i in 1:ndims(system)
+            du_system[i, particle] = duprev_system[i, particle] +
+                                     dt * kdu_system[i, particle]
+        end
+    end
+end
+
+@inline function update_density!(du_system, kdu_system, duprev_system,
+                                 system::WeaklyCompressibleSPHSystem, dt)
+    update_density!(du_system, kdu_system, duprev_system,
+                    system.density_calculator, system, dt)
+end
+
+@inline function update_density!(du_system, kdu_system, duprev_system,
+                                 density_calculator, system, dt)
+    # Don't do anything when the density is not integrated.
+    # This scheme is then equivalent to the `LeapfrogDriftKickDrift` scheme.
+    return du_system
+end
+
+@muladd function update_density!(du_system, kdu_system, duprev_system,
+                                 ::ContinuityDensity, system, dt)
+    @threaded system for particle in each_moving_particle(system)
+        density_prev = duprev_system[end, particle]
+        density_half = du_system[end, particle]
+        epsilon = -kdu_system[end, particle] / density_half * dt
+        du_system[end, particle] = density_prev * (2 - epsilon) / (2 + epsilon)
+    end
+end
+
+end # module

src/TrixiParticles.jl

@@ -89,7 +89,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure,
        max_density, min_density, avg_density
 export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d,
        interpolate_plane_2d_vtk
-export SurfaceTensionAkinci, CohesionForceAkinci
-export ColorfieldSurfaceNormal
+export SurfaceTensionAkinci, CohesionForceAkinci, ColorfieldSurfaceNormal
+export SymplecticPositionVerlet
 
 end # module

src/general/general.jl

@@ -9,3 +9,4 @@ include("buffer.jl")
 include("interpolation.jl")
 include("custom_quantities.jl")
 include("neighborhood_search.jl")
+include("time_integration.jl")

src/general/semidiscretization.jl

@@ -242,8 +242,9 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 A `DynamicalODEProblem` (see [the OrdinaryDiffEq.jl docs](https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/))
 to be integrated with [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
 Note that this is not a true `DynamicalODEProblem` where the acceleration does not depend on the velocity.
-Therefore, not all integrators designed for `DynamicalODEProblems` will work properly.
-However, all integrators designed for `ODEProblems` can be used.
+Therefore, not all integrators designed for `DynamicalODEProblem`s will work properly.
+However, all integrators designed for `ODEProblem`s can be used.
+See [time integration](@ref time_integration) for more details.
 
 # Examples
 ```jldoctest; output = false, filter = r"u0: .*", setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); ref_system = fluid_system)