Merge branch 'more-robust-obc' into fsi-aorta

Author: LasNikas

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -440,11 +440,19 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
             end
         end
 
-        # This only occurs if `face_normal` is not exactly normal to the `boundary_face`.
-        # In such cases, particles that are actually outside the simulation domain (outflow particles)
-        # may be incorrectly kept active as inflow particles and therefore cannot be assigned to any boundary zone.
-        # This issue should not occur and has been fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926 .
-        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
+        # Assert that every active buffer particle is assigned to a boundary zone.
+        # This should always be true if the boundary zone geometry is set up correctly.
+        # However, rare edge cases during particle conversion (`convert_particle!`)
+        # may leave a particle unassigned. Potential causes for failure:
+        # - `face_normal` is not exactly normal to the `boundary_face`
+        #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926).
+        # - Large downstream expansion can shift an inflow particle to the zone edge;
+        #   even after upstream adjustment it may remain outside
+        #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
+        # - Floating-point rounding when a particle lies almost exactly on the `boundary_face`
+        #   during transition, causing a reset just outside the zone
+        #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
+        @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle"
     end
 
     return system

src/schemes/boundary/open_boundary/system.jl

@@ -421,9 +421,22 @@ end
     # Activate a new particle in simulation domain
     transfer_particle!(fluid_system, system, particle, particle_new, v_fluid, u_fluid, v, u)
 
-    # Reset position of boundary particle
+    # Reset position of boundary particle back into the boundary zone.
+    # Particles lying almost exactly on the `boundary_face` might be reset
+    # slightly outside the boundary zone.
+    # Thus, we use a shortened vector to account for numerical precision issues.
+    reset_dist = boundary_zone.zone_width - sqrt(eps(boundary_zone.zone_width))
+    reset_vector = -boundary_zone.face_normal * reset_dist
     for dim in 1:ndims(system)
-        u[dim, particle] += boundary_zone.spanning_set[1][dim]
+        u[dim, particle] += reset_vector[dim]
+    end
+
+    # Verify the particle remains inside the boundary zone after the reset; deactivate it if not.
+    particle_coords = current_coords(u, system, particle)
+    if !is_in_boundary_zone(boundary_zone, particle_coords)
+        deactivate_particle!(system, particle, u)
+
+        return system
     end
 
     impose_rest_density!(v, system, particle, system.boundary_model)