Skip to content

Commit

Permalink
Merge branch 'main'
Browse files Browse the repository at this point in the history 
Conflicts:
	src/schemes/boundary/dummy_particles/dummy_particles.jl
	src/schemes/fluid/weakly_compressible_sph/rhs.jl
  • Loading branch information
@svchb
svchb committed Nov 24, 2023
2 parents ad42394 + 2307c50 commit 8612c82
Show file tree
Hide file tree
Showing 10 changed files with 233 additions and 74 deletions.
25 changes: 13 additions & 12 deletions examples/fluid/dam_break_2d_corrections.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,11 +54,12 @@ smoothing_kernel_dict = Dict(

begin
mod = trixi_include_safe(joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length,
boundary_density_calculator=ContinuityDensity(),
fluid_density_calculator=ContinuityDensity(),
correction=Nothing(), use_reinit=true,
prefix="continuity_reinit", tspan=tspan)
fluid_particle_spacing=particle_spacing,
smoothing_length=smoothing_length,
boundary_density_calculator=ContinuityDensity(),
fluid_density_calculator=ContinuityDensity(),
correction=Nothing(), use_reinit=true,
prefix="continuity_reinit", tspan=tspan)
end

# Clip negative pressure to be able to use `SummationDensity`
Expand All @@ -75,11 +76,11 @@ for correction_name in keys(correction_dict)
println("fluid/dam_break_2d.jl with ", correction_name)

trixi_include_safe(joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
fluid_particle_spacing=particle_spacing,
smoothing_length=smoothing_length,
boundary_density_calculator=boundary_density_calculator,
fluid_density_calculator=fluid_density_calculator,
correction=correction, use_reinit=false,
state_equation=state_equation, smoothing_kernel=smoothing_kernel,
prefix="$(correction_name)", tspan=tspan)
fluid_particle_spacing=particle_spacing,
smoothing_length=smoothing_length,
boundary_density_calculator=boundary_density_calculator,
fluid_density_calculator=fluid_density_calculator,
correction=correction, use_reinit=false,
state_equation=state_equation, smoothing_kernel=smoothing_kernel,
prefix="$(correction_name)", tspan=tspan)
end
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, correction)
(; 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, correction)
Expand Down Expand Up @@ -242,18 +240,14 @@ 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,
correction)
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 @@ -262,18 +256,14 @@ 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,
correction)
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 @@ -283,18 +273,14 @@ 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,
correction)
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 @@ -303,18 +289,14 @@ 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,
correction)
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, correction)
Expand Down
49 changes: 31 additions & 18 deletions src/schemes/fluid/weakly_compressible_sph/rhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ function interact!(dv, v_particle_system, u_particle_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 @@ -39,10 +39,9 @@ function interact!(dv, v_particle_system, u_particle_system,
m_b = hydrodynamic_mass(neighbor_system, neighbor)

# if variable smoothing_length is used this uses the particle smoothing length
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,
correction)

Expand Down Expand Up @@ -77,10 +76,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, correction)
return (-m_b *
(particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) /
Expand All @@ -91,10 +90,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 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, correction)
return (-m_b *
(particle_system.pressure[particle] / rho_a^2 +
Expand All @@ -116,17 +115,31 @@ end
# pressure_correction
# end

@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
v_particle_system, neighbor,
# # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
# # formulations" by Bonet and Lok (1999), for a corrected gradient of the smoothing kernel
# # this simplifies to the formulation below.
# @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
# v_particle_system, neighbor,
# neighbor_system::WeaklyCompressibleSPHSystem,
# v_neighbor_system, rho_a, rho_b, pos_diff, distance,
# grad_kernel, density_calculator,
# ::Union{KernelGradientCorrection,
# MixedKernelGradientCorrection})
# return -m_b / rho_b * neighbor_system.pressure[neighbor] * grad_kernel *
# pressure_correction
# end

@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, correction)
(; 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, correction)
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 @@ -99,9 +99,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
1 change: 0 additions & 1 deletion src/util.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,6 @@ function trixi_include_safe(file_path; kwargs...)
return isolated_mod
end


# Helper methods used in the functions defined above, also copied from Trixi.jl

# Apply the function `f` to `expr` and all sub-expressions recursively.
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 8612c82

Please sign in to comment.