diff --git a/NEWS.md b/NEWS.md index f1940cbc46..d4e2ccdd7d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,14 +4,27 @@ 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.4 + +### API Changes + +- Combined transport velocity formulation (TVF) and particle shifting technique (PST) into + one unified framework. + The keyword argument `transport_velocity` now changed to `shifting_technique`. + The `ParticleShiftingCallback` has been removed. To use PST, use the `UpdateCallback` + instead, and pass `shifting_technique=ParticleShiftingTechnique()` to the system. + +- Renamed the keyword argument `tlsph` to `place_on_shell` for `ParticlePackingSystem`, + `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape`. + ## Version 0.3.1 ### Features -- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, +- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, which implement a simplified Smagorinsky-type sub-grid-scale viscosity. (#753) -- **Multithreaded Integration Array**: Introduced a new array type for CPU backends +- **Multithreaded Integration Array**: Introduced a new array type for CPU backends that enables multithreaded broadcasting, delivering speed-ups of up to 5× on systems with many threads when combined with thread pinning. (#722) @@ -21,17 +34,17 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - **DXF file format support**: Import complex geometries using the DXF file format. (#821) - **Improved Plane interpolation**: Massively improved interpolation performance for planes (#763). - + ### GPU - Make PST GPU-compatible (#813). - + - Make open boundaries GPU-compatible (#773). - + - Make interpolation GPU-compatible (#812). ### Important Bugfixes - + - Fix validation setups (#801). - Calculate interpolated density instead of computed density when using interpolation (#808). diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 0db8d275e2..d035fb6fc1 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -45,64 +45,3 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")] ``` - -## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) -Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow. -To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation. -The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF -also for the [EDAC](@ref edac) scheme. - -The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by - -```math -\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a -``` - -and is obtained at every time-step ``\Delta t`` from - -```math -\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), -``` - -where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. -The tilde in the second term of the right hand side indicates that the material derivative has an advection part. - -The discretized form of the last term is - -```math - -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, -``` - -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. -Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, -which means that there is a non-vanishing contribution only when particles are disordered. -That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted. - -The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is - -```math -\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, -``` - - where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified - advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. - -The discretized form of the momentum equation for a particle ``a`` reads as - -```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. -``` - -Here, ``\tilde{p}_{ab}`` is the density-weighted pressure - -```math -\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, -``` - -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")] -``` diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 51355961f6..62a1bc8d0c 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -152,8 +152,74 @@ as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9. The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as ``\delta^+``-SPH. -The Particle Shifting Technique can be applied in form -of the [`ParticleShiftingCallback`](@ref). +To apply particle shifting, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + + +## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) + +An alternative formulation is the so-called Transport Velocity Formulation (TVF) +by [Adami (2013)](@cite Adami2013). +[Ramachandran (2019)](@cite Ramachandran2019) applied the TVF also for the [EDAC](@ref edac) +scheme. + +The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position +of the particle ``r_a`` from one time step to the next by +```math +\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a +``` +and is obtained at every time step ``\Delta t`` from +```math +\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), +``` +where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` +is a constant background pressure field. +The tilde in the second term of the right-hand side indicates that the material derivative +has an advection part. + +The discretized form of the last term is +```math + -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, +``` +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, +the discretization is not 0th-order consistent for **non**-uniform particle distribution, +which means that there is a non-vanishing contribution only when particles are disordered. +That also means that ``p_{\text{background}}`` occurs as pre-factor to correct +the trajectory of a particle resulting in uniform pressure distributions. +Suggested is a background pressure which is in the order of the reference pressure, +but it can be chosen arbitrarily large when the time-step criterion is adjusted. + +The inviscid momentum equation with an additional convection term for a particle +moving with ``\tilde{v}`` is +```math +\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, +``` +where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence +of the modified advection velocity and can be interpreted as the convection of momentum +with the relative velocity ``\tilde{v}-v``. + +The discretized form of the momentum equation for a particle ``a`` reads as +```math +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. +``` +Here, ``\tilde{p}_{ab}`` is the density-weighted pressure +```math +\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, +``` +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` +and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors +for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``, +as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. + +To apply the TVF, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "shifting_techniques.jl")] +``` + ## [Tensile Instability Control](@id tic) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index b6582a1701..1b8b15cda4 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -65,7 +65,7 @@ if wcsph state_equation, smoothing_kernel, pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) else state_equation = nothing density_calculator = ContinuityDensity() @@ -73,7 +73,7 @@ else smoothing_length, density_calculator=density_calculator, sound_speed, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) end # ========================================================================================== diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 850c9f2ff9..648a98be11 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -60,7 +60,7 @@ smoothing_length = 1.2 * particle_spacing smoothing_kernel = SchoenbergQuarticSplineKernel{2}() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(pressure), + shifting_technique=TransportVelocityAdami(pressure), acceleration=(acceleration_x, 0.0)) # ========================================================================================== diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 74770c318b..9cd4bd95a5 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -28,7 +28,7 @@ initial_fluid_size = tank_size initial_velocity = (1.0, 0.0) fluid_density = 1000.0 -sound_speed = initial_velocity[1] +sound_speed = 10 * initial_velocity[1] state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7) @@ -48,6 +48,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + shifting_technique=nothing, pressure_acceleration=nothing) # ========================================================================================== diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 3544ac7a8d..85e3461ccc 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -74,6 +74,7 @@ viscosity = ViscosityAdami(nu=kinematic_viscosity) fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=viscosity, density_calculator=fluid_density_calculator, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) # Alternatively the WCSPH scheme can be used @@ -86,6 +87,7 @@ if wcsph fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) end @@ -159,12 +161,10 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -particle_shifting = ParticleShiftingCallback() extra_callback = nothing -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), - particle_shifting, extra_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) sol = solve(ode, RDPK3SpFSAL35(), abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index e6e74fbc4b..8a99c455c6 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -86,13 +86,13 @@ if wcsph pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(background_pressure)) + shifting_technique=TransportVelocityAdami(background_pressure)) else density_calculator = SummationDensity() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, density_calculator=density_calculator, - transport_velocity=TransportVelocityAdami(background_pressure), + shifting_technique=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 987e208750..2a53025d46 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,10 +65,9 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - ParticleShiftingCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback export ContinuityDensity, SummationDensity -export PenaltyForceGanzenmueller, TransportVelocityAdami +export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 062a16e3a9..e9eb048c3f 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -32,4 +32,3 @@ include("post_process.jl") include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") -include("particle_shifting.jl") diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl deleted file mode 100644 index 921370f3b1..0000000000 --- a/src/callbacks/particle_shifting.jl +++ /dev/null @@ -1,148 +0,0 @@ -@doc raw""" - ParticleShiftingCallback() - -Callback to apply the Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). -Following the original paper, the callback is applied in every time step and not -in every stage of a multi-stage time integration method to reduce the computational -cost and improve the stability of the scheme. - -See [Callbacks](@ref Callbacks) for more information on how to use this callback. -See [Particle Shifting Technique](@ref shifting) for more information on the method itself. - -!!! warning - The Particle Shifting Technique needs to be disabled close to the free surface - and therefore requires a free surface detection method. This is not yet implemented. - **This callback cannot be used in a free surface simulation.** -""" -function ParticleShiftingCallback() - # The first one is the `condition`, the second the `affect!` - return DiscreteCallback((particle_shifting_condition), particle_shifting!, - save_positions=(false, false)) -end - -# `condition` -function particle_shifting_condition(u, t, integrator) - return true -end - -# `affect!` -function particle_shifting!(integrator) - t = integrator.t - semi = integrator.p - v_ode, u_ode = integrator.u.x - dt = integrator.dt - # Internal cache vector, which is safe to use as temporary array - vu_cache = first(get_tmp_cache(integrator)) - - @trixi_timeit timer() "particle shifting callback" begin - # Update quantities that are stored in the systems. These quantities (e.g. pressure) - # still have the values from the last stage of the previous step if not updated here. - @trixi_timeit timer() "update systems and nhs" begin - # Don't create sub-timers here to avoid cluttering the timer output - @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) - end - - @trixi_timeit timer() "particle shifting" foreach_system(semi) do system - u = wrap_u(u_ode, system, semi) - v = wrap_v(v_ode, system, semi) - particle_shifting!(u, v, system, v_ode, u_ode, semi, vu_cache, dt) - end - end - - # Tell OrdinaryDiffEq that `u` has been modified - u_modified!(integrator, true) - - return integrator -end - -function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) - return u -end - -function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, - vu_cache, dt) - # Wrap the cache vector to an NDIMS x NPARTICLES matrix. - # We need this buffer because we cannot safely update `u` while iterating over it. - _, u_cache = vu_cache.x - delta_r = wrap_u(u_cache, system, semi) - set_zero!(delta_r) - - # This has similar performance to `maximum(..., eachparticle(system))`, - # but is GPU-compatible. - v_max = maximum(x -> sqrt(dot(x, x)), - reinterpret(reshape, SVector{ndims(system), eltype(v)}, - current_velocity(v, system))) - - # TODO this needs to be adapted to multi-resolution. - # Section 3.2 explains what else needs to be changed. - dx = particle_spacing(system, 1) - Wdx = smoothing_kernel(system, dx, 1) - h = smoothing_length(system, 1) - - foreach_system(semi) do neighbor_system - u_neighbor = wrap_u(u_ode, neighbor_system, semi) - v_neighbor = wrap_v(v_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - rho_a = current_density(v, system, particle) - rho_b = current_density(v_neighbor, neighbor_system, neighbor) - - kernel = smoothing_kernel(system, distance, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # According to p. 29 below Eq. 9 - R = 2 // 10 - n = 4 - - # Eq. 7 in Sun et al. (2017). - # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h - # (see p. 29, right above Eq. 9), but this does not work when scaling h. - # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected - # for both small and large smoothing length factors. - # We need to scale - # - quadratically with the smoothing length, - # - linearly with the particle spacing, - # - linearly with the time step. - # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. - delta_r_ = -dt * v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * - m_b / (rho_a + rho_b) * grad_kernel - - # Write into the buffer - for i in eachindex(delta_r_) - @inbounds delta_r[i, particle] += delta_r_[i] - end - end - end - - # Add δ_r from the buffer to the current coordinates - @threaded semi for particle in eachparticle(system) - for i in axes(delta_r, 1) - @inbounds u[i, particle] += delta_r[i, particle] - end - end - - return u -end - -function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - print(io, "ParticleShiftingCallback()") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - summary_box(io, "ParticleShiftingCallback") - end -end diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 77496a46c5..8eecb8380a 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -74,22 +74,32 @@ function (update_callback!::UpdateCallback)(integrator) semi = integrator.p v_ode, u_ode = integrator.u.x - # Update quantities that are stored in the systems. These quantities (e.g. pressure) - # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t) - - # Update open boundaries first, since particles might be activated or deactivated - @trixi_timeit timer() "update open boundary" foreach_system(semi) do system - update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) - end - - @trixi_timeit timer() "update particle packing" foreach_system(semi) do system - update_particle_packing(system, v_ode, u_ode, semi, integrator) - end - - # This is only used by the particle packing system and should be removed in the future - @trixi_timeit timer() "update TVF" foreach_system(semi) do system - update_transport_velocity!(system, v_ode, semi) + @trixi_timeit timer() "update callback" begin + # Update quantities that are stored in the systems. These quantities (e.g. pressure) + # still have the values from the last stage of the previous step if not updated here. + @trixi_timeit timer() "update systems and nhs" begin + # Don't create sub-timers here to avoid cluttering the timer output + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) + end + + # Update open boundaries first, since particles might be activated or deactivated + @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + end + + @trixi_timeit timer() "update particle packing" foreach_system(semi) do system + update_particle_packing(system, v_ode, u_ode, semi, integrator) + end + + # This is only used by the particle packing system and should be removed in the future + @trixi_timeit timer() "update TVF" foreach_system(semi) do system + update_transport_velocity!(system, v_ode, semi) + end + + @trixi_timeit timer() "particle shifting" foreach_system(semi) do system + particle_shifting_from_callback!(u_ode, shifting_technique(system), system, + v_ode, semi, integrator.dt) + end end # Tell OrdinaryDiffEq that `u` has been modified diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9c78b3349..9532d812ab 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -485,7 +485,7 @@ end @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du @inline function add_velocity!(du, v, particle, system::FluidSystem) - # This is zero unless a transport velocity is used + # This is zero unless a shifting technique is used delta_v_ = delta_v(system, particle) for i in 1:ndims(system) diff --git a/src/general/system.jl b/src/general/system.jl index 929fb8b425..5b4d5aecb8 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -147,9 +147,6 @@ function update_final!(system, v, u, v_ode, u_ode, semi, t) return system end -# Only for systems requiring the use of the `UpdateCallback` -@inline requires_update_callback(system) = false - @inline initial_smoothing_length(system) = smoothing_length(system, nothing) @inline function smoothing_length(system, particle) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index b93e25f324..1aca732098 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -349,9 +349,6 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) else vtk["solver"] = "EDAC" vtk["sound_speed"] = system.sound_speed - vtk["background_pressure_TVF"] = system.transport_velocity isa Nothing ? - "-" : - system.transport_velocity.background_pressure end end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 08139e2ca6..a4d7d280a4 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -51,13 +51,12 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Add convection term (only when using `TransportVelocityAdami`) - dv_tvf = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = dv_shifting(shifting_technique(particle_system), + particle_system, neighbor_system, particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 19e50c441a..89d8a923ff 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -3,7 +3,7 @@ smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, + shifting_technique=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, surface_normal_method=nothing, buffer_size=nothing, @@ -30,10 +30,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles - from the local pressure (default: `true` when using TVF, `false` otherwise). + from the local pressure (default: `true` when using shifting, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -67,7 +68,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, acceleration :: SVector{NDIMS, ELTYPE} correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: TV average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT @@ -83,8 +84,8 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, - average_pressure_reduction=(!isnothing(transport_velocity)), + shifting_technique=nothing, + average_pressure_reduction=(!isnothing(shifting_technique)), alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -137,7 +138,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., create_cache_avg_pressure_reduction(initial_condition, avg_pressure_reduction)..., create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, @@ -161,14 +162,14 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), typeof(transport_velocity), + typeof(pressure_acceleration), typeof(shifting_technique), typeof(avg_pressure_reduction), typeof(source_terms), typeof(surface_tension), typeof(surface_normal_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, acceleration_, correction, - pressure_acceleration, transport_velocity, + pressure_acceleration, shifting_technique, avg_pressure_reduction, source_terms, surface_tension, surface_normal_method, buffer, @@ -218,8 +219,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) + summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "average pressure reduction", typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") summary_line(io, "acceleration", system.acceleration) @@ -271,7 +271,7 @@ end @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed -@inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity +@inline shifting_technique(system::EntropicallyDampedSPHSystem) = system.shifting_technique @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) average_pressure(system, system.average_pressure_reduction, particle) @@ -321,7 +321,7 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end # No average pressure reduction is used diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index d3e6adf971..59703f9b6d 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -193,7 +193,7 @@ end include("pressure_acceleration.jl") include("viscosity.jl") -include("transport_velocity.jl") +include("shifting_techniques.jl") include("surface_tension.jl") include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl new file mode 100644 index 0000000000..98775b3be0 --- /dev/null +++ b/src/schemes/fluid/shifting_techniques.jl @@ -0,0 +1,291 @@ +abstract type AbstractShiftingTechnique end + +# No shifting for a system by default +@inline shifting_technique(system) = nothing + +# WARNING: Be careful if defining this function for a specific system type. +# The version for a specific system type will override this generic version. +requires_update_callback(system) = requires_update_callback(shifting_technique(system)) +requires_update_callback(::Nothing) = false +requires_update_callback(::AbstractShiftingTechnique) = true + +# This is called from the `UpdateCallback` +particle_shifting_from_callback!(u_ode, shifting, system, v_ode, semi, dt) = u_ode + +create_cache_shifting(initial_condition, ::Nothing) = (;) + +function create_cache_shifting(initial_condition, ::AbstractShiftingTechnique) + delta_v = zeros(eltype(initial_condition), ndims(initial_condition), + nparticles(initial_condition)) + + return (; delta_v) +end + +# `δv` is the correction to the particle velocity due to the shifting. +# Particles are advected with the velocity `v + δv`. +@propagate_inbounds function delta_v(system, particle) + return delta_v(system, shifting_technique(system), particle) +end + +# Zero when no shifting is used +@inline function delta_v(system, shifting, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +@propagate_inbounds function delta_v(system, ::AbstractShiftingTechnique, particle) + return extract_svector(system.cache.delta_v, system, particle) +end + +function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi) + return system +end + +# Additional term in the momentum equation due to the shifting technique +@inline function dv_shifting(shifting, system, neighbor_system, + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + return zero(grad_kernel) +end + +@doc raw""" + ParticleShiftingTechnique() + +Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). +Following the original paper, the callback is applied in every time step and not +in every stage of a multi-stage time integration method to reduce the computational +cost and improve the stability of the scheme. + +See [Particle Shifting Technique](@ref shifting) for more information on the method. + +!!! warning + The Particle Shifting Technique needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** +""" +struct ParticleShiftingTechnique <: AbstractShiftingTechnique end + +# Zero because PST is applied in a callback +@inline function delta_v(system, ::ParticleShiftingTechnique, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +function particle_shifting_from_callback!(u_ode, shifting::ParticleShiftingTechnique, + system, v_ode, semi, dt) + @trixi_timeit timer() "particle shifting" begin + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + # Update the shifting velocity + update_shifting_from_callback!(system, shifting, v, u, v_ode, u_ode, semi) + + # Update the particle positions with the shifting velocity + particle_shifting!(u_ode, shifting, system, semi, dt) + end +end + +function update_shifting_from_callback!(system, ::ParticleShiftingTechnique, + v, u, v_ode, u_ode, semi) + (; cache) = system + (; delta_v) = cache + + set_zero!(delta_v) + + # This has similar performance to `maximum(..., eachparticle(system))`, + # but is GPU-compatible. + v_max = maximum(x -> sqrt(dot(x, x)), + reinterpret(reshape, SVector{ndims(system), eltype(v)}, + current_velocity(v, system))) + + # TODO this needs to be adapted to multi-resolution. + # Section 3.2 explains what else needs to be changed. + dx = particle_spacing(system, 1) + Wdx = smoothing_kernel(system, dx, 1) + h = smoothing_length(system, 1) + + foreach_system(semi) do neighbor_system + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = current_density(v, system, particle) + rho_b = current_density(v_neighbor, neighbor_system, neighbor) + + kernel = smoothing_kernel(system, distance, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # According to p. 29 below Eq. 9 + R = 2 // 10 + n = 4 + + # Eq. 7 in Sun et al. (2017). + # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h + # (see p. 29, right above Eq. 9), but this does not work when scaling h. + # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected + # for both small and large smoothing length factors. + # We need to scale + # - quadratically with the smoothing length, + # - linearly with the particle spacing, + # - linearly with the time step. + # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. + delta_v_ = -v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * + m_b / (rho_a + rho_b) * grad_kernel + + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end + + return system +end + +function particle_shifting!(u_ode, ::ParticleShiftingTechnique, system, semi, dt) + (; cache) = system + (; delta_v) = cache + + u = wrap_u(u_ode, system, semi) + + # Add δr from the cache to the current coordinates + @threaded semi for particle in eachparticle(system) + for i in axes(delta_v, 1) + @inbounds u[i, particle] += dt * delta_v[i, particle] + end + end + + return u +end + +""" + TransportVelocityAdami(background_pressure::Real) + +Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013) +to suppress pairing and tensile instability. +See [TVF](@ref transport_velocity_formulation) for more details of the method. + +# Arguments +- `background_pressure`: Background pressure. Suggested is a background pressure which is + on the order of the reference pressure. + +!!! warning + The Transport Velocity Formulation needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** +""" +struct TransportVelocityAdami{T <: Real} <: AbstractShiftingTechnique + background_pressure::T +end + +@inline function dv_shifting(::TransportVelocityAdami, system, neighbor_system, + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + v_a = current_velocity(v_system, system, particle) + delta_v_a = delta_v(system, particle) + + v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + delta_v_b = delta_v(neighbor_system, neighbor) + + A_a = rho_a * v_a * delta_v_a' + A_b = rho_b * v_b * delta_v_b' + + # The following term depends on the pressure acceleration formulation. + # See the large comment below. In the original paper (Adami et al., 2013), this is + # (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab. + # With the most common pressure acceleration formulation, this is + # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. + # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the + # `pressure_acceleration` function. + return pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) +end + +function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode, + u_ode, semi) + (; cache, correction) = system + (; delta_v) = cache + (; background_pressure) = shifting + + sound_speed = system_sound_speed(system) + + set_zero!(delta_v) + + foreach_system(semi) do neighbor_system + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_a = @inbounds hydrodynamic_mass(system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + rho_a = @inbounds current_density(v, system, particle) + rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor) + + h = smoothing_length(system, particle) + + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # In the original paper (Adami et al., 2013), the transport velocity is applied + # as follows: + # v_{1/2} = v_0 + Δt/2 * a, + # where a is the regular SPH acceleration term (pressure, viscosity, etc.). + # r_1 = r_0 + Δt * (v_{1/2}, + # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ] + # is the transport velocity. + # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity. + # We will call δv_{1/2} the shifting velocity, which is given by + # δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ], + # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. + # This term depends on the pressure acceleration formulation. + # In Zhang et al. (2017), the pressure acceleration term + # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab + # is used. They consequently changed the shifting velocity to + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ]. + # We therefore use the function `pressure_acceleration` to compute the + # shifting velocity according to the used pressure acceleration formulation. + # In most cases, this will be + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. + # + # In these papers, the shifting velocity is scaled by the time step Δt. + # We generally want the spatial discretization to be independent of the time step. + # Scaling the shifting velocity by the time step would lead to less shifting + # when very small time steps are used for testing/debugging purposes. + # This is especially problematic in TrixiParticles.jl, as the time step can vary + # significantly between different time integration methods (low vs high order). + # In order to eliminate the time step from the shifting velocity, we apply the + # CFL condition used in Adami et al. (2013): + # Δt <= 0.25 * h / c, + # where h is the smoothing length and c is the sound speed. + # Applying this equation as equality yields the shifting velocity + # δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. + # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the + # `pressure_acceleration` function. + delta_v_ = background_pressure / 8 * h / sound_speed * + pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) + + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end + + return system +end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl deleted file mode 100644 index fce7375dd5..0000000000 --- a/src/schemes/fluid/transport_velocity.jl +++ /dev/null @@ -1,158 +0,0 @@ -""" - TransportVelocityAdami(background_pressure::Real) - -Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013) -to suppress pairing and tensile instability. -See [TVF](@ref transport_velocity_formulation) for more details of the method. - -# Arguments -- `background_pressure`: Background pressure. Suggested is a background pressure which is - on the order of the reference pressure. -""" -struct TransportVelocityAdami{T <: Real} - background_pressure::T -end - -# No TVF for a system by default -@inline transport_velocity(system) = nothing - -create_cache_tvf(initial_condition, ::Nothing) = (;) - -function create_cache_tvf(initial_condition, ::TransportVelocityAdami) - delta_v = zeros(eltype(initial_condition), ndims(initial_condition), - nparticles(initial_condition)) - - return (; delta_v) -end - -# `δv` is the correction to the particle velocity due to the TVF. -# Particles are advected with the velocity `v + δv`. -@propagate_inbounds function delta_v(system, particle) - return delta_v(system, transport_velocity(system), particle) -end - -@propagate_inbounds function delta_v(system, ::TransportVelocityAdami, particle) - return extract_svector(system.cache.delta_v, system, particle) -end - -# Zero when no TVF is used -@inline function delta_v(system, transport_velocity, particle) - return zero(SVector{ndims(system), eltype(system)}) -end - -@inline function dv_transport_velocity(::Nothing, system, neighbor_system, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - return zero(grad_kernel) -end - -@inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - v_a = current_velocity(v_system, system, particle) - delta_v_a = delta_v(system, particle) - - v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - delta_v_b = delta_v(neighbor_system, neighbor) - - A_a = rho_a * v_a * delta_v_a' - A_b = rho_b * v_b * delta_v_b' - - # The following term depends on the pressure acceleration formulation. - # See the large comment below. In the original paper (Adami et al., 2013), this is - # (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab. - # With the most common pressure acceleration formulation, this is - # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. - # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the - # `pressure_acceleration` function. - return pressure_acceleration(system, neighbor_system, particle, neighbor, - m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) -end - -function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t) - return system -end - -function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode, - u_ode, semi, t) - (; cache, correction) = system - (; delta_v) = cache - (; background_pressure) = transport_velocity - - sound_speed = system_sound_speed(system) - - set_zero!(delta_v) - - foreach_system(semi) do neighbor_system - v_neighbor = wrap_v(v_ode, neighbor_system, semi) - u_neighbor = wrap_u(u_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - m_a = @inbounds hydrodynamic_mass(system, particle) - m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) - - rho_a = @inbounds current_density(v, system, particle) - rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor) - - h = smoothing_length(system, particle) - - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # In the original paper (Adami et al., 2013), the transport velocity is applied - # as follows: - # v_{1/2} = v_0 + Δt/2 * a, - # where a is the regular SPH acceleration term (pressure, viscosity, etc.). - # r_1 = r_0 + Δt * (v_{1/2}, - # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ] - # is the transport velocity. - # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity. - # We will call δv_{1/2} the shifting velocity, which is given by - # δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ], - # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. - # This term depends on the pressure acceleration formulation. - # In Zhang et al. (2017), the pressure acceleration term - # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab - # is used. They consequently changed the shifting velocity to - # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ]. - # We therefore use the function `pressure_acceleration` to compute the - # shifting velocity according to the used pressure acceleration formulation. - # In most cases, this will be - # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. - # - # In these papers, the shifting velocity is scaled by the time step Δt. - # We generally want the spatial discretization to be independent of the time step. - # Scaling the shifting velocity by the time step would lead to less shifting - # when very small time steps are used for testing/debugging purposes. - # This is especially problematic in TrixiParticles.jl, as the time step can vary - # significantly between different time integration methods (low vs high order). - # In order to eliminate the time step from the shifting velocity, we apply the - # CFL condition used in Adami et al. (2013): - # Δt <= 0.25 * h / c, - # where h is the smoothing length and c is the sound speed. - # Applying this equation as equality yields the shifting velocity - # δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. - # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the - # `pressure_acceleration` function. - delta_v_ = background_pressure / 8 * h / sound_speed * - pressure_acceleration(system, neighbor_system, particle, neighbor, - m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) - - # Write into the buffer - for i in eachindex(delta_v_) - @inbounds delta_v[i, particle] += delta_v_[i] - end - end - end - - return system -end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index af088e0a1b..b9cd27f079 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -70,13 +70,12 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Add convection term (only when using `TransportVelocityAdami`) - dv_tvf = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = dv_shifting(shifting_technique(particle_system), + particle_system, neighbor_system, particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index a5f56499b8..dda2ead5f1 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -5,7 +5,7 @@ acceleration=ntuple(_ -> 0.0, NDIMS), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -36,8 +36,9 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. density calculator and the correction method. To use [Tensile Instability Control](@ref tic), pass [`tensile_instability_control`](@ref) here. -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -59,7 +60,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. - `color_value`: The value used to initialize the color of particles in the system. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, - PF, TV, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} + PF, SC, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -71,7 +72,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: SC source_terms :: ST surface_tension :: SRFT surface_normal_method :: SRFN @@ -89,7 +90,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, ndims(smoothing_kernel)), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -147,7 +148,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, smoothing_length)..., - create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., color=Int(color_value)) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` @@ -162,7 +163,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, - transport_velocity, source_terms, surface_tension, + shifting_technique, source_terms, surface_tension, surface_normal_method, buffer, particle_refinement, cache) end @@ -177,6 +178,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) + print(io, ", ", system.shifting_technique) print(io, ", ", system.surface_tension) print(io, ", ", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal @@ -207,9 +209,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "state equation", system.state_equation |> typeof |> nameof) summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) summary_line(io, "density diffusion", system.density_diffusion) + summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "surface tension", system.surface_tension) summary_line(io, "surface normal method", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal @@ -278,7 +279,7 @@ end @inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed -@inline transport_velocity(system::WeaklyCompressibleSPHSystem) = system.transport_velocity +@inline shifting_technique(system::WeaklyCompressibleSPHSystem) = system.shifting_technique function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) @@ -316,7 +317,7 @@ function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, # Surface normal of neighbor and boundary needs to have been calculated already compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index ede9e6bbf9..a84a69317d 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -273,7 +273,19 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback()) + shifting_technique=ParticleShiftingTechnique(), + extra_callback=UpdateCallback()) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/periodic_channel_2d.jl with TVF" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "periodic_channel_2d.jl"), + tspan=(0.0, 0.2), + shifting_technique=TransportVelocityAdami(50_000.0), + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @@ -283,8 +295,9 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback(), - pressure_acceleration=tensile_instability_control) + shifting_technique=ParticleShiftingTechnique(), + pressure_acceleration=tensile_instability_control, + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 3d8cf9babb..ed37b9ea55 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -32,7 +32,7 @@ @test system.mass == mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -87,7 +87,7 @@ @test system.mass == setup.mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -140,7 +140,7 @@ │ viscosity: …………………………………………………… Nothing │ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ - │ tansport velocity formulation: Nothing │ + │ shifting technique: …………………………… nothing │ │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ @@ -225,7 +225,7 @@ names = ["No TVF", "TransportVelocityAdami"] @testset "$(names[i])" for i in eachindex(transport_velocity) system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, - transport_velocity=transport_velocity[i], + shifting_technique=transport_velocity[i], average_pressure_reduction=true, smoothing_length, 0.0) semi = Semidiscretization(system) diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index cf6020d068..99f7152db6 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -199,7 +199,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -211,8 +211,8 @@ │ state equation: ……………………………………… Val │ │ smoothing kernel: ………………………………… Val │ │ viscosity: …………………………………………………… nothing │ - │ tansport velocity formulation: Nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ + │ shifting technique: …………………………… nothing │ │ surface tension: …………………………………… nothing │ │ surface normal method: …………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │