Skip to content

Fix Adami boundary viscosity #853

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
LasNikas merged 6 commits into from
Jul 23, 2025
Merged

Fix Adami boundary viscosity #853

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
63ebd28
Fix Adami boundary viscosity
efaulhaber Jul 1, 2025
94a50cb
Merge branch 'main' into fix-adami-boundary-viscosity
efaulhaber Jul 1, 2025
3de404e
Merge branch 'main' into fix-adami-boundary-viscosity
efaulhaber Jul 22, 2025
b915fff
Add tests
efaulhaber Jul 22, 2025
52e114c
Reformat
efaulhaber Jul 22, 2025
fcf780c
Merge branch 'main' into fix-adami-boundary-viscosity
svchb Jul 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -427,14 +427,14 @@ function compute_pressure!(boundary_model,

@trixi_timeit timer() "inverse state equation" @threaded semi for particle in
eachparticle(system)
compute_adami_density!(boundary_model, system, system_coords, particle)
compute_adami_density!(boundary_model, system, v, particle)
end
end

# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`).
# Otherwise, `@threaded` does not work here with Julia ARM on macOS.
# See https://github.com/JuliaSIMD/Polyester.jl/issues/88.
function compute_adami_density!(boundary_model, system, system_coords, particle)
function compute_adami_density!(boundary_model, system, v, particle)
(; pressure, state_equation, cache, viscosity) = boundary_model
(; volume, density) = cache

Expand All @@ -445,7 +445,7 @@ function compute_adami_density!(boundary_model, system, system_coords, particle)
@inbounds pressure[particle] /= volume[particle]

# To impose no-slip condition
compute_wall_velocity!(viscosity, system, system_coords, particle)
compute_wall_velocity!(viscosity, system, v, particle)
end

# Limit pressure to be non-negative to avoid attractive forces between fluid and
Expand Down Expand Up @@ -599,14 +599,14 @@ end
return viscosity
end

@inline function compute_wall_velocity!(viscosity, system, system_coords, particle)
@inline function compute_wall_velocity!(viscosity, system, v, particle)
(; boundary_model) = system
(; cache) = boundary_model
(; volume, wall_velocity) = cache

# Prescribed velocity of the boundary particle.
# This velocity is zero when not using moving boundaries.
v_boundary = current_velocity(system_coords, system, particle)
v_boundary = current_velocity(v, system, particle)

for dim in eachindex(v_boundary)
# The second term is the precalculated smoothed velocity field of the fluid
Expand Down
107 changes: 65 additions & 42 deletions test/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
@test repr(boundary_model) == expected_repr
end

@testset "Viscosity Adami/Bernoulli: Wall Velocity" begin
@testset verbose=true "Viscosity Adami/Bernoulli: Wall Velocity" begin
particle_spacing = 0.1

# Boundary particles in fluid compact support
Expand All @@ -30,6 +30,7 @@
density=257.0)

boundary = union(boundary_1, boundary_2, boundary_3)
v_boundary = zeros(2, TrixiParticles.nparticles(boundary))

particles_in_compact_support = length(boundary_1.mass) + length(boundary_2.mass)

Expand All @@ -49,62 +50,82 @@
exponent=7)

# Define pressure extrapolation methods to test
pressure_extrapolations = [
AdamiPressureExtrapolation(),
BernoulliPressureExtrapolation()
boundary_model_adami = BoundaryModelDummyParticles(boundary.density,
boundary.mass,
state_equation=state_equation,
AdamiPressureExtrapolation(),
smoothing_kernel,
smoothing_length,
viscosity=viscosity)
boundary_model_bernoulli = BoundaryModelDummyParticles(boundary.density,
Comment thread
svchb marked this conversation as resolved.
boundary.mass,
state_equation=state_equation,
BernoulliPressureExtrapolation(),
smoothing_kernel,
smoothing_length,
viscosity=viscosity)

boundary_systems = [
BoundarySPHSystem(boundary, boundary_model_adami),
BoundarySPHSystem(boundary, boundary_model_bernoulli),
TotalLagrangianSPHSystem(boundary, smoothing_kernel,
smoothing_length, 1e6, 0.3;
boundary_model=boundary_model_adami),
TotalLagrangianSPHSystem(boundary, smoothing_kernel,
smoothing_length, 1e6, 0.3;
boundary_model=boundary_model_bernoulli)
]

semi = DummySemidiscretization()
# Create fluid system
fluid_system = WeaklyCompressibleSPHSystem(fluid,
SummationDensity(),
state_equation,
smoothing_kernel,
smoothing_length)

velocities = [
[0.0; -1.0],
[1.0; 1.0],
[-1.0; 0.0],
[0.7; 0.2],
[0.3; 0.8]
]

for pressure_extrapolation in pressure_extrapolations
@testset "Pressure Extrapolation: $(typeof(pressure_extrapolation))" begin
# Create boundary and fluid systems
boundary_model = BoundaryModelDummyParticles(boundary.density,
boundary.mass,
state_equation=state_equation,
pressure_extrapolation,
smoothing_kernel,
smoothing_length,
viscosity=viscosity)
boundary_system = BoundarySPHSystem(boundary, boundary_model)
fluid_system = WeaklyCompressibleSPHSystem(fluid,
SummationDensity(),
state_equation,
smoothing_kernel,
smoothing_length)
semi = DummySemidiscretization()

velocities = [
[0.0; -1.0],
[1.0; 1.0],
[-1.0; 0.0],
[0.7; 0.2],
[0.3; 0.8]
]
for boundary_system in boundary_systems
system_name = boundary_system |> typeof |> nameof
density_calculator = boundary_system.boundary_model.density_calculator |>
typeof |> nameof
@testset "$system_name with $density_calculator" begin
boundary_model = boundary_system.boundary_model

@testset "Wall Velocity $v_fluid" for v_fluid in velocities
viscosity = boundary_system.boundary_model.viscosity
volume = boundary_system.boundary_model.cache.volume
viscosity = boundary_model.viscosity
pressure = boundary_model.pressure
volume = boundary_model.cache.volume

# Reset cache and perform pressure extrapolation
TrixiParticles.reset_cache!(boundary_system.boundary_model.cache,
boundary_system.boundary_model.viscosity)
TrixiParticles.reset_cache!(boundary_model.cache,
boundary_model.viscosity)
TrixiParticles.boundary_pressure_extrapolation!(Val(true),
boundary_model,
boundary_system,
fluid_system,
boundary.coordinates,
fluid.coordinates,
v_fluid,
v_boundary,
v_fluid .*
ones(size(fluid.coordinates)),
semi)

# Compute wall velocities
for particle in TrixiParticles.eachparticle(boundary_system)
if volume[particle] > eps()
pressure[particle] /= volume[particle]
TrixiParticles.compute_wall_velocity!(viscosity,
boundary_system,
boundary.coordinates,
v_boundary,
particle)
end
end
Expand All @@ -113,7 +134,7 @@
v_wall = zeros(size(boundary.coordinates))
v_wall[:, 1:particles_in_compact_support] .= -v_fluid

@test isapprox(boundary_system.boundary_model.cache.wall_velocity,
@test isapprox(boundary_model.cache.wall_velocity,
v_wall)
end

Expand All @@ -124,8 +145,9 @@
# With a staggered velocity profile, we can test the smoothed velocity field
# for a variable velocity profile.
@testset "Wall Velocity Staggered: Factor $scale" for scale in scale_factors
viscosity = boundary_system.boundary_model.viscosity
volume = boundary_system.boundary_model.cache.volume
viscosity = boundary_model.viscosity
pressure = boundary_model.pressure
volume = boundary_model.cache.volume

# Create a staggered velocity profile
v_fluid = zeros(size(fluid.coordinates))
Expand All @@ -136,23 +158,24 @@
end

# Reset cache and perform pressure extrapolation
TrixiParticles.reset_cache!(boundary_system.boundary_model.cache,
boundary_system.boundary_model.viscosity)
TrixiParticles.reset_cache!(boundary_model.cache,
boundary_model.viscosity)
TrixiParticles.boundary_pressure_extrapolation!(Val(true),
boundary_model,
boundary_system,
fluid_system,
boundary.coordinates,
fluid.coordinates,
v_fluid, v_fluid, semi)
v_boundary, v_fluid,
semi)

# Compute wall velocities
for particle in TrixiParticles.eachparticle(boundary_system)
if volume[particle] > eps()
pressure[particle] /= volume[particle]
TrixiParticles.compute_wall_velocity!(viscosity,
boundary_system,
boundary.coordinates,
particle)
v_boundary, particle)
end
end

Expand Down
Loading