Fix tafuni extrapolation

examples/fluid/pipe_flow_2d.jl

@@ -155,7 +155,8 @@ saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
 
 extra_callback = nothing
 
-callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(),
+                        ParticleShiftingCallback(), extra_callback)
 
 sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)

src/io/write_vtk.jl

@@ -403,8 +403,10 @@ end
 function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true)
     vtk["velocity"] = [current_velocity(v, system, particle)
                        for particle in active_particles(system)]
-    vtk["density"] = current_density(v, system)
-    vtk["pressure"] = current_pressure(v, system)
+    vtk["density"] = [current_density(v, system, particle)
+                      for particle in active_particles(system)]
+    vtk["pressure"] = [current_pressure(v, system, particle)
+                       for particle in active_particles(system)]
 
     if write_meta_data
         vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters))

src/schemes/boundary/open_boundary/mirroring.jl

@@ -74,18 +74,20 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary,
             m_b = hydrodynamic_mass(fluid_system, neighbor)
             rho_b = current_density(v_fluid, fluid_system, neighbor)
             pressure_b = current_pressure(v_fluid, fluid_system, neighbor)
-            v_b = current_velocity(v_fluid, fluid_system, neighbor)
+            v_b_ = current_velocity(v_fluid, fluid_system, neighbor)
 
-            # Project `v_b` on the normal direction of the boundary zone
+            # Project `v_b_` on the normal direction of the boundary zone (only for inflow boundaries).
             # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.:
             # "Because flow from the inlet interface occurs perpendicular to the boundary,
             # only this component of interpolated velocity is kept [...]"
-            v_b = dot(v_b, boundary_zone.plane_normal) * boundary_zone.plane_normal
+            v_b = project_velocity_on_plane_normal(v_b_, boundary_zone)
 
             kernel_value = smoothing_kernel(fluid_system, distance, particle)
             grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance,
                                                 particle)
 
+            # `pos_diff` corresponds to `x_{kl} = x_k - x_l` in the paper (Tafuni et al., 2018),
+            # where `x_k` is the position of the ghost node and `x_l` is the position of the neighbor particle
             L, R = correction_arrays(kernel_value, grad_kernel, pos_diff, rho_b, m_b)
 
             correction_matrix[] += L
@@ -164,20 +166,23 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary,
 end
 
 function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{3}, rho_b, m_b)
-    x_ab = pos_diff[1]
-    y_ab = pos_diff[2]
-    z_ab = pos_diff[3]
+    # `pos_diff` corresponds to `x_{kl} = x_k - x_l` in the paper (Tafuni et al., 2018),
+    # where `x_k` is the position of the ghost node and `x_l` is the position of the neighbor particle
+    # Note that in eq. (16) and (17) the indices are swapped, i.e. `x_{lk}` is used instead of `x_{kl}`.
+    x_ba = -pos_diff[1]
+    y_ba = -pos_diff[2]
+    z_ba = -pos_diff[3]
 
     grad_W_ab_x = grad_W_ab[1]
     grad_W_ab_y = grad_W_ab[2]
     grad_W_ab_z = grad_W_ab[3]
 
     V_b = m_b / rho_b
 
-    M = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab W_ab*z_ab;
-                  grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab grad_W_ab_x*z_ab;
-                  grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab grad_W_ab_y*z_ab;
-                  grad_W_ab_z grad_W_ab_z*x_ab grad_W_ab_z*y_ab grad_W_ab_z*z_ab]
+    M = @SMatrix [W_ab W_ab*x_ba W_ab*y_ba W_ab*z_ba;
+                  grad_W_ab_x grad_W_ab_x*x_ba grad_W_ab_x*y_ba grad_W_ab_x*z_ba;
+                  grad_W_ab_y grad_W_ab_y*x_ba grad_W_ab_y*y_ba grad_W_ab_y*z_ba;
+                  grad_W_ab_z grad_W_ab_z*x_ba grad_W_ab_z*y_ba grad_W_ab_z*z_ba]
 
     L = V_b * M
 
@@ -187,17 +192,20 @@ function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{3}, rho_b, m_b)
 end
 
 function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{2}, rho_b, m_b)
-    x_ab = pos_diff[1]
-    y_ab = pos_diff[2]
+    # `pos_diff` corresponds to `x_{kl} = x_k - x_l` in the paper (Tafuni et al., 2018),
+    # where `x_k` is the position of the ghost node and `x_l` is the position of the neighbor particle
+    # Note that in eq. (16) and (17) the indices are swapped, i.e. `x_{lk}` is used instead of `x_{kl}`.
+    x_ba = -pos_diff[1]
+    y_ba = -pos_diff[2]
 
     grad_W_ab_x = grad_W_ab[1]
     grad_W_ab_y = grad_W_ab[2]
 
     V_b = m_b / rho_b
 
-    M = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab;
-                  grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab;
-                  grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab]
+    M = @SMatrix [W_ab W_ab*x_ba W_ab*y_ba;
+                  grad_W_ab_x grad_W_ab_x*x_ba grad_W_ab_x*y_ba;
+                  grad_W_ab_y grad_W_ab_y*x_ba grad_W_ab_y*y_ba]
     L = V_b * M
 
     R = V_b * SVector(W_ab, grad_W_ab_x, grad_W_ab_y)
@@ -211,3 +219,13 @@ function mirror_position(particle_coords, boundary_zone)
 
     return particle_coords - 2 * dist * boundary_zone.plane_normal
 end
+
+project_velocity_on_plane_normal(vel, boundary_zone) = vel
+
+function project_velocity_on_plane_normal(vel, boundary_zone::BoundaryZone{InFlow})
+    # Project `v_b` on the normal direction of the boundary zone
+    # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.:
+    # "Because flow from the inlet interface occurs perpendicular to the boundary,
+    # only this component of interpolated velocity is kept [...]"
+    return dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal
+end

test/examples/examples_fluid.jl

@@ -339,7 +339,8 @@
     end
 
     @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin
-        steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=10)
+        steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=10,
+                                                          reltol=1e-3)
 
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(examples_dir(), "fluid",
@@ -353,13 +354,12 @@
     end
 
     @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin
-        steady_state_reached = SteadyStateReachedCallback(; interval=1,
-                                                          interval_size=10,
-                                                          abstol=1.0e-5, reltol=1.0e-4)
+        steady_state_reached = SteadyStateReachedCallback(; interval=1, interval_size=10,
+                                                          reltol=1e-3)
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(examples_dir(), "fluid",
                                                   "pipe_flow_2d.jl"),
-                                         extra_callback=steady_state_reached,
+                                         extra_callback=steady_state_reached, dtmax=2e-3,
                                          tspan=(0.0, 1.5), viscosity_boundary=nothing)
 
         # Make sure that the simulation is terminated after a reasonable amount of time