fix interpolation

Author: LasNikas

src/callbacks/particle_shifting.jl

@@ -79,9 +79,6 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi,
     Wdx = smoothing_kernel(system, particle_spacing(system, 1), 1)
     h = smoothing_length(system, 1)
 
-    k_factor(::WendlandC2Kernel) = 2
-    k_factor(kernel) = 4
-
     foreach_system(semi) do neighbor_system
         u_neighbor = wrap_u(u_ode, neighbor_system, semi)
         v_neighbor = wrap_v(v_ode, neighbor_system, semi)
@@ -106,8 +103,8 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi,
 
             # Eq. 7 in Sun et al. (2017).
             # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9).
-            delta_r_ = -dt * v_max * k_factor(system_smoothing_kernel(system)) * h *
-                       (1 + R * (kernel / Wdx)^n) * m_b / (rho_a + rho_b) * grad_kernel
+            delta_r_ = -dt * v_max * 4 * h * (1 + R * (kernel / Wdx)^n) *
+                       m_b / (rho_a + rho_b) * grad_kernel
 
             # Write into the buffer
             for i in eachindex(delta_r_)

validation/vortex_street_2d/validation_vortex_street_2d.jl

@@ -3,15 +3,17 @@ using FFTW
 using CSV, DataFrames
 using Test
 
-# Results in [90k particles, 220k particles, 1.2M particles, 5M particles]
+# Results in [90k particles, 340k particles, 1.2M particles, 5M particles]
 # In the Tafuni et al. (2018), the resolution is `0.01` (5M particles).
-resolution_factor = 0.08 # [0.08, 0.05, 0.02, 0.01]
+resolution_factor = 0.08 # [0.08, 0.04, 0.02, 0.01]
 
 # ======================================================================================
 # ==== Run the simulation
 trixi_include(joinpath(validation_dir(), "vortex_street_2d", "vortex_street_2d.jl"),
               parallelization_backend=PolyesterBackend(),
-              factor_d=resolution_factor, saving_callback=nothing, tspan=(0.0, 20.0))
+              sound_speed=10.0, reynolds_number=200, saving_callback=nothing,
+              smoothing_kernel=SchoenbergQuinticSplineKernel{2}(), h_factor=1.2,
+              factor_d=resolution_factor, tspan=(0.0, 20.0))
 
 # ======================================================================================
 # ==== Read results and compute the Strouhal number

validation/vortex_street_2d/vortex_street_2d.jl

@@ -142,7 +142,7 @@ circle = SphereShape(particle_spacing, (cylinder_diameter + particle_spacing) /
                      sphere_type=RoundSphere())
 
 # Points for pressure interpolation, located at the wall interface
-const data_points =  copy(circle.coordinates)
+const data_points = copy(circle.coordinates)
 const center = SVector(cylinder_center)
 
 calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing
@@ -185,8 +185,8 @@ end
 
 # ==========================================================================================
 # ==== Simulation
-min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2)
-max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2)
+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,
@@ -203,7 +203,8 @@ info_callback = InfoCallback(interval=100)
 
 output_directory = joinpath(validation_dir(), "vortex_street_2d",
                             "out_vortex_street_dp_$(factor_d)D_c_$(sound_speed)_h_factor_$(h_factor)_" *
-                            TrixiParticles.type2string(smoothing_kernel))
+                            TrixiParticles.type2string(smoothing_kernel) *
+                            "_Re_$reynolds_number")
 
 saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory)