diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 9912b91726..1d30a5274a 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -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 @@ -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 @@ -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 diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index ddd99010ad..c3e0d69ee3 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,71 @@ 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 +122,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 +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 @@ -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)) @@ -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