From 18195a54a344c82934365d8c876482349af96e4b Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 27 May 2026 15:02:51 +0200 Subject: [PATCH 1/5] Add vortex street example file Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- examples/fluid/vortex_street_2d.jl | 160 +++++++++++++++++++++++++++++ test/examples/examples_fluid.jl | 9 ++ 2 files changed, 169 insertions(+) create mode 100644 examples/fluid/vortex_street_2d.jl diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl new file mode 100644 index 0000000000..7e57041262 --- /dev/null +++ b/examples/fluid/vortex_street_2d.jl @@ -0,0 +1,160 @@ +# ========================================================================================== +# 2D Vortex Street +# +# Flow past a circular cylinder (vortex street), Tafuni et al. (2018). +# Other literature using this validation: +# Vacondio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998) +# ========================================================================================== + +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +factor_d = 0.1 # Resolution in the paper is `0.01` (5M particles) + +cylinder_diameter = 0.1 +particle_spacing = factor_d * cylinder_diameter + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 8 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 4.0) + +# Boundary geometry and initial fluid particle positions +# For better results the domain_size should be increased to +# domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter) +domain_size = (20 * cylinder_diameter, 10 * cylinder_diameter) +open_boundary_size = (particle_spacing * open_boundary_layers, domain_size[2]) + +flow_direction = SVector(1.0, 0.0) +reynolds_number = 200 +prescribed_velocity = 1.0 +fluid_density = 1000.0 +# Maximum velocity observed in the vortex street is typically around 1.5 +v_max = 1.5 +sound_speed = 10 * v_max + +pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density; + n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0) +inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density; n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +min_coords_outlet = (pipe.fluid_size[1], 0.0) +outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density; n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +n_buffer_particles = 20 * pipe.n_particles_per_dimension[2] + +cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2) +cylinder = SphereShape(particle_spacing, cylinder_diameter / 2, + cylinder_center, fluid_density; sphere_type=RoundSphere()) + +fluid = setdiff(pipe.fluid, cylinder) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 1.5 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = maximum(prescribed_velocity) * cylinder_diameter / reynolds_number + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) +viscosity = ViscosityAdami(nu=kinematic_viscosity) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +shifting_technique = ParticleShiftingTechnique(; sound_speed_factor=0.2, v_max_factor=0) + +fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=fluid_density_calculator, + state_equation, density_diffusion, viscosity, + pressure_acceleration=tensile_instability_control, + shifting_technique, buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()) + +face_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing, + reference_velocity=prescribed_velocity * flow_direction, + initial_condition=inlet.fluid, boundary_type=InFlow()) + +face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]]) +outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + initial_condition=outlet.fluid, boundary_type=OutFlow()) + +open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, + buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Boundary +wall = union(pipe.boundary, inlet.boundary, outlet.boundary) +boundary_model_wall = BoundaryModelDummyParticles(wall.density, wall.mass, + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length; + state_equation) + +boundary_system_wall = WallBoundarySystem(wall, boundary_model_wall) + +boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass, + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length; + state_equation, viscosity) + +boundary_system_cylinder = WallBoundarySystem(cylinder, boundary_model_cylinder) + +# ========================================================================================== +# ==== Simulation +boundary = union(wall, cylinder) +min_corner = minimum(boundary.coordinates .- particle_spacing, dims=2) +max_corner = maximum(boundary.coordinates .+ particle_spacing, dims=2) + +nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner), + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, open_boundary, boundary_system_wall, + boundary_system_cylinder; neighborhood_search=nhs, + parallelization_backend=PolyesterBackend()) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +extra_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # May need tuning to prevent boundary penetration + reltol=1e-4, # May need tuning to prevent boundary penetration + save_everystep=false, callback=callbacks); diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index b5e2b97199..f50a67fc79 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -525,6 +525,15 @@ @test count_rhs_allocations(sol) == 0 end + @trixi_testset "fluid/vortex_street_2d.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "vortex_street_2d.jl"), + tspan=(0.0, 0.2)) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol) == 0 + end + @trixi_testset "fluid/lid_driven_cavity_2d.jl (EDAC)" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", From 0b86be4ce7f37b55fb6ba69fad9bedf36c202f4a Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 27 May 2026 15:12:33 +0200 Subject: [PATCH 2/5] Reformat --- examples/fluid/vortex_street_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl index 7e57041262..2622d7bc78 100644 --- a/examples/fluid/vortex_street_2d.jl +++ b/examples/fluid/vortex_street_2d.jl @@ -93,7 +93,8 @@ fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_le density_calculator=fluid_density_calculator, state_equation, density_diffusion, viscosity, pressure_acceleration=tensile_instability_control, - shifting_technique, buffer_size=n_buffer_particles) + shifting_technique, + buffer_size=n_buffer_particles) # ========================================================================================== # ==== Open Boundary From f373295da76523ab299904c1eac7d1d8b0f8f1a0 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 28 May 2026 16:07:38 +0200 Subject: [PATCH 3/5] Implement copilot suggestion --- examples/fluid/vortex_street_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl index 2622d7bc78..a154f7177f 100644 --- a/examples/fluid/vortex_street_2d.jl +++ b/examples/fluid/vortex_street_2d.jl @@ -80,7 +80,7 @@ smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() -kinematic_viscosity = maximum(prescribed_velocity) * cylinder_diameter / reynolds_number +kinematic_viscosity = prescribed_velocity * cylinder_diameter / reynolds_number state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=1) From f56b8b94d02071d4caa2ba8c7d2824ec122f552f Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 2 Jun 2026 10:22:55 +0200 Subject: [PATCH 4/5] Implement suggestions --- examples/fluid/vortex_street_2d.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl index a154f7177f..912c168c77 100644 --- a/examples/fluid/vortex_street_2d.jl +++ b/examples/fluid/vortex_street_2d.jl @@ -22,15 +22,14 @@ boundary_layers = 4 # Make sure that the kernel support of fluid particles at an open boundary is always # fully sampled. # Note: Due to the dynamics at the inlets and outlets of open boundaries, -# it is recommended to use `open_boundary_layers > boundary_layers` +# it is recommended to use `open_boundary_layers > boundary_layers`. open_boundary_layers = 8 # ========================================================================================== # ==== Experiment Setup tspan = (0.0, 4.0) -# Boundary geometry and initial fluid particle positions -# For better results the domain_size should be increased to +# We use a smaller domain for a faster runtime. The original domain size in the paper is # domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter) domain_size = (20 * cylinder_diameter, 10 * cylinder_diameter) open_boundary_size = (particle_spacing * open_boundary_layers, domain_size[2]) @@ -65,7 +64,7 @@ outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_siz faces=(false, false, true, true), coordinates_eltype=Float64) -n_buffer_particles = 20 * pipe.n_particles_per_dimension[2] +n_buffer_particles = nparticles(inlet.fluid) cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2) cylinder = SphereShape(particle_spacing, cylinder_diameter / 2, @@ -104,12 +103,12 @@ face_in = ([0.0, 0.0], [0.0, domain_size[2]]) inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, open_boundary_layers, density=fluid_density, particle_spacing, reference_velocity=prescribed_velocity * flow_direction, - initial_condition=inlet.fluid, boundary_type=InFlow()) + initial_condition=inlet.fluid) face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]]) outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), open_boundary_layers, density=fluid_density, particle_spacing, - initial_condition=outlet.fluid, boundary_type=OutFlow()) + initial_condition=outlet.fluid) open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, boundary_model=open_boundary_model, @@ -119,6 +118,7 @@ open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, # ========================================================================================== # ==== Boundary wall = union(pipe.boundary, inlet.boundary, outlet.boundary) +# Free-slip boundary condition for the wall. boundary_model_wall = BoundaryModelDummyParticles(wall.density, wall.mass, AdamiPressureExtrapolation(), smoothing_kernel, smoothing_length; @@ -126,6 +126,7 @@ boundary_model_wall = BoundaryModelDummyParticles(wall.density, wall.mass, boundary_system_wall = WallBoundarySystem(wall, boundary_model_wall) +# No-slip boundary condition for the cylinder. boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass, AdamiPressureExtrapolation(), smoothing_kernel, smoothing_length; @@ -135,9 +136,8 @@ boundary_system_cylinder = WallBoundarySystem(cylinder, boundary_model_cylinder) # ========================================================================================== # ==== Simulation -boundary = union(wall, cylinder) -min_corner = minimum(boundary.coordinates .- particle_spacing, dims=2) -max_corner = maximum(boundary.coordinates .+ particle_spacing, dims=2) +min_corner = minimum(wall.coordinates .- particle_spacing, dims=2) +max_corner = maximum(wall.coordinates .+ particle_spacing, dims=2) nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner), update_strategy=ParallelUpdate()) From cd64dd2842fe222c4d602595f4118834faee3861 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 8 Jun 2026 14:24:31 +0200 Subject: [PATCH 5/5] Implement suggestions --- examples/fluid/vortex_street_2d.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl index 912c168c77..93ff39a752 100644 --- a/examples/fluid/vortex_street_2d.jl +++ b/examples/fluid/vortex_street_2d.jl @@ -1,9 +1,12 @@ # ========================================================================================== # 2D Vortex Street # -# Flow past a circular cylinder (vortex street), Tafuni et al. (2018). -# Other literature using this validation: -# Vacondio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998) +# Based on: +# A. Tafuni, J. M. Domínguez, R. Vacondio, A. J. C. Crespo. +# "A versatile algorithm for the treatment of open boundary conditions in smoothed +# particle hydrodynamics GPU models". +# Computer Methods in Applied Mechanics and Engineering, Volume 342 (2018), pages 604-624. +# https://doi.org/10.1016/j.cma.2018.08.004 # ========================================================================================== using TrixiParticles @@ -11,10 +14,10 @@ using OrdinaryDiffEqLowStorageRK # ========================================================================================== # ==== Resolution -factor_d = 0.1 # Resolution in the paper is `0.01` (5M particles) +particle_spacing_factor = 0.1 # Resolution in the paper is `0.01` (5M particles) cylinder_diameter = 0.1 -particle_spacing = factor_d * cylinder_diameter +particle_spacing = particle_spacing_factor * cylinder_diameter # Make sure that the kernel support of fluid particles at a boundary is always fully sampled boundary_layers = 4