Implement the time integration method used in DualSPHysics and add docs for time integration #716
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
| @@ -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. | ||||||
|
Comment on lines
+13
to
+14
|
||||||
| 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}``, | ||||||
|
There was a problem hiding this comment. Function kick below is (u,v, t) here it is (u,t) this difference should be explained |
||||||
| 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. | ||||||
|
Comment on lines
+84
to
+85
There was a problem hiding this comment. Weird jump. Please include the explanation in the previous text. Also redundant because of the warning below. |
||||||
|
|
||||||
| !!! 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 | ||||||
|
There was a problem hiding this comment.
Suggested change
|
||||||
| 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. | ||||||
|
There was a problem hiding this comment.
Suggested change
|
||||||
|
|
||||||
| 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 | ||||||
| ``` | ||||||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change | ||||||
|---|---|---|---|---|---|---|---|---|
| @@ -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. | ||||||||
svchb marked this conversation as resolved.
Show resolved
Hide resolved
|
||||||||
| 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, | ||||||||
|
Comment on lines
+49
to
+50
There was a problem hiding this comment.
Suggested change
|
||||||||
| ::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) | ||||||||
|
There was a problem hiding this comment.
Suggested change
|
||||||||
| 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") | ||||||||
|
There was a problem hiding this comment.
Suggested change
|
||||||||
| 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) | ||||||||
|
Comment on lines
+83
to
+84
There was a problem hiding this comment.
Suggested change
|
||||||||
| 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 | ||||||||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also cite a DualSPHysics paper in which this is explained.