FSI aorta

src/general/neighborhood_search.jl

@@ -9,3 +9,29 @@ function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system,
     foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search;
                            points, parallelization_backend)
 end
+
+deactivate_out_of_bounds_particles!(system, ::Nothing, nhs, u, semi) = system
+
+function deactivate_out_of_bounds_particles!(system, ::SystemBuffer,
+                                             cell_list::FullGridCellList, u, semi)
+    (; min_corner, max_corner) = cell_list
+    (; cell_size) = get_neighborhood_search(system, semi)
+
+    # Remove the padding layer (see comment in PointNeighbors: full_grid.jl line 60)
+    min_corner_ = min_corner .+ 1001 // 1000 .* cell_size
+    max_corner_ = max_corner .- 1001 // 1000 .* cell_size
+
+    @threaded semi for particle in each_integrated_particle(system)
+        particle_position = current_coords(u, system, particle)
+
+        if !all(min_corner_ .<= particle_position .<= max_corner_)
+            deactivate_particle!(system, particle, u)
+        end
+    end
+
+    if count(system.buffer.active_particle) != system.buffer.active_particle_count[]
+        update_system_buffer!(system.buffer, semi)
+    end
+
+    return system
+end

src/io/read_vtk.jl

@@ -31,9 +31,10 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
 │ eltype: …………………………………………………………… Float64                                                          │
 │ coordinate eltype: ……………………………… Float64                                                          │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
-````
+```
 """
-function vtk2trixi(file)
+function vtk2trixi(file; element_type=Float64)
+    ELTYPE = element_type
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -42,7 +43,7 @@ function vtk2trixi(file)
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
-    coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :]
+    coordinates = convert.(ELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :])
 
     fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
     results = Dict{String, Array{Float64}}()
@@ -52,11 +53,11 @@ function vtk2trixi(file)
         all_keys = keys(point_data)
         idx = findfirst(k -> occursin(field, k), all_keys)
         if idx !== nothing
-            results[field] = ReadVTK.get_data(point_data[all_keys[idx]])
+            results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
         else
             # Use zeros as default values when a field is missing
             results[field] = field in ["mass"] ?
-                             zeros(size(coordinates, 2)) : zero(coordinates)
+                             zeros(ELTYPE, size(coordinates, 2)) : zero(coordinates)
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
@@ -65,7 +66,8 @@ function vtk2trixi(file)
                        first(results["particle_spacing"]) :
                        results["particle_spacing"]
 
-    return InitialCondition(; coordinates, particle_spacing=particle_spacing,
+    return InitialCondition(; coordinates,
+                            particle_spacing=convert(ELTYPE, particle_spacing),
                             velocity=results["velocity"],
                             mass=results["mass"],
                             density=results["density"],

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -441,10 +441,18 @@ 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 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 domain 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
 

src/schemes/boundary/open_boundary/pressure_model.jl

@@ -90,7 +90,6 @@ end
 
 function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system,
                                            boundary_zone, v, u, dt)
-    (; particle_spacing) = system.initial_condition
     (; characteristic_resistance, peripheral_resistance, compliance,
      flow_rate, pressure) = pressure_model
     (; face_normal) = boundary_zone
@@ -100,19 +99,17 @@ function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, s
                                      current_boundary_zone(system, particle),
                          each_integrated_particle(system))
 
-    # Assuming negligible transverse velocity gradients within the boundary zone,
-    # the full area of the zone is taken as the representative cross-sectional
-    # area for volumetric flow-rate estimation.
-    cross_sectional_area = length(candidates) * particle_spacing^(ndims(system) - 1)
+    # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
+    current_flow_rate = sum(candidates) do particle
+        vn = dot(current_velocity(v, system, particle), -face_normal)
+        V_particle = hydrodynamic_mass(system, particle) /
+                     current_density(v, system, particle)
+        # Cross-sectional area represented by the particle in the boundary zone
+        A_particle = V_particle / boundary_zone.zone_width
 
-    # Division inside the `sum` closure to maintain GPU compatibility
-    velocity_avg = sum(candidates) do particle
-        return dot(current_velocity(v, system, particle), -face_normal) / length(candidates)
+        return vn * A_particle
     end
 
-    # Compute volumetric flow rate: Q = v * A
-    current_flow_rate = velocity_avg * cross_sectional_area
-
     previous_pressure = pressure[]
     previous_flow_rate = flow_rate[]
     flow_rate[] = current_flow_rate

src/schemes/boundary/open_boundary/system.jl

@@ -64,6 +64,7 @@ end
 function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                             fluid_system::AbstractFluidSystem, buffer_size::Integer,
                             boundary_model,
+                            deactivate_isolated_particles=false,
                             pressure_acceleration=boundary_model isa
                                                   BoundaryModelDynamicalPressureZhang ?
                                                   fluid_system.pressure_acceleration_formulation :
@@ -82,7 +83,9 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
     mass = copy(initial_conditions.mass)
     volume = similar(initial_conditions.density)
 
-    cache = (;
+    neighbor_counter = deactivate_isolated_particles ?
+                       zeros(Int, nparticles(initial_conditions)) : nothing
+    cache = (; neighbor_counter=neighbor_counter,
              create_cache_shifting(initial_conditions, shifting_technique)...,
              create_cache_open_boundary(boundary_model, fluid_system, initial_conditions,
                                         boundary_zones_)...)
@@ -290,6 +293,12 @@ end
     return du
 end
 
+function update_positions!(system::OpenBoundarySystem, v, u, v_ode, u_ode, semi, t)
+    cell_list = get_neighborhood_search(system, semi).cell_list
+    deactivate_out_of_bounds_particles!(system, buffer(system), cell_list, u, semi)
+end
+
+
 function update_boundary_interpolation!(system::OpenBoundarySystem, v, u, v_ode, u_ode,
                                         semi, t)
     update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t)
@@ -305,6 +314,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode
         u = wrap_u(u_ode, system, semi)
         v = wrap_v(v_ode, system, semi)
 
+        calculate_neighbor_count!(system, shifting_technique(system), u, semi)
+
         @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
 
         update_pressure_model!(system, v, u, semi, integrator.dt)
@@ -418,12 +429,31 @@ end
         return system
     end
 
+    if is_isolated(system, shifting_technique(system), particle)
+        deactivate_particle!(system, particle, u)
+
+        return system
+    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 - 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)
@@ -538,16 +568,64 @@ end
         dist_free_surface = boundary_zone.zone_width - dist_to_transition
 
         if dist_free_surface < compact_support(fluid_system, fluid_system)
-            # Disable shifting for this particle.
-            # Note that Sun et al. 2017 propose a more sophisticated approach with a transition phase
-            # where only the component orthogonal to the surface normal is kept and the tangential
-            # component is set to zero. However, we assume laminar flow in the boundary zone,
-            # so we simply disable shifting completely.
+            # Ramp shifting velocity near the free surface using a kernel-weighted transition.
+            # According to our experiments, the proposed alternative approaches lead to particle disorder:
+            # - Sun et al. 2017: only use surface-tangential component
+            # - Zhang et al. 2025: disable shifting entirely
+            kernel_max = smoothing_kernel(system, 0, particle)
+            dist_from_cutoff = compact_support(fluid_system, fluid_system) -
+                               dist_free_surface
+            shifting_weight = smoothing_kernel(system, dist_from_cutoff, particle) /
+                              kernel_max
+            delta_v_ramped = delta_v(system, particle) * shifting_weight
             for dim in 1:ndims(system)
-                cache.delta_v[dim, particle] = zero(eltype(system))
+                cache.delta_v[dim, particle] = delta_v_ramped[dim]
             end
         end
     end
 
     return system
 end
+
+@inline calculate_neighbor_count!(system, ::Nothing, u, semi) = system
+
+@inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi)
+    (; neighbor_counter) = system.cache
+
+    # Skip when `deactivate_isolated_particles` is not activated
+    isnothing(neighbor_counter) && return system
+
+    set_zero!(neighbor_counter)
+
+    # Loop over all pairs of particles and neighbors within the kernel cutoff.
+    system_coords = current_coordinates(u, system)
+    foreach_point_neighbor(system, system, system_coords, system_coords, semi;
+                           points=each_integrated_particle(system)) do particle,
+                                                                       neighbor,
+                                                                       pos_diff,
+                                                                       distance
+        (particle == neighbor) && return
+        neighbor_counter[particle] += 1
+    end
+
+    return system
+end
+
+@inline is_isolated(system, ::Nothing, particle) = false
+
+@inline function is_isolated(system, ::AbstractShiftingTechnique, particle)
+    (; particle_spacing) = system.initial_condition
+    (; neighbor_counter) = system.cache
+
+    # Skip when `deactivate_isolated_particles` is not activated
+    isnothing(neighbor_counter) && return false
+
+    min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing,
+                                         compact_support(system, system)) / 10
+
+    if neighbor_counter[particle] < min_neighbors
+        return true
+    end
+
+    return false
+end

src/schemes/fluid/fluid.jl

@@ -115,6 +115,11 @@ end
 
 @inline acceleration_source(system::AbstractFluidSystem) = system.acceleration
 
+function update_positions!(system::AbstractFluidSystem, v, u, v_ode, u_ode, semi, t)
+    cell_list = get_neighborhood_search(system, semi).cell_list
+    deactivate_out_of_bounds_particles!(system, buffer(system), cell_list, u, semi)
+end
+
 function compute_density!(system, u, u_ode, semi, ::ContinuityDensity)
     # No density update with `ContinuityDensity`
     return system

test/schemes/boundary/open_boundary/pressure_model.jl

@@ -1,8 +1,7 @@
 @testset verbose=true "`RCRWindkesselModel`" begin
     @testset verbose=true "Show" begin
         pressure_model = RCRWindkesselModel(; peripheral_resistance=1.2,
-                                            compliance=0.5,
-                                            characteristic_resistance=2.3)
+                                            compliance=0.5, characteristic_resistance=2.3)
 
         show_box = """
             ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
@@ -75,7 +74,7 @@
                                          particle_spacing=0.1, face_normal=(-1.0, 0.0),
                                          density=1000.0, reference_velocity,
                                          reference_pressure=pressure_model,
-                                         open_boundary_layers=1, rest_pressure=p_0)
+                                         open_boundary_layers=17, rest_pressure=p_0)
 
             system = OpenBoundarySystem(boundary_zone; buffer_size=0,
                                         boundary_model=nothing,