undo interpolation

trixi-framework · Jul 3, 2024 · e4a8c06 · e4a8c06
1 parent 1a602f7
commit e4a8c06
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 2 deletions.
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
@@ -216,7 +216,7 @@ function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi,
         throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`"))
     end
 
-    interpolate_reference_values!(system, system.fluid_system, v, u, v_ode, u_ode, semi, t)
+    # interpolate_reference_values!(system, system.fluid_system, v, u, v_ode, u_ode, semi, t)
 
     @trixi_timeit timer() "evaluate characteristics" evaluate_characteristics!(system, v, u,
                                                                                v_ode, u_ode,
@@ -574,7 +574,7 @@ function reference_values(v, system, particle, position, t)
         p_ref = initial_condition.pressure[particle]
     end
     if isnothing(v_ref)
-        v_ref = current_velocity(initial_condition.velocity, system, particle)
+        v_ref = current_velocity(v, system, particle)
     end
 
     return rho_ref, p_ref, v_ref

diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl
@@ -7,6 +7,51 @@ function interact!(dv, v_particle_system, u_particle_system,
     return dv
 end
 
+# Interaction of outlet open boundary with other systems
+function interact!(dv, v_particle_system, u_particle_system,
+                   v_neighbor_system, u_neighbor_system, neighborhood_search,
+                   particle_system::OpenBoundarySPHSystem{<:OutFlow},
+                   neighbor_system::FluidSystem)
+    (; fluid_system, sound_speed) = particle_system
+
+    system_coords = current_coordinates(u_particle_system, particle_system)
+    neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+
+    # Loop over all pairs of particles and neighbors within the kernel cutoff.
+    for_particle_neighbor(particle_system, neighbor_system,
+                          system_coords, neighbor_coords,
+                          neighborhood_search) do particle, neighbor, pos_diff, distance
+        # Only consider particles with a distance > 0.
+        distance < sqrt(eps()) && return
+
+        rho_a = particle_density(v_particle_system, particle_system, particle)
+        rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
+
+        p_a = particle_pressure(v_particle_system, particle_system, particle)
+        p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)
+
+        m_a = hydrodynamic_mass(particle_system, particle)
+        m_b = hydrodynamic_mass(neighbor_system, neighbor)
+
+        grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance)
+
+        dv_pressure = pressure_acceleration(fluid_system, fluid_system, neighbor,
+                                            m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
+                                            distance, grad_kernel, nothing)
+
+        dv_viscosity_ = dv_viscosity(particle_system, neighbor_system,
+                                     v_particle_system, v_neighbor_system,
+                                     particle, neighbor, pos_diff, distance,
+                                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+
+        for i in 1:ndims(particle_system)
+            dv[i, particle] += dv_pressure[i] + dv_viscosity_[i]
+        end
+    end
+
+    return dv
+end
+
 # For dummy particles with `ContinuityDensity`, solve the continuity equation
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system, neighborhood_search,