Validate open boundaries

examples/fluid/vortex_street_2d.jl

@@ -0,0 +1,158 @@
+# ==========================================================================================
+# 2D Vortex Street
+#
+# Flow past a circular cylinder (vortex street), Tafuni et al. (2018).
+# Other literature using this validation:
+# Vacandio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998)
+# ==========================================================================================
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== 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)
+
+flow_direction = [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
+
+boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+                 domain_size[2])
+
+pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
+                       n_layers=boundary_layers, velocity=[prescribed_velocity, 0.0],
+                       faces=(false, false, true, true))
+
+# Shift pipe walls in negative x-direction for the inflow
+pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
+
+n_buffer_particles = 20 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1)
+
+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 = 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 = TransportVelocityAdami(background_pressure=5 * fluid_density *
+#                                                                 sound_speed^2)
+
+shifting_technique = ParticleShiftingTechnique(; sound_speed_factor=0.2, v_max_factor=0)
+
+fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           density_diffusion=density_diffusion,
+                                           smoothing_length, viscosity=viscosity,
+                                           pressure_acceleration=tensile_instability_control,
+                                           shifting_technique=shifting_technique,
+                                           buffer_size=n_buffer_particles)
+
+# ==========================================================================================
+# ==== Open Boundary
+open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring())
+# open_boundary_model = BoundaryModelDynamicalPressureZhang()
+
+inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, domain_size[2]]),
+                      face_normal=flow_direction, open_boundary_layers,
+                      reference_velocity=(pos, t) -> SVector(prescribed_velocity, 0.0),
+                      density=fluid_density, particle_spacing)
+
+outflow = BoundaryZone(;
+                       boundary_face=([pipe.fluid_size[1], 0.0],
+                                      [pipe.fluid_size[1], pipe.fluid_size[2]]),
+                       face_normal=(-flow_direction), particle_spacing,
+                       density=fluid_density, open_boundary_layers)
+
+# While the outlet velocity is not explicitly prescribed,
+# initializing it in the x-direction helps to suppress initial pressure waves.
+outflow.initial_condition.velocity[1, :] .= prescribed_velocity
+
+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
+boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             state_equation=state_equation,
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system_wall = WallBoundarySystem(pipe.boundary, boundary_model)
+
+boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass,
+                                                      AdamiPressureExtrapolation(),
+                                                      state_equation=state_equation,
+                                                      viscosity=viscosity,
+                                                      smoothing_kernel, smoothing_length)
+
+boundary_system_cylinder = WallBoundarySystem(cylinder, boundary_model_cylinder)
+
+# ==========================================================================================
+# ==== Simulation
+min_corner = minimum(pipe.boundary.coordinates .- 5 * particle_spacing, dims=2)
+max_corner = maximum(pipe.boundary.coordinates .+ 5 * particle_spacing, dims=2)
+cell_list = FullGridCellList(; min_corner, max_corner)
+
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list,
+                                                update_strategy=ParallelUpdate())
+
+semi = Semidiscretization(fluid_system, open_boundary, boundary_system_wall,
+                          boundary_system_cylinder; neighborhood_search=neighborhood_search,
+                          parallelization_backend=PolyesterBackend())
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+
+saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory="out")
+
+extra_callback = nothing
+
+callbacks = CallbackSet(info_callback, UpdateCallback(), saving_callback, extra_callback)
+
+sol = solve(ode, RDPK3SpFSAL35(),
+            abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

test/Project.toml

@@ -2,6 +2,7 @@
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"

test/examples/examples_fluid.jl

@@ -486,6 +486,15 @@
         @test count_rhs_allocations(sol, semi) == 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, semi) == 0
+    end
+
     @trixi_testset "fluid/lid_driven_cavity_2d.jl (EDAC)" begin
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(examples_dir(), "fluid",

test/validation/validation.jl

@@ -112,4 +112,32 @@
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
     end
+
+    @trixi_testset "vortex_street_2d" begin
+        @trixi_testset "Validation" begin
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(validation_dir(),
+                                                      "vortex_street_2d",
+                                                      "validation_vortex_street_2d.jl"),
+                                             output_directory=joinpath(validation_dir(),
+                                                                       "vortex_street_2d",
+                                                                       "out_test"),
+                                             tspan=(0.0, 1.0))
+            @test sol.retcode == ReturnCode.Success
+            @test count_rhs_allocations(sol, semi) == 0
+        end
+
+        @trixi_testset "Plot" begin
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(validation_dir(),
+                                                      "vortex_street_2d",
+                                                      "plot_vortex_street_2d.jl")) [
+                # exact info outputs
+                r"┌ Info: Strouhal number\n└\s+round\(strouhal_number, digits = 3\) = 0\.228\n",
+                r"┌ Info: Fraction of the dominant frequency band in the total spectrum\n└\s+round\(integral_peak / integral_total, digits = 3\) = 0\.39\n",
+                r"┌ Info: C_L_max for the periodic shedding\n└\s+round\(maximum\(f_lift_cut\), digits = 3\) = 0\.7\n",
+                r"┌ Info: C_D_max for the periodic shedding\n└\s+round\(maximum\(f_drag\[start_index:end\]\), digits = 3\) = 1\.631\n"
+            ]
+        end
+    end
 end

validation/vortex_street_2d/plot_vortex_street_2d.jl

@@ -0,0 +1,82 @@
+using TrixiParticles
+using CSV, DataFrames, Plots
+using FFTW
+
+# Results in 1.2M particles.
+# In the Tafuni et al. (2018), the resolution is `0.01` (5M particles).
+resolution_factor = 0.02
+cylinder_diameter = 0.1
+prescribed_velocity = 1.0
+
+reynolds_number = 200
+
+output_directory = joinpath(validation_dir(), "vortex_street_2d")
+
+# ======================================================================================
+# ==== Read results
+data = CSV.read(joinpath(output_directory, "resulting_force.csv"), DataFrame)
+
+step = 1
+times = data[!, "time"][1:step:end]
+
+f_lift = data[!, "f_l_fluid_1"][1:step:end]
+f_drag = data[!, "f_d_fluid_1"][1:step:end]
+
+# ======================================================================================
+# ==== Compute strouhal number
+t_start = 6.0
+start_index = findfirst(t -> t ≥ t_start, times)
+times_cut = times[start_index:end]
+dt = times_cut[2] - times_cut[1]
+
+f_lift_cut = f_lift[start_index:end]
+
+# Compute the frequency for the FFT.
+# For N time samples with uniform time steps dt, the corresponding frequencies are:
+# f_k = k / (N * dt), where k = 0, 1, ..., N-1.
+# This gives the frequency bins in Hz, matching the order of FFT.
+N = length(f_lift_cut)
+frequencies = (0:(N - 1)) / (N * dt)
+
+spectrum = abs.(fft(f_lift_cut))
+spectrum_half = spectrum[1:div(N, 2)]
+
+# For real-valued signals, the FFT output is symmetric.
+# Only the first half (up to the Nyquist frequency) contains unique, physically meaningful frequency components.
+# We therefore analyze only the first N/2 values of the frequency spectrum.
+f_dominant = frequencies[argmax(spectrum_half)]
+
+# Verify whether the dominant frequency is indeed unique.
+# In theory, for a purely harmonic oscillation, the spectrum should exhibit only a single dominant frequency component.
+delta = 2 * (frequencies[2] - frequencies[1])
+frequency_band = (abs.(frequencies[1:div(N, 2)] .- f_dominant) .< delta)
+
+strouhal_number = f_dominant * cylinder_diameter / prescribed_velocity
+
+# ======================================================================================
+# ==== Validate the frequency spectrum
+integral_total = sum(spectrum_half)
+integral_peak = sum(spectrum_half[frequency_band])
+
+@info "Strouhal number" round(strouhal_number, digits=3)
+@info "Fraction of the dominant frequency band in the total spectrum" round(integral_peak /
+                                                                            integral_total,
+                                                                            digits=3)
+@info "C_L_max for the periodic shedding" round(maximum(f_lift_cut), digits=3)
+@info "C_D_max for the periodic shedding" round(maximum(f_drag[start_index:end]), digits=3)
+
+dp = round(Int, 1 / resolution_factor)
+plot_title = "Drag and lift force coefficients (Δx = d/$(dp))"
+pC = plot(times, f_lift, ylims=(-1, 3), xlims=(0, 20), label="C_L", color=:red, linewidth=2)
+plot!(pC, times, f_drag, ylims=(-1, 3), xlims=(0, 20), label="C_D", color=:blue,
+      linewidth=2, title=plot_title, xlabel="t (s)")
+plot!(pC, top_margin=2Plots.mm)
+
+pS = plot(frequencies[1:div(N, 2)], spectrum_half, xlabel="Frequency (Hz)", size=(400, 200),
+          ylabel="Amplitude",
+          title="Frequency Spectrum (St = $(round(strouhal_number, digits=4)))",
+          label=nothing, linewidth=2)
+plot!(pC, top_margin=2Plots.mm)
+
+p = plot(pC, pS, layout=@layout([a; b{0.3h}]), size=(800, 800))
+plot!(p, right_margin=5Plots.mm)