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/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 +``` 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..b473aa6cf --- /dev/null +++ b/src/general/time_integration.jl @@ -0,0 +1,13 @@ +# 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). +This scheme is used by DualSPHysics: +https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#372-symplectic-position-verlet-scheme +""" +SymplecticPositionVerlet(_) = nothing