Skip to content

Commit

Permalink
Add tests for wcsph-wcsph pressure_acceleration and interact! (#293)
Browse files Browse the repository at this point in the history 
* Change arguments of `pressure_acceleration`

* Add tests for `pressure_acceleration`

* Move `using` to `test_util.jl`

* Remove redundant unpack

* Add momentum conservation tests for fluid-fluid interaction

* Fix typo

* Add source for momentum conservation

* Implement suggestions

* Change density and pressure, so that the factor `p/rho^2` isn't too large and increases the error
  • Loading branch information
@efaulhaber
efaulhaber authored Nov 22, 2023
1 parent fe6eedc commit 2307c50
Show file tree
Hide file tree
Showing 8 changed files with 207 additions and 62 deletions.
52 changes: 17 additions & 35 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,18 +203,16 @@ function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing
return zero(initial_density)
end

@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system,
v_boundary_system,
@inline function pressure_acceleration(pressure_correction, m_b,
particle, boundary_particle,
particle_system, boundary_system,
boundary_model::BoundaryModelDummyParticles,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator)
(; density_calculator) = boundary_model

pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system, v_boundary_system,
pressure_acceleration(pressure_correction, m_b, particle, boundary_particle,
particle_system, boundary_system,
boundary_model, density_calculator,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator)
Expand All @@ -223,17 +221,13 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with ContinuityDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system,
v_boundary_system,
@inline function pressure_acceleration(pressure_correction, m_b,
particle, boundary_particle,
particle_system, boundary_system,
boundary_model::BoundaryModelDummyParticles,
boundary_density_calculator,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator::ContinuityDensity)
rho_a = particle_density(v_particle_system, particle_system, particle)
rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)

return -m_b *
(particle_system.pressure[particle] + boundary_model.pressure[boundary_particle]) /
(rho_a * rho_b) * grad_kernel
Expand All @@ -242,17 +236,13 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with SummationDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system,
v_boundary_system,
@inline function pressure_acceleration(pressure_correction, m_b,
particle, boundary_particle,
particle_system, boundary_system,
boundary_model::BoundaryModelDummyParticles,
boundary_density_calculator,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator::SummationDensity)
rho_a = particle_density(v_particle_system, particle_system, particle)
rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)

return -m_b *
(particle_system.pressure[particle] / rho_a^2 +
boundary_model.pressure[boundary_particle] / rho_b^2) *
Expand All @@ -262,17 +252,13 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with ContinuityDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system,
v_boundary_system,
@inline function pressure_acceleration(pressure_correction, m_b,
particle, boundary_particle,
particle_system, boundary_system,
boundary_model::BoundaryModelDummyParticles,
::PressureMirroring,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator::ContinuityDensity)
rho_a = particle_density(v_particle_system, particle_system, particle)
rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)

return -m_b *
(particle_system.pressure[particle] + particle_system.pressure[particle]) /
(rho_a * rho_b) * grad_kernel
Expand All @@ -281,17 +267,13 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with SummationDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, boundary_particle,
boundary_system,
v_boundary_system,
@inline function pressure_acceleration(pressure_correction, m_b,
particle, boundary_particle,
particle_system, boundary_system,
boundary_model::BoundaryModelDummyParticles,
::PressureMirroring,
rho_a, rho_b, pos_diff, distance, grad_kernel,
fluid_density_calculator::SummationDensity)
rho_a = particle_density(v_particle_system, particle_system, particle)
rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)

return -m_b *
(particle_system.pressure[particle] / rho_a^2 +
particle_system.pressure[particle] / rho_b^2) *
Expand Down
5 changes: 2 additions & 3 deletions src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,9 +76,8 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar)
print(io, ")")
end

@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor, neighbor_system,
v_neighbor_system,
@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system, neighbor_system,
boundary_model::BoundaryModelMonaghanKajtar, rho_a,
rho_b, pos_diff, distance, grad_kernel,
density_calculator)
Expand Down
37 changes: 18 additions & 19 deletions src/schemes/fluid/weakly_compressible_sph/rhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,15 @@ function interact!(dv, v_particle_system, u_particle_system,
v_neighbor_system, u_neighbor_system, neighborhood_search,
particle_system::WeaklyCompressibleSPHSystem,
neighbor_system)
(; density_calculator, state_equation, correction, smoothing_length) = particle_system
(; density_calculator, state_equation, correction) = particle_system
(; sound_speed) = state_equation

viscosity = viscosity_model(neighbor_system)
system_coords = current_coordinates(u_particle_system, particle_system)
neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system)

# In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code
# and the two other lines below that are marked as "debug example".
# In order to visualize quantities like pressure forces or viscosity forces, uncomment
# the following code and the two other lines below that are marked as "debug example".
# debug_array = zeros(ndims(particle_system), nparticles(particle_system))

# Loop over all pairs of particles and neighbors within the kernel cutoff.
Expand All @@ -38,10 +38,9 @@ function interact!(dv, v_particle_system, u_particle_system,
m_a = hydrodynamic_mass(particle_system, particle)
m_b = hydrodynamic_mass(neighbor_system, neighbor)

dv_pressure = pressure_acceleration(pressure_correction, m_b, particle,
particle_system, v_particle_system,
neighbor, neighbor_system,
v_neighbor_system, rho_a, rho_b, pos_diff,
dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system, neighbor_system,
rho_a, rho_b, pos_diff,
distance, grad_kernel, density_calculator)

dv_viscosity = viscosity_correction *
Expand Down Expand Up @@ -74,10 +73,10 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with ContinuityDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor,
@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system,
neighbor_system::WeaklyCompressibleSPHSystem,
v_neighbor_system, rho_a, rho_b, pos_diff, distance,
rho_a, rho_b, pos_diff, distance,
grad_kernel, ::ContinuityDensity)
return (-m_b *
(particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) /
Expand All @@ -88,28 +87,28 @@ end
# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
# used with SummationDensity.
@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor,
@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system,
neighbor_system::WeaklyCompressibleSPHSystem,
v_neighbor_system, rho_a, rho_b, pos_diff, distance,
rho_a, rho_b, pos_diff, distance,
grad_kernel, ::SummationDensity)
return (-m_b *
(particle_system.pressure[particle] / rho_a^2 +
neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) *
pressure_correction
end

@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor,
@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system,
neighbor_system::Union{BoundarySPHSystem,
TotalLagrangianSPHSystem},
v_neighbor_system, rho_a, rho_b, pos_diff, distance,
rho_a, rho_b, pos_diff, distance,
grad_kernel, density_calculator)
(; boundary_model) = neighbor_system

return pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor, neighbor_system,
v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff,
return pressure_acceleration(pressure_correction, m_b, particle, neighbor,
particle_system, neighbor_system,
boundary_model, rho_a, rho_b, pos_diff,
distance, grad_kernel, density_calculator)
end

Expand Down
5 changes: 2 additions & 3 deletions src/schemes/solid/total_lagrangian_sph/rhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,9 +98,8 @@ function interact!(dv, v_particle_system, u_particle_system,

# Boundary forces
# Note: neighbor and particle are switched in this call and `pressure_correction` is set to `1.0` (no correction)
dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system,
v_neighbor_system, particle,
particle_system, v_particle_system,
dv_boundary = pressure_acceleration(1.0, m_b, neighbor, particle,
neighbor_system, particle_system,
boundary_model, rho_a,
rho_b, pos_diff, distance,
grad_kernel, density_calculator)
Expand Down
2 changes: 0 additions & 2 deletions test/general/smoothing_kernels.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
using QuadGK

@testset verbose=true "Smoothing Kernels" begin
# Don't show all kernel tests in the final overview
@testset verbose=false "Integral" begin
Expand Down
Loading

0 comments on commit 2307c50

Please sign in to comment.