Skip to content

Commit

Permalink
add new bnd density calculator
Browse files Browse the repository at this point in the history
  • Loading branch information
@svchb
svchb committed Jul 15, 2024
1 parent a6183e2 commit 0f3536c
Show file tree
Hide file tree
Showing 5 changed files with 106 additions and 46 deletions.
14 changes: 9 additions & 5 deletions docs/src/systems/boundary.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,24 +55,28 @@ of the boundary particle ``b``.

### Hydrodynamic density of dummy particles

We provide five options to compute the boundary density and pressure, determined by the `density_calculator`:
We provide six options to compute the boundary density and pressure, determined by the `density_calculator`:
1. (Recommended) With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the
fluid according to (Adami et al., 2012), and the density is obtained by applying the inverse of the state equation.
This option usually yields the best results of the options listed here.
2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
2. (Only relevant for FSI) With `BernoulliPressureExtrapolation`, the pressure is extrapolated from the
pressure similar to the `AdamiPressureExtrapolation`, but a relative velocity dependent pressure part
is calculated between moving solids and fluids, which increases the boundary pressure in areas prone to
penetrations.
3. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
and the pressure is computed from the density with the state equation.
3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
4. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
and the pressure is computed from the density with the state equation.
Note that this causes a gap between fluid and boundary where the boundary is initialized
without any contact to the fluid. This is due to overestimation of the boundary density
as soon as the fluid comes in contact with boundary particles that initially did not have
contact to the fluid.
Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat.
See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13).
4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure
5. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure
is computed from the density with the state equation.
This option is not recommended. The other options yield significantly better results.
5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure
6. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure
is not used. Instead, the fluid pressure is mirrored as boundary pressure in the
momentum equation.
This option is not recommended due to stability issues. See [`PressureMirroring`](@ref)
Expand Down
2 changes: 1 addition & 1 deletion examples/fsi/falling_spheres_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,

# ==========================================================================================
# ==== Boundary
boundary_density_calculator = AdamiPressureExtrapolation()
boundary_density_calculator = BernoulliPressureExtrapolation()
boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
state_equation=state_equation,
boundary_density_calculator,
Expand Down
2 changes: 1 addition & 1 deletion src/TrixiParticles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris
export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
DensityDiffusionAntuono
export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
PressureMirroring, PressureZeroing
PressureMirroring, PressureZeroing, BernoulliPressureExtrapolation
export BoundaryMovement
export examples_dir, validation_dir, trixi_include
export trixi2vtk
Expand Down
5 changes: 5 additions & 0 deletions src/general/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,11 @@ end

@inline current_velocity(v, system, particle) = extract_svector(v, system, particle)

@inline function current_acceleration(system, particle)
# TODO: Return `dv` of solid particles
return zero(SVector{ndims(system), eltype(system)})
end

@inline function smoothing_kernel(system, distance)
(; smoothing_kernel, smoothing_length) = system
return kernel(smoothing_kernel, distance, smoothing_length)
Expand Down
129 changes: 90 additions & 39 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,24 @@ struct AdamiPressureExtrapolation{ELTYPE}
end
end

@doc raw"""
BernoulliPressureExtrapolation(; pressure_offset=0.0)
`density_calculator` for `BoundaryModelDummyParticles`.
# Keywords
- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure
to prevent penetration, which is possible by increasing this value.
"""
struct BernoulliPressureExtrapolation{ELTYPE}
pressure_offset::ELTYPE

function BernoulliPressureExtrapolation(; pressure_offset=0.0)
return new{eltype(pressure_offset)}(pressure_offset)
end
end

@doc raw"""
PressureMirroring()
Expand Down Expand Up @@ -145,7 +163,9 @@ end

@inline create_cache_model(initial_density, ::ContinuityDensity) = (; initial_density)

function create_cache_model(initial_density, ::AdamiPressureExtrapolation)
function create_cache_model(initial_density,
::Union{AdamiPressureExtrapolation,
BernoulliPressureExtrapolation})
density = copy(initial_density)
volume = similar(initial_density)

Expand Down Expand Up @@ -205,7 +225,7 @@ end
# Note that the other density calculators are dispatched in `density_calculators.jl`
@inline function particle_density(v,
::Union{AdamiPressureExtrapolation, PressureMirroring,
PressureZeroing},
PressureZeroing, BernoulliPressureExtrapolation},
boundary_model, particle)
(; cache) = boundary_model

Expand All @@ -227,10 +247,11 @@ end

function compute_density!(boundary_model,
::Union{ContinuityDensity, AdamiPressureExtrapolation,
BernoulliPressureExtrapolation,
PressureMirroring, PressureZeroing},
system, v, u, v_ode, u_ode, semi)
# No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`.
# For `AdamiPressureExtrapolation`, the density is updated in `compute_pressure!`.
# For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`.
return boundary_model
end

Expand Down Expand Up @@ -310,7 +331,10 @@ end
boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0)
end

function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u,
function compute_pressure!(boundary_model,
bnd_density_calc::Union{AdamiPressureExtrapolation,
BernoulliPressureExtrapolation}, system,
v, u,
v_ode, u_ode, semi)
(; pressure, cache, viscosity) = boundary_model

Expand Down Expand Up @@ -338,16 +362,18 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system,
nhs = get_neighborhood_search(neighbor_system, system, semi)

# Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries
adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, nhs)
bnd_pressure_extrapolation_neighbor!(boundary_model, bnd_density_calc, system,
neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, nhs)
else
nhs = get_neighborhood_search(system, neighbor_system, semi)

# Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries
adami_pressure_extrapolation!(boundary_model, system, neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, nhs)
bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, system,
neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, nhs)
end

@threaded system for particle in eachparticle(system)
Expand Down Expand Up @@ -406,43 +432,43 @@ end
v_neighbor_system[1:ndims(system), neighbor])^2
end

@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system,
neighbor_system, system_coords,
neighbor_coords, v_neighbor_system,
neighborhood_search)
@inline function bnd_pressure_extrapolation_neighbor!(boundary_model, system,
neighbor_system, system_coords,
neighbor_coords, v_neighbor_system,
neighborhood_search)
return boundary_model
end

@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system,
neighbor_system::FluidSystem,
system_coords, neighbor_coords, v,
v_neighbor_system,
neighborhood_search)
@inline function bnd_pressure_extrapolation_neighbor!(boundary_model, bnd_density_calc,
system, neighbor_system::FluidSystem,
system_coords, neighbor_coords, v,
v_neighbor_system,
neighborhood_search)
(; pressure, cache, viscosity) = boundary_model
(; pressure_offset) = density_calculator


foreach_point_neighbor(neighbor_system, system, neighbor_coords, system_coords,
neighborhood_search; points=eachparticle(neighbor_system),
parallel=false) do neighbor, particle, pos_diff, distance
# Since neighbor and particle are switched
pos_diff = -pos_diff
adami_pressure_inner!(boundary_model, system, neighbor_system, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure, pressure_offset)
bnd_pressure_inner!(boundary_model, bnd_density_calc, system, neighbor_system, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure, pressure_offset)
end
end

@inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, neighborhood_search)
@inline function bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, system,
neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, neighborhood_search)
return boundary_model
end

@inline function adami_pressure_extrapolation!(boundary_model, system,
neighbor_system::FluidSystem,
system_coords, neighbor_coords, v,
v_neighbor_system, neighborhood_search)
@inline function bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, system,
neighbor_system::FluidSystem,
system_coords, neighbor_coords, v,
v_neighbor_system, neighborhood_search)
(; pressure, cache, viscosity, density_calculator) = boundary_model
(; pressure_offset) = density_calculator

Expand All @@ -451,19 +477,21 @@ end
neighborhood_search;
points=eachparticle(system)) do particle, neighbor,
pos_diff, distance
adami_pressure_inner!(boundary_model, system, neighbor_system, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure, pressure_offset)
bnd_pressure_inner!(boundary_model, bnd_density_calc, system, neighbor_system, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure, pressure_offset)
end
end

@inline function adami_pressure_inner!(boundary_model, system,
neighbor_system::FluidSystem, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure,
pressure_offset)
@inline function bnd_pressure_inner!(boundary_model, ::BernoulliPressureExtrapolation,
system, neighbor_system::FluidSystem, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure,
pressure_offset)
density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor)

resulting_acc = neighbor_system.acceleration - current_acceleration(system, particle)

kernel_weight = smoothing_kernel(boundary_model, distance)

pressure[particle] += (pressure_offset +
Expand All @@ -472,7 +500,30 @@ end
dynamic_pressure(density_neighbor, v, v_neighbor_system,
particle, neighbor, system)
+
dot(neighbor_system.acceleration,
dot(resulting_acc,
density_neighbor * pos_diff)) * kernel_weight

cache.volume[particle] += kernel_weight

compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
kernel_weight, particle, neighbor)
end

@inline function bnd_pressure_inner!(boundary_model, ::AdamiPressureExtrapolation,
system, neighbor_system::FluidSystem, v,
v_neighbor_system, particle, neighbor, pos_diff,
distance, viscosity, cache, pressure,
pressure_offset)
density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor)

resulting_acc = neighbor_system.acceleration - current_acceleration(system, particle)

kernel_weight = smoothing_kernel(boundary_model, distance)

pressure[particle] += (pressure_offset +
particle_pressure(v_neighbor_system, neighbor_system,
neighbor) +
dot(resulting_acc,
density_neighbor * pos_diff)) * kernel_weight

cache.volume[particle] += kernel_weight
Expand Down

0 comments on commit 0f3536c

Please sign in to comment.