From 63ebd286161d75db987869039b819dcbf49caa1d Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 1 Jul 2025 23:33:51 +0200 Subject: [PATCH 1/3] Fix Adami boundary viscosity --- .../boundary/dummy_particles/dummy_particles.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 3e2870628e..6b67beddb2 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -433,14 +433,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 @@ -451,7 +451,7 @@ function compute_adami_density!(boundary_model, system, system_coords, particle) 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 # Apply inverse state equation to compute density (not used with EDAC) @@ -597,14 +597,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 1:ndims(system) # The second term is the precalculated smoothed velocity field of the fluid. From b915ffffab5a73f66ee9fa27e663d6525e8f6996 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 22 Jul 2025 14:07:48 +0200 Subject: [PATCH 2/3] Add tests --- .../dummy_particles/dummy_particles.jl | 105 +++++++++++------- 1 file changed, 63 insertions(+), 42 deletions(-) diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index ddd99010ad..8dc4a92766 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -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 @@ -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) @@ -49,52 +50,70 @@ 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, + 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) @@ -102,9 +121,10 @@ # 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 @@ -113,7 +133,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 @@ -124,8 +144,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)) @@ -136,23 +157,23 @@ 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 From 52e114cb6b01d28246c32c5bda27ef03d8140fe5 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 22 Jul 2025 14:20:06 +0200 Subject: [PATCH 3/3] Reformat --- test/schemes/boundary/dummy_particles/dummy_particles.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 8dc4a92766..c3e0d69ee3 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -95,7 +95,8 @@ for boundary_system in boundary_systems system_name = boundary_system |> typeof |> nameof - density_calculator = boundary_system.boundary_model.density_calculator |> 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 @@ -165,7 +166,8 @@ fluid_system, boundary.coordinates, fluid.coordinates, - v_boundary, v_fluid, semi) + v_boundary, v_fluid, + semi) # Compute wall velocities for particle in TrixiParticles.eachparticle(boundary_system)