diff --git a/NEWS.md b/NEWS.md index fe05cfcf50..f891c9c374 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,13 +9,10 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### API Changes - `DensityDiffusionAntuono` now only takes only one kwarg `delta` (#1142). +- Clipping of negative pressure values in the `DummyParticleBoundaryModel` is now disabled + by default and can be enabled with the keyword argument `clip_negative_pressure=true` (#1143). - Return type of `vtk2trixi` changed to `NamedTuple` including an optional - `:initial_condition` field if `create_initial_condition=true` is passed. (#959) - -## Version 0.4.4 - -### API Changes - + `:initial_condition` field if `create_initial_condition=true` is passed (#959). - Custom quantities called in the `PostprocessCallback` are now passed CPU arrays when the simulation is run on a GPU (#1065). diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 1f36fac5f1..def3729c31 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -79,13 +79,16 @@ boundary_density_calculator = AdamiPressureExtrapolation() viscosity_wall = nothing # For a no-slip boundary condition, define a wall viscosity: # viscosity_wall = viscosity_fluid + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, smoothing_kernel, smoothing_length, correction=nothing, reference_particle_spacing=0, - viscosity=viscosity_wall) + viscosity=viscosity_wall, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl index a073c001ec..8f521696f3 100644 --- a/examples/fluid/dam_break_3d.jl +++ b/examples/fluid/dam_break_3d.jl @@ -56,10 +56,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl index f712f88d93..a185556673 100644 --- a/examples/fluid/falling_water_column_2d.jl +++ b/examples/fluid/falling_water_column_2d.jl @@ -53,10 +53,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 1e952316d7..79a0e86246 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -77,12 +77,15 @@ sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() wall_viscosity = nu + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, fluid_smoothing_kernel, fluid_smoothing_length, viscosity=ViscosityAdami(nu=wall_viscosity), - reference_particle_spacing=fluid_particle_spacing) + reference_particle_spacing=fluid_particle_spacing, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model, adhesion_coefficient=1.0) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 4667b1af42..bc0e48e575 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -118,15 +118,19 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model_tank = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) boundary_model_gate = BoundaryModelDummyParticles(gate.density, gate.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) boundary_system_tank = WallBoundarySystem(tank.boundary, boundary_model_tank) boundary_system_gate = WallBoundarySystem(gate, boundary_model_gate, @@ -145,7 +149,8 @@ boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites, hydrodynamic_masses, state_equation=state_equation, AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index 44d86c5317..dfc0b8a78a 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -88,10 +88,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) @@ -121,7 +124,8 @@ boundary_model_structure = BoundaryModelMonaghanKajtar(k_structure, spacing_rati # hydrodynamic_masses, # state_equation=state_equation, # boundary_density_calculator, -# smoothing_kernel, smoothing_length) +# smoothing_kernel, smoothing_length, +# clip_negative_pressure=true) structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index 026e35793e..ce9226285d 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -67,10 +67,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length) + fluid_smoothing_kernel, fluid_smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/falling_rotating_rigid_squares_2d.jl b/examples/fsi/falling_rotating_rigid_squares_2d.jl index 5eb1160ea0..a5cfc2d7a2 100644 --- a/examples/fsi/falling_rotating_rigid_squares_2d.jl +++ b/examples/fsi/falling_rotating_rigid_squares_2d.jl @@ -83,10 +83,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length) + fluid_smoothing_kernel, fluid_smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index e938cbf6f8..cdd6eb822b 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -71,10 +71,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = BernoulliPressureExtrapolation() + +# Clip negative boundary pressure values to avoid sticking artifacts at the boundary. boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length) + fluid_smoothing_kernel, fluid_smoothing_length, + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) @@ -93,7 +96,8 @@ structure_boundary_model_1 = BoundaryModelDummyParticles(hydrodynamic_densites_1 state_equation=state_equation, boundary_density_calculator, fluid_smoothing_kernel, - fluid_smoothing_length) + fluid_smoothing_length, + clip_negative_pressure=true) hydrodynamic_densites_2 = fluid_density * ones(size(sphere2.density)) hydrodynamic_masses_2 = hydrodynamic_densites_2 * @@ -104,7 +108,8 @@ structure_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2 state_equation=state_equation, boundary_density_calculator, fluid_smoothing_kernel, - fluid_smoothing_length) + fluid_smoothing_length, + clip_negative_pressure=true) structure_system_1 = TotalLagrangianSPHSystem(sphere1, structure_smoothing_kernel, diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index 6927e0d7e5..4a22db8e4f 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -3,6 +3,7 @@ density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, state_equation=nothing, correction=nothing, + clip_negative_pressure=false, reference_particle_spacing=0.0) Boundary model for [`WallBoundarySystem`](@ref). @@ -10,9 +11,9 @@ Boundary model for [`WallBoundarySystem`](@ref). # Arguments - `initial_density`: Vector holding the initial density of each boundary particle. - `hydrodynamic_mass`: Vector holding the "hydrodynamic mass" of each boundary particle. - See description above for more information. + See [the docs](@ref boundary_models) for more information. - `density_calculator`: Strategy to compute the hydrodynamic density of the boundary particles. - See description below for more information. + See [the docs](@ref boundary_models) for more information. - `smoothing_kernel`: Smoothing kernel should be the same as for the adjacent fluid system. - `smoothing_length`: Smoothing length should be the same as for the adjacent fluid system. @@ -22,6 +23,15 @@ Boundary model for [`WallBoundarySystem`](@ref). - `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). - `viscosity`: Slip (default) or no-slip condition. See description below for further information. +- `clip_negative_pressure=false`: Clip negative boundary pressures to avoid sticking + artifacts from attractive fluid-boundary forces at free + surfaces. Note that this is not a correct formulation + at interior boundaries (away from free surfaces). + Disable this option when simulating closed systems without + free surfaces to avoid artificially increased boundary + pressures that cause larger gaps between fluid and boundary + in areas of low pressure, against which the particle + shifting technique is fighting. - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, which currently is only needed when using surface tension. # Examples @@ -39,7 +49,7 @@ boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExt BoundaryModelDummyParticles(AdamiPressureExtrapolation, ViscosityAdami) ``` """ -struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C} +struct BoundaryModelDummyParticles{DC, SE, CLIP, ELTYPE <: Real, VECTOR, K, V, COR, C} pressure :: VECTOR # Vector{ELTYPE} hydrodynamic_mass :: VECTOR # Vector{ELTYPE} state_equation :: SE @@ -49,14 +59,28 @@ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C} viscosity :: V correction :: COR cache :: C + # Store this both as field and type parameter to avoid annoying hand-written + # `adapt_structure` functions. + clip_negative_pressure::Bool + + function BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation, + density_calculator, smoothing_kernel, + smoothing_length, viscosity, correction, + cache, clip_negative_pressure) + return new{typeof(density_calculator), typeof(state_equation), + clip_negative_pressure, eltype(pressure), typeof(pressure), + typeof(smoothing_kernel), typeof(viscosity), typeof(correction), + typeof(cache)}(pressure, hydrodynamic_mass, state_equation, + density_calculator, smoothing_kernel, smoothing_length, + viscosity, correction, cache, clip_negative_pressure) + end end -# The default constructor needs to be accessible for Adapt.jl to work with this struct. -# See the comments in general/gpu.jl for more details. function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, state_equation=nothing, correction=nothing, + clip_negative_pressure=false, reference_particle_spacing=0.0) pressure = initial_boundary_pressure(initial_density, density_calculator, state_equation) @@ -82,10 +106,17 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation, density_calculator, smoothing_kernel, - smoothing_length, viscosity, correction, cache) + smoothing_length, viscosity, correction, cache, + clip_negative_pressure) +end + +@inline function Base.ndims(boundary_model::BoundaryModelDummyParticles) + return ndims(boundary_model.smoothing_kernel) end -@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel) +@inline function clip_negative_pressure(::BoundaryModelDummyParticles{<:Any, <:Any, CLIP}) where {CLIP} + return CLIP +end @doc raw""" AdamiPressureExtrapolation(; pressure_offset=0, allow_loop_flipping=true) @@ -108,7 +139,6 @@ end loop is not thread parallelizable. This can cause error variations between simulations with different numbers of threads. - """ struct AdamiPressureExtrapolation{ELTYPE} pressure_offset :: ELTYPE @@ -423,33 +453,33 @@ end function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityDensity}, system, v, u, v_ode, u_ode, semi) - - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). @threaded semi for particle in eachparticle(system) - apply_state_equation!(boundary_model, current_density(v, system, particle), - particle) + @inbounds apply_state_equation!(boundary_model, + current_density(v, system, particle), particle) end return boundary_model 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. -@inline function apply_state_equation!(boundary_model, density, particle) - boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0) +@propagate_inbounds function apply_state_equation!(boundary_model, density, particle) + pressure = boundary_model.state_equation(density) + + # This is determined statically and has therefore no overhead. + if clip_negative_pressure(boundary_model) + # Clip pressure to avoid sticking artifacts from negative pressures at free surfaces. + pressure = max(pressure, 0) + end + + boundary_model.pressure[particle] = pressure end @propagate_inbounds function apply_state_equation!(boundary_model::BoundaryModelDummyParticles{<:Any, - <:Any, - <:Any, Nothing}, density, particle) # Contact-only wall setups can reuse dummy particles for rigid contact without # configuring a hydrodynamic state equation. In that case, keep the auxiliary wall # pressure at zero because it is only meaningful for fluid-coupled updates. - boundary_model.pressure[particle] = zero(eltype(boundary_model.pressure)) + boundary_model.pressure[particle] = 0 end function compute_pressure!(boundary_model, @@ -518,9 +548,11 @@ function compute_adami_density!(boundary_model, system, v, particle) compute_wall_velocity!(viscosity, system, v, particle) end - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). - @inbounds pressure[particle] = max(pressure[particle], 0) + # Clip pressure to avoid sticking artifacts from negative pressures at free surfaces. + # This is determined statically and has therefore no overhead. + if clip_negative_pressure(boundary_model) + @inbounds pressure[particle] = max(pressure[particle], 0) + end # Apply inverse state equation to compute density (not used with EDAC) inverse_state_equation!(density, state_equation, pressure, particle) diff --git a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl index b85a0d5e52..bfa7affd2b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl +++ b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl @@ -25,7 +25,14 @@ The speed of sound is initialized as `min_sound_speed`. - `max_sound_speed=100.0f0`: The maximum permissible speed of sound. - `exponent`: An exponent, typically 7 for water simulations. - `background_pressure=0.0f0`: A constant background pressure. -- `clip_negative_pressure=false`: When true, negative pressure values are clipped to 0.0. This can prevent spurious surface tension effects but might allow for unphysical fluid rarefaction. +- `clip_negative_pressure=false`: When true, negative pressure values are clipped to 0.0. + This can prevent spurious surface tension effects but allows for unphysical fluid + rarefaction and causes incorrect fluid behavior in regions of low pressure. + Note that negative pressure values of boundary particles in the + [`BoundaryModelDummyParticles`](@ref) can lead to sticking artifacts at the boundary. + The `BoundaryModelDummyParticles` itself also provides an option to clip negative + pressure of boundary particles, which is a more targeted way to prevent sticking artifacts + without affecting the fluid behavior in the rest of the domain. """ struct StateEquationAdaptiveCole{ELTYPE, CLIP, SR} # Boolean to clip negative pressure sound_speed_ref :: SR diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 125f3bd38e..bbb5cfcac6 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -39,16 +39,20 @@ viscosity_fluid=ViscosityMorris(nu=0.0015), fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), - "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3,), - "WCSPH with SchoenbergQuarticSplineKernel" => (smoothing_length=1.1, + "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3 * + fluid_particle_spacing,), + "WCSPH with SchoenbergQuarticSplineKernel" => (smoothing_length=1.1 * + fluid_particle_spacing, smoothing_kernel=SchoenbergQuarticSplineKernel{2}()), - "WCSPH with SchoenbergQuinticSplineKernel" => (smoothing_length=1.1, + "WCSPH with SchoenbergQuinticSplineKernel" => (smoothing_length=1.1 * + fluid_particle_spacing, smoothing_kernel=SchoenbergQuinticSplineKernel{2}()), - "WCSPH with WendlandC2Kernel" => (smoothing_length=1.5, + "WCSPH with WendlandC2Kernel" => (smoothing_length=1.5 * fluid_particle_spacing, smoothing_kernel=WendlandC2Kernel{2}()), - "WCSPH with WendlandC4Kernel" => (smoothing_length=1.75, + "WCSPH with WendlandC4Kernel" => (smoothing_length=1.75 * + fluid_particle_spacing, smoothing_kernel=WendlandC4Kernel{2}()), - "WCSPH with WendlandC6Kernel" => (smoothing_length=2.0, + "WCSPH with WendlandC6Kernel" => (smoothing_length=2.0 * fluid_particle_spacing, smoothing_kernel=WendlandC6Kernel{2}()), "EDAC with source term damping" => (source_terms=SourceTermDamping(damping_coefficient=1e-4), fluid_system=EntropicallyDampedSPHSystem(tank.fluid, @@ -426,6 +430,8 @@ @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=false, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelMirroringTafuni(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success @@ -433,9 +439,11 @@ end @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (WCSPH)" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=true, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelMirroringTafuni(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 3117a67ca5..afdf95b0e3 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -508,14 +508,12 @@ end joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + wcsph=false, coordinates_eltype=Float32, - open_boundary_model=BoundaryModelMirroringTafuni(), + open_boundary_model=BoundaryModelMirroringTafuni(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_velocity_out=nothing, parallelization_backend=Main.parallelization_backend) [ r"\[ Info: To move data to the GPU, `semidiscretize` creates a deep copy.*\n" ] @@ -530,17 +528,12 @@ end joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - wcsph=true, sound_speed=20.0f0, + wcsph=true, coordinates_eltype=Float32, open_boundary_model=BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_pressure_out=nothing, - reference_velocity_out=nothing, parallelization_backend=Main.parallelization_backend) [ r"\[ Info: To move data to the GPU, `semidiscretize` creates a deep copy.*\n" ] diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 239f992e67..3c096f83b7 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -10,6 +10,63 @@ @test repr(boundary_model) == expected_repr end + @testset "Pressure clipping" begin + state_equation = StateEquationCole(sound_speed=10.0, + reference_density=1000.0, + exponent=7, + clip_negative_pressure=false) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.1 + + boundary_model_clip = BoundaryModelDummyParticles([1000.0], [1.0], + SummationDensity(), + smoothing_kernel, + smoothing_length, + state_equation=state_equation, + clip_negative_pressure=true) + boundary_model_no_clip = BoundaryModelDummyParticles([1000.0], [1.0], + SummationDensity(), + smoothing_kernel, + smoothing_length, + state_equation=state_equation) + + @test TrixiParticles.clip_negative_pressure(boundary_model_clip) + @test !TrixiParticles.clip_negative_pressure(boundary_model_no_clip) + + TrixiParticles.apply_state_equation!(boundary_model_clip, 900.0, 1) + TrixiParticles.apply_state_equation!(boundary_model_no_clip, 900.0, 1) + + @test boundary_model_clip.pressure[1] == 0 + @test boundary_model_no_clip.pressure[1] < 0 + + # Test that clipping behavior is preserved through `adapt_structure`. + adapted_model_clip = TrixiParticles.Adapt.adapt(Array, boundary_model_clip) + adapted_model_no_clip = TrixiParticles.Adapt.adapt(Array, boundary_model_no_clip) + @test TrixiParticles.clip_negative_pressure(adapted_model_clip) + @test !TrixiParticles.clip_negative_pressure(adapted_model_no_clip) + + boundary_model_adami_clip = BoundaryModelDummyParticles([1000.0], [1.0], + AdamiPressureExtrapolation(), + smoothing_kernel, + smoothing_length, + state_equation=state_equation, + clip_negative_pressure=true) + boundary_model_adami_no_clip = BoundaryModelDummyParticles([1000.0], [1.0], + AdamiPressureExtrapolation(), + smoothing_kernel, + smoothing_length, + state_equation=state_equation) + + for model in (boundary_model_adami_clip, boundary_model_adami_no_clip) + model.cache.volume[1] = 1 + model.pressure[1] = -1 + TrixiParticles.compute_adami_density!(model, nothing, nothing, 1) + end + + @test boundary_model_adami_clip.pressure[1] == 0 + @test boundary_model_adami_no_clip.pressure[1] == -1 + end + @testset verbose=true "Viscosity Adami/Bernoulli: Wall Velocity" begin particle_spacing = 0.1