From b00d796dfed78292a65ea35e20edf3bcb8a2d333 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 3 Jan 2025 12:01:44 +0100 Subject: [PATCH 1/8] Implement the time integration method used in DualSPHysics --- Project.toml | 7 + ext/TrixiParticlesOrdinaryDiffEqExt.jl | 175 +++++++++++++++++++++++++ src/TrixiParticles.jl | 1 + src/general/general.jl | 1 + src/general/time_integration.jl | 11 ++ 5 files changed, 195 insertions(+) create mode 100644 ext/TrixiParticlesOrdinaryDiffEqExt.jl create mode 100644 src/general/time_integration.jl diff --git a/Project.toml b/Project.toml index 398268160..fa7f8d0ba 100644 --- a/Project.toml +++ b/Project.toml @@ -31,6 +31,12 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +[weakdeps] +OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" + +[extensions] +TrixiParticlesOrdinaryDiffEqExt = "OrdinaryDiffEqCore" + [compat] Adapt = "3, 4" CSV = "0.10" @@ -44,6 +50,7 @@ GPUArraysCore = "0.1, 0.2" JSON = "0.21" KernelAbstractions = "0.9" MuladdMacro = "0.2" +OrdinaryDiffEqCore = "1.14.1" PointNeighbors = "0.4.2" Polyester = "0.7.10" RecipesBase = "1" diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl new file mode 100644 index 000000000..51e1c3aae --- /dev/null +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -0,0 +1,175 @@ +module TrixiParticlesOrdinaryDiffEqExt + +using TrixiParticles + +# This is needed because `TrixiParticles.@threaded` translates +# to `PointNeighbors.parallel_foreach`, so `PointNeighbors` must be available. +const PointNeighbors = TrixiParticles.PointNeighbors + +using OrdinaryDiffEqCore: @.., @muladd, @cache, OrdinaryDiffEqCore, + OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache + +struct SymplecticPositionVerlet <: OrdinaryDiffEqPartitionedAlgorithm end + +TrixiParticles.SymplecticPositionVerlet() = SymplecticPositionVerlet() + +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 + +# Copied from OrdinaryDiffEqSymplecticRK +# +# provide the mutable uninitialized objects to keep state and derivative in case of mutable caches +# no such objects are required for constant caches +function alloc_symp_state(integrator) + (integrator.u.x..., integrator.cache.tmp.x...) +end + +# load state and derivatives at begin of symplectic iteration steps +function load_symp_state(integrator) + (integrator.uprev.x..., integrator.fsallast.x...) +end + +# store state and derivatives at the end of symplectic iteration steps +function store_symp_state!(integrator, cache, kdu, ku) + copyto!(integrator.k[1].x[1], integrator.k[2].x[1]) + copyto!(integrator.k[1].x[2], integrator.k[2].x[2]) + copyto!(integrator.k[2].x[2], ku) + copyto!(integrator.k[2].x[1], kdu) + nothing +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} + 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) + TrixiParticles.@threaded system for particle in TrixiParticles.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) + return du_system +end + +@muladd function update_density!(du_system, kdu_system, duprev_system, + ::ContinuityDensity, system, dt) + TrixiParticles.@threaded system for particle in TrixiParticles.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 diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d85096860..7f24cab2d 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -88,5 +88,6 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk export SurfaceTensionAkinci, CohesionForceAkinci +export SymplecticPositionVerlet end # module diff --git a/src/general/general.jl b/src/general/general.jl index dfd709861..9307671cd 100644 --- a/src/general/general.jl +++ b/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") diff --git a/src/general/time_integration.jl b/src/general/time_integration.jl new file mode 100644 index 000000000..f5c6f9bed --- /dev/null +++ b/src/general/time_integration.jl @@ -0,0 +1,11 @@ +# Time integration is handles by the package OrdinaryDiffEq.jl. +# See the docs for more details. +# In this file, we define the structs for extra time integration schemes that +# are implemented in the package extension TrixiParticlesOrdinaryDiffEqExt.jl. +""" + SymplecticPositionVerlet() + +Modified leapfrog integration scheme for Weakly Compressible SPH (WCSPH) when integrating +the density with [`ContinuityDensity`](@ref). +""" +SymplecticPositionVerlet(_) = nothing From 73382283ea58e48480076cca6d9ea2c6d525be5a Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 3 Jan 2025 12:02:52 +0100 Subject: [PATCH 2/8] Add source --- src/general/time_integration.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/general/time_integration.jl b/src/general/time_integration.jl index f5c6f9bed..b473aa6cf 100644 --- a/src/general/time_integration.jl +++ b/src/general/time_integration.jl @@ -7,5 +7,7 @@ Modified leapfrog integration scheme for Weakly Compressible SPH (WCSPH) when integrating the density with [`ContinuityDensity`](@ref). +This scheme is used by DualSPHysics: +https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme """ SymplecticPositionVerlet(_) = nothing From 479874e63782c2a9366de0685e15d141d2dd1598 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 3 Jan 2025 12:12:02 +0100 Subject: [PATCH 3/8] Add docstring to docs --- docs/src/time_integration.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index e69de29bb..9e8c5d59d 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -0,0 +1,3 @@ +```@docs + SymplecticPositionVerlet +``` From 12e7dbd7eb68f8e512caffe2cd6747edd6050077 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 11 Feb 2025 12:43:41 +0100 Subject: [PATCH 4/8] Revise package extension --- ext/TrixiParticlesOrdinaryDiffEqExt.jl | 54 ++++++++++---------------- src/general/time_integration.jl | 8 +++- 2 files changed, 27 insertions(+), 35 deletions(-) diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl index 51e1c3aae..9ed44d2c1 100644 --- a/ext/TrixiParticlesOrdinaryDiffEqExt.jl +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -1,18 +1,25 @@ module TrixiParticlesOrdinaryDiffEqExt -using TrixiParticles +using TrixiParticles: TrixiParticles, @threaded, each_moving_particle -# This is needed because `TrixiParticles.@threaded` translates +# This is needed because `@threaded` translates # to `PointNeighbors.parallel_foreach`, so `PointNeighbors` must be available. -const PointNeighbors = TrixiParticles.PointNeighbors +using TrixiParticles.PointNeighbors: PointNeighbors -using OrdinaryDiffEqCore: @.., @muladd, @cache, OrdinaryDiffEqCore, - OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache +using OrdinaryDiffEq.OrdinaryDiffEqCore: @.., @muladd, @cache, OrdinaryDiffEqCore, + OrdinaryDiffEqPartitionedAlgorithm, + OrdinaryDiffEqMutableCache +using OrdinaryDiffEq: alloc_symp_state, load_symp_state, store_symp_state! + +# 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} <: @@ -29,33 +36,10 @@ function OrdinaryDiffEqCore.get_fsalfirstlast(cache::SymplecticPositionVerletCac return (cache.fsalfirst, cache.k) end -# Copied from OrdinaryDiffEqSymplecticRK -# -# provide the mutable uninitialized objects to keep state and derivative in case of mutable caches -# no such objects are required for constant caches -function alloc_symp_state(integrator) - (integrator.u.x..., integrator.cache.tmp.x...) -end - -# load state and derivatives at begin of symplectic iteration steps -function load_symp_state(integrator) - (integrator.uprev.x..., integrator.fsallast.x...) -end - -# store state and derivatives at the end of symplectic iteration steps -function store_symp_state!(integrator, cache, kdu, ku) - copyto!(integrator.k[1].x[1], integrator.k[2].x[1]) - copyto!(integrator.k[1].x[2], integrator.k[2].x[2]) - copyto!(integrator.k[2].x[2], ku) - copyto!(integrator.k[2].x[1], kdu) - nothing -end - function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, - uprev, - uprev2, f, t, + uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, @@ -70,12 +54,13 @@ end function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, - uprev, - uprev2, f, t, + 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 @@ -112,6 +97,7 @@ end # 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 @@ -143,7 +129,7 @@ end @muladd function update_velocity!(du_system, kdu_system, duprev_system, system::WeaklyCompressibleSPHSystem, dt) - TrixiParticles.@threaded system for particle in TrixiParticles.each_moving_particle(system) + @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] @@ -159,12 +145,14 @@ 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) - TrixiParticles.@threaded system for particle in TrixiParticles.each_moving_particle(system) + @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 diff --git a/src/general/time_integration.jl b/src/general/time_integration.jl index b473aa6cf..1e620c3d3 100644 --- a/src/general/time_integration.jl +++ b/src/general/time_integration.jl @@ -1,4 +1,4 @@ -# Time integration is handles by the package OrdinaryDiffEq.jl. +# Time integration is handled by the package OrdinaryDiffEq.jl. # See the docs for more details. # In this file, we define the structs for extra time integration schemes that # are implemented in the package extension TrixiParticlesOrdinaryDiffEqExt.jl. @@ -9,5 +9,9 @@ Modified leapfrog integration scheme for Weakly Compressible SPH (WCSPH) when in the density with [`ContinuityDensity`](@ref). This scheme is used by DualSPHysics: https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme + +See [time integration](@ref time_integration) for more details. """ -SymplecticPositionVerlet(_) = nothing +function SymplecticPositionVerlet(_...) + error("the package OrdinaryDiffEq.jl needs to be loaded to use this scheme.") +end From 1b5b4e26a797e09a1ff7675e48f44a4575da2757 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 11 Feb 2025 17:39:09 +0100 Subject: [PATCH 5/8] Add docs for time integration --- Project.toml | 5 +- docs/src/refs.bib | 18 ++++ docs/src/time_integration.md | 141 +++++++++++++++++++++++++ ext/TrixiParticlesOrdinaryDiffEqExt.jl | 21 ++-- src/general/semidiscretization.jl | 5 +- src/general/time_integration.jl | 3 +- 6 files changed, 181 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index fa7f8d0ba..0224f4f34 100644 --- a/Project.toml +++ b/Project.toml @@ -32,10 +32,11 @@ 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 = "OrdinaryDiffEqCore" +TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"] [compat] Adapt = "3, 4" @@ -50,7 +51,7 @@ GPUArraysCore = "0.1, 0.2" JSON = "0.21" KernelAbstractions = "0.9" MuladdMacro = "0.2" -OrdinaryDiffEqCore = "1.14.1" +OrdinaryDiffEq = "6.91" PointNeighbors = "0.4.2" Polyester = "0.7.10" RecipesBase = "1" diff --git a/docs/src/refs.bib b/docs/src/refs.bib index e1b287b07..957f8ee88 100644 --- a/docs/src/refs.bib +++ b/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}, diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index 9e8c5d59d..96f7b2fbe 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -1,3 +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 ``` diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl index 9ed44d2c1..0046d5861 100644 --- a/ext/TrixiParticlesOrdinaryDiffEqExt.jl +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -1,16 +1,23 @@ module TrixiParticlesOrdinaryDiffEqExt -using TrixiParticles: TrixiParticles, @threaded, each_moving_particle +# 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. -# This is needed because `@threaded` translates +# We need to load the name `PointNeighbors` because `@threaded` translates # to `PointNeighbors.parallel_foreach`, so `PointNeighbors` must be available. -using TrixiParticles.PointNeighbors: PointNeighbors +using TrixiParticles: TrixiParticles, @threaded, each_moving_particle, + WeaklyCompressibleSPHSystem, ContinuityDensity, + PointNeighbors -using OrdinaryDiffEq.OrdinaryDiffEqCore: @.., @muladd, @cache, OrdinaryDiffEqCore, - OrdinaryDiffEqPartitionedAlgorithm, - OrdinaryDiffEqMutableCache +using OrdinaryDiffEq.OrdinaryDiffEqSymplecticRK: alloc_symp_state, load_symp_state, + store_symp_state! -using OrdinaryDiffEq: 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 diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..741e02e74 100644 --- a/src/general/semidiscretization.jl +++ b/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) diff --git a/src/general/time_integration.jl b/src/general/time_integration.jl index 1e620c3d3..f7da46d07 100644 --- a/src/general/time_integration.jl +++ b/src/general/time_integration.jl @@ -8,7 +8,8 @@ Modified leapfrog integration scheme for Weakly Compressible SPH (WCSPH) when integrating the density with [`ContinuityDensity`](@ref). This scheme is used by DualSPHysics: -https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme +[https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme] +(https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme). See [time integration](@ref time_integration) for more details. """ From 429b0a17fe23f14e3cbb45576ef48c7c8b3c9211 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 11 Feb 2025 18:19:34 +0100 Subject: [PATCH 6/8] Add tests --- test/examples/examples.jl | 56 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index b6d3b9f2e..0fc624f64 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -366,6 +366,62 @@ end include("dam_break_2d_corrections.jl") + + @testset "`SymplecticPositionVerlet` 2D" begin + @testset "2D unstable" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d.jl"), + tspan=(0, 0.1), sol=nothing) [ + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n"] + + sol = solve(ode, SymplecticPositionVerlet(), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + + # Unstable with this CFL + @test maximum(abs, sol.u[end]) > 2^15 + end + + @testset "2D stable" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d.jl"), + tspan=(0, 0.1), sol=nothing, + cfl=0.25) [ + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n"] + + sol = solve(ode, SymplecticPositionVerlet(), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + + # Stable with this CFL + @test maximum(abs, sol.u[end]) < 2^15 + end + + @testset "3D" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_3d.jl"), + fluid_particle_spacing=0.1, + tspan=(0, 0.1), sol=nothing) + stepsize_callback = StepsizeCallback(cfl=0.65) + callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback) + + sol = solve(ode, SymplecticPositionVerlet(), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + @test maximum(abs, sol.u[end]) < 2^15 + end + end end @testset verbose=true "Solid" begin From 4b6f1467c49770d49d2bd92970fe6eba654c806c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 18 Feb 2025 10:40:35 +0100 Subject: [PATCH 7/8] Implement suggestions --- docs/src/refs.bib | 12 ++++++++++++ docs/src/time_integration.md | 15 +++++++-------- ext/TrixiParticlesOrdinaryDiffEqExt.jl | 10 ++++------ src/general/time_integration.jl | 6 ++++-- 4 files changed, 27 insertions(+), 16 deletions(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 957f8ee88..1b0eb5ca0 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -264,6 +264,18 @@ @Article{DiRenzo2004 } +@article{Dominguez2022, + title={DualSPHysics: from fluid dynamics to multiphysics problems}, + author={Domínguez, Jose M and Fourtakas, Georgios and others}, + journal={Computational Particle Mechanics}, + volume={9}, + pages={867--895}, + year={2022}, + publisher={Springer}, + doi={10.1007/s40571-021-00404-2} +} + + @Article{Ferrari2009, author = {Ferrari, Angela and Dumbser, Michael and Toro, Eleuterio F. and Armanini, Aronne}, journal = {Computers \& Fluids}, diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index 96f7b2fbe..62c7b5ca3 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -10,8 +10,8 @@ In particular, a [`DynamicalODEProblem`](https://docs.sciml.ai/DiffEqDocs/stable 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. +This approach allows us to use [specialized time integration methods that do not work with +general `ODEProblem`s](@ref symplectic_schemes). 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)). @@ -56,7 +56,7 @@ We commonly use the following three schemes: - `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](@id symplectic_schemes) Symplectic schemes like the leapfrog method are often used for SPH. @@ -81,12 +81,10 @@ 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. + They will require very small time steps due to first-order convergence in the viscid case. ### Leapfrog drift-kick-drift @@ -115,12 +113,13 @@ 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 +to ``v`` as an 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. +See [Domínguez et al. 2022, Section 2.5.2](@cite Dominguez2022). 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. +even though it is actually included in the ``\operatorname{kick}`` as an additional dimension. This scheme reads ```math diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl index 0046d5861..181bf412c 100644 --- a/ext/TrixiParticlesOrdinaryDiffEqExt.jl +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -46,8 +46,7 @@ end function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, - uprev, uprev2, f, t, - dt, reltol, p, calck, + uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} @@ -55,7 +54,7 @@ function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_pro k = zero(rate_prototype) fsalfirst = zero(rate_prototype) half = uEltypeNoUnits(1 // 2) - SymplecticPositionVerletCache(u, uprev, k, tmp, fsalfirst, half) + SymplecticPositionVerletCache(u, uprev, tmp, k, fsalfirst, half) end function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_prototype, @@ -68,7 +67,7 @@ function OrdinaryDiffEqCore.alg_cache(alg::SymplecticPositionVerlet, u, rate_pro 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") + error("`SymplecticPositionVerlet` supports only in-place functions") end function OrdinaryDiffEqCore.initialize!(integrator, @@ -80,8 +79,7 @@ function OrdinaryDiffEqCore.initialize!(integrator, 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) + integrator.f.f2(integrator.k[2].x[2], duprev, uprev, integrator.p, integrator.t) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) integrator.stats.nf2 += 1 end diff --git a/src/general/time_integration.jl b/src/general/time_integration.jl index f7da46d07..2d1827b1a 100644 --- a/src/general/time_integration.jl +++ b/src/general/time_integration.jl @@ -7,9 +7,11 @@ Modified leapfrog integration scheme for Weakly Compressible SPH (WCSPH) when integrating the density with [`ContinuityDensity`](@ref). -This scheme is used by DualSPHysics: +This scheme is used by the SPH code [DualSPHysics](https://github.com/DualSPHysics/DualSPHysics). +See [https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme] -(https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme). +(https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme) +and [Domínguez et al. 2022, Section 2.5.2](@cite Dominguez2022). See [time integration](@ref time_integration) for more details. """ From 5927b8899f60c2c3a9b12845efc4a2c3ca73f5cd Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 20 Feb 2025 10:37:48 +0100 Subject: [PATCH 8/8] Update NEWS.md --- NEWS.md | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index e1504a4ea..ba57c5694 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,16 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.2.7 + +### Features + +- The symplectic time integration scheme from DualSPHysics was added (#716) + +### Documentation + +- Documentation for time integration was added (#716) + ## Version 0.2.6 ### Features @@ -11,7 +21,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Refactored -- Surface normal calculation was moved from surface_tension.jl to surface_normal_sph.jl (#539) +- Surface normal calculation was moved from `surface_tension.jl` to `surface_normal_sph.jl` (#539) ## Version 0.2.5