Skip to content

Commit ff5403a

Browse files
efaulhaberefaulhaber
committed
Make pressure clipping in BoundaryModelDummyParticles optional
1 parent 379db88 commit ff5403a

File tree

2 files changed

+122
-33
lines changed

2 files changed

+122
-33
lines changed

‎src/schemes/boundary/wall_boundary/dummy_particles.jl‎

Lines changed: 65 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
density_calculator, smoothing_kernel,
44
smoothing_length; viscosity=nothing,
55
state_equation=nothing, correction=nothing,
6+
clip_negative_pressure=true,
67
reference_particle_spacing=0.0)
78
89
Boundary model for [`WallBoundarySystem`](@ref).
@@ -22,6 +23,15 @@ Boundary model for [`WallBoundarySystem`](@ref).
2223
- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)).
2324
- `viscosity`: Slip (default) or no-slip condition. See description below for further
2425
information.
26+
- `clip_negative_pressure=true`: Clip negative boundary pressures to avoid sticking
27+
artifacts from attractive fluid-boundary forces at free
28+
surfaces. Note that this is not a correct formulation
29+
at interior boundaries (away from free surfaces).
30+
Disable this option when simulating closed systems without
31+
free surfaces to avoid artificially increased boundary
32+
pressures that cause larger gaps between fluid and boundary
33+
in areas of low pressure, against which the particle
34+
shifting technique is fighting.
2535
- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
2636
which currently is only needed when using surface tension.
2737
# Examples
@@ -39,24 +49,38 @@ boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExt
3949
BoundaryModelDummyParticles(AdamiPressureExtrapolation, ViscosityAdami)
4050
```
4151
"""
42-
struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C}
43-
pressure :: VECTOR # Vector{ELTYPE}
44-
hydrodynamic_mass :: VECTOR # Vector{ELTYPE}
45-
state_equation :: SE
46-
density_calculator :: DC
47-
smoothing_kernel :: K
48-
smoothing_length :: ELTYPE
49-
viscosity :: V
50-
correction :: COR
51-
cache :: C
52-
end
53-
54-
# The default constructor needs to be accessible for Adapt.jl to work with this struct.
55-
# See the comments in general/gpu.jl for more details.
52+
struct BoundaryModelDummyParticles{DC, SE, CLIP, ELTYPE <: Real, VECTOR, K, V, COR, C}
53+
pressure :: VECTOR # Vector{ELTYPE}
54+
hydrodynamic_mass :: VECTOR # Vector{ELTYPE}
55+
state_equation :: SE
56+
density_calculator :: DC
57+
smoothing_kernel :: K
58+
smoothing_length :: ELTYPE
59+
viscosity :: V
60+
correction :: COR
61+
cache :: C
62+
# Store this both as field and type parameter to avoid annoying hand-written
63+
# `adapt_structure` functions.
64+
clip_negative_pressure :: Bool
65+
66+
function BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation,
67+
density_calculator, smoothing_kernel,
68+
smoothing_length, viscosity, correction,
69+
cache, clip_negative_pressure)
70+
return new{typeof(density_calculator), typeof(state_equation),
71+
clip_negative_pressure, eltype(pressure), typeof(pressure),
72+
typeof(smoothing_kernel), typeof(viscosity), typeof(correction),
73+
typeof(cache)}(pressure, hydrodynamic_mass, state_equation,
74+
density_calculator, smoothing_kernel, smoothing_length,
75+
viscosity, correction, cache, clip_negative_pressure)
76+
end
77+
end
78+
5679
function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
5780
density_calculator, smoothing_kernel,
5881
smoothing_length; viscosity=nothing,
5982
state_equation=nothing, correction=nothing,
83+
clip_negative_pressure=true,
6084
reference_particle_spacing=0.0)
6185
pressure = initial_boundary_pressure(initial_density, density_calculator,
6286
state_equation)
@@ -82,10 +106,17 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
82106

83107
return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation,
84108
density_calculator, smoothing_kernel,
85-
smoothing_length, viscosity, correction, cache)
109+
smoothing_length, viscosity, correction, cache,
110+
clip_negative_pressure)
111+
end
112+
113+
@inline function Base.ndims(boundary_model::BoundaryModelDummyParticles)
114+
return ndims(boundary_model.smoothing_kernel)
86115
end
87116

88-
@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel)
117+
@inline function clip_negative_pressure(::BoundaryModelDummyParticles{<:Any, <:Any, CLIP}) where {CLIP}
118+
return CLIP
119+
end
89120

90121
@doc raw"""
91122
AdamiPressureExtrapolation(; pressure_offset=0, allow_loop_flipping=true)
@@ -108,7 +139,6 @@ end
108139
loop is not thread parallelizable.
109140
This can cause error variations between simulations with
110141
different numbers of threads.
111-
112142
"""
113143
struct AdamiPressureExtrapolation{ELTYPE}
114144
pressure_offset :: ELTYPE
@@ -423,33 +453,33 @@ end
423453

424454
function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityDensity},
425455
system, v, u, v_ode, u_ode, semi)
426-
427-
# Limit pressure to be non-negative to avoid attractive forces between fluid and
428-
# boundary particles at free surfaces (sticking artifacts).
429456
@threaded semi for particle in eachparticle(system)
430-
apply_state_equation!(boundary_model, current_density(v, system, particle),
431-
particle)
457+
@inbounds apply_state_equation!(boundary_model,
458+
current_density(v, system, particle), particle)
432459
end
433460

434461
return boundary_model
435462
end
436463

437-
# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`).
438-
# Otherwise, `@threaded` does not work here with Julia ARM on macOS.
439-
# See https://github.com/JuliaSIMD/Polyester.jl/issues/88.
440-
@inline function apply_state_equation!(boundary_model, density, particle)
441-
boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0)
464+
@propagate_inbounds function apply_state_equation!(boundary_model, density, particle)
465+
pressure = boundary_model.state_equation(density)
466+
467+
# This is determined statically and has therefore no overhead.
468+
if clip_negative_pressure(boundary_model)
469+
# Clip pressure to avoid sticking artifacts from negative pressures at free surfaces.
470+
pressure = max(pressure, 0)
471+
end
472+
473+
boundary_model.pressure[particle] = pressure
442474
end
443475

444476
@propagate_inbounds function apply_state_equation!(boundary_model::BoundaryModelDummyParticles{<:Any,
445-
<:Any,
446-
<:Any,
447477
Nothing},
448478
density, particle)
449479
# Contact-only wall setups can reuse dummy particles for rigid contact without
450480
# configuring a hydrodynamic state equation. In that case, keep the auxiliary wall
451481
# pressure at zero because it is only meaningful for fluid-coupled updates.
452-
boundary_model.pressure[particle] = zero(eltype(boundary_model.pressure))
482+
boundary_model.pressure[particle] = 0
453483
end
454484

455485
function compute_pressure!(boundary_model,
@@ -518,9 +548,11 @@ function compute_adami_density!(boundary_model, system, v, particle)
518548
compute_wall_velocity!(viscosity, system, v, particle)
519549
end
520550

521-
# Limit pressure to be non-negative to avoid attractive forces between fluid and
522-
# boundary particles at free surfaces (sticking artifacts).
523-
@inbounds pressure[particle] = max(pressure[particle], 0)
551+
# Clip pressure to avoid sticking artifacts from negative pressures at free surfaces.
552+
# This is determined statically and has therefore no overhead.
553+
if clip_negative_pressure(boundary_model)
554+
@inbounds pressure[particle] = max(pressure[particle], 0)
555+
end
524556

525557
# Apply inverse state equation to compute density (not used with EDAC)
526558
inverse_state_equation!(density, state_equation, pressure, particle)

‎test/schemes/boundary/dummy_particles/dummy_particles.jl‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,63 @@
1010
@test repr(boundary_model) == expected_repr
1111
end
1212

13+
@testset "Pressure clipping" begin
14+
state_equation = StateEquationCole(sound_speed=10.0,
15+
reference_density=1000.0,
16+
exponent=7,
17+
clip_negative_pressure=false)
18+
smoothing_kernel = SchoenbergCubicSplineKernel{2}()
19+
smoothing_length = 0.1
20+
21+
boundary_model_clip = BoundaryModelDummyParticles([1000.0], [1.0],
22+
SummationDensity(),
23+
smoothing_kernel,
24+
smoothing_length,
25+
state_equation=state_equation)
26+
boundary_model_no_clip = BoundaryModelDummyParticles([1000.0], [1.0],
27+
SummationDensity(),
28+
smoothing_kernel,
29+
smoothing_length,
30+
state_equation=state_equation,
31+
clip_negative_pressure=false)
32+
33+
@test TrixiParticles.clip_negative_pressure(boundary_model_clip)
34+
@test !TrixiParticles.clip_negative_pressure(boundary_model_no_clip)
35+
36+
TrixiParticles.apply_state_equation!(boundary_model_clip, 900.0, 1)
37+
TrixiParticles.apply_state_equation!(boundary_model_no_clip, 900.0, 1)
38+
39+
@test boundary_model_clip.pressure[1] == 0
40+
@test boundary_model_no_clip.pressure[1] < 0
41+
42+
# Test that clipping behavior is preserved through `adapt_structure`.
43+
adapted_model_clip = TrixiParticles.Adapt.adapt(Array, boundary_model_clip)
44+
adapted_model_no_clip = TrixiParticles.Adapt.adapt(Array, boundary_model_no_clip)
45+
@test TrixiParticles.clip_negative_pressure(adapted_model_clip)
46+
@test !TrixiParticles.clip_negative_pressure(adapted_model_no_clip)
47+
48+
boundary_model_adami_clip = BoundaryModelDummyParticles([1000.0], [1.0],
49+
AdamiPressureExtrapolation(),
50+
smoothing_kernel,
51+
smoothing_length,
52+
state_equation=state_equation)
53+
boundary_model_adami_no_clip = BoundaryModelDummyParticles([1000.0], [1.0],
54+
AdamiPressureExtrapolation(),
55+
smoothing_kernel,
56+
smoothing_length,
57+
state_equation=state_equation,
58+
clip_negative_pressure=false)
59+
60+
for model in (boundary_model_adami_clip, boundary_model_adami_no_clip)
61+
model.cache.volume[1] = 1
62+
model.pressure[1] = -1
63+
TrixiParticles.compute_adami_density!(model, nothing, nothing, 1)
64+
end
65+
66+
@test boundary_model_adami_clip.pressure[1] == 0
67+
@test boundary_model_adami_no_clip.pressure[1] == -1
68+
end
69+
1370
@testset verbose=true "Viscosity Adami/Bernoulli: Wall Velocity" begin
1471
particle_spacing = 0.1
1572

0 commit comments

Comments
 (0)