Fix open boundary shifting (#1195)

* Fix open boundary shifting

* Add docs and test

* Update NEWS.md

* Reformat

* Implement suggestions

* Implement suggestions

* Fix

* Add check for boundary zone width

Author: efaulhaber

NEWS.md

@@ -19,6 +19,10 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - Fixed interpolation with `cut_off_bnd=false` and `interpolate_wall_velocity=true` (#1185).
 - Fixed interpolation from a `PostprocessCallback` on GPUs (#1192).
 - Fixed interpolation with different smoothing length on the GPU (#952).
+- Disabled shifting in the free-surface region of `OpenBoundarySystem` with
+  `BoundaryModelDynamicalPressureZhang` and ramp it only between one and two compact
+  supports away from the free surface (#1195).
+
 ## Version 0.5
 
 ### API Changes

docs/src/systems/boundary.md

@@ -344,6 +344,14 @@ To further handle incomplete kernel support, for example in the viscous term of
 the updated velocity of particles within the [`BoundaryZone`](@ref) is projected onto the face normal,
 so that only the component in flow direction is kept.
 
+When a shifting technique or transport velocity formulation is used with this boundary model,
+the shifting velocity is modified near the free surface of the [`BoundaryZone`](@ref).
+The free-surface region is defined as the layer within one compact support of the free surface;
+inside this region, shifting is set to zero because the kernel support is incomplete.
+From one to two compact supports away from the free surface, the shifting velocity is ramped
+from zero to the unmodified value using a kernel-weighted transition.
+Particles farther away from the free surface use the unmodified shifting velocity.
+
 # Pressure Models
 ```@autodocs
 Modules = [TrixiParticles]

src/schemes/boundary/open_boundary/system.jl

@@ -81,6 +81,16 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                                               density_diffusion(fluid_system) : nothing)
     boundary_zones_ = filter(bz -> !isnothing(bz), boundary_zones)
 
+    if boundary_model isa BoundaryModelDynamicalPressureZhang &&
+       !isnothing(shifting_technique)
+        # When dynamical pressure is used with shifting, the shifting is ramped up until
+        # it reaches the full value at two compact supports away from the free surface
+        # of the boundary zone.
+        # In order to prevent discontinuities in the shifting velocity, the boundary zone
+        # must be wide enough to accommodate this ramping region.
+        check_boundary_zone_widths(boundary_zones_, fluid_system)
+    end
+
     initial_conditions = union((bz.initial_condition for bz in boundary_zones_)...)
 
     buffer = SystemBuffer(nparticles(initial_conditions), buffer_size)
@@ -605,15 +615,30 @@ end
                                  -boundary_zone.face_normal)
         dist_free_surface = boundary_zone.zone_width - dist_to_transition
 
-        if dist_free_surface < compact_support(fluid_system, fluid_system)
-            # Ramp shifting velocity near the free surface using a kernel-weighted transition.
+        # Width of the free-surface region where shifting is disabled.
+        free_surface_region_width = compact_support(fluid_system, fluid_system)
+        # Width of the ramp region behind the free-surface region where shifting is ramped
+        # from zero to the unmodified value.
+        ramp_width = free_surface_region_width
+
+        if dist_free_surface <= free_surface_region_width
+            # The free-surface region is the layer within one compact support of the
+            # free surface. Shifting is disabled there to avoid shifting particles with
+            # incomplete support.
+            for dim in 1:ndims(system)
+                cache.delta_v[dim, particle] = zero(eltype(system))
+            end
+        elseif dist_free_surface < free_surface_region_width + ramp_width
+            # Between one and two compact supports from the free surface, shifting is
+            # ramped from zero to the unmodified shifting velocity with 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) /
+            # Distance from the area where full shifting is applied.
+            dist_from_ramp_end = free_surface_region_width + ramp_width - dist_free_surface
+            shifting_weight = smoothing_kernel(system, dist_from_ramp_end, particle) /
                               kernel_max
             delta_v_ramped = delta_v(system, particle) * shifting_weight
             for dim in 1:ndims(system)
@@ -689,8 +714,7 @@ end
 
 @inline use_open_boundary_interpolation_neighbor(system) = false
 
-function check_configuration(system::OpenBoundarySystem, systems,
-                             neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch)
+function check_configuration(system::OpenBoundarySystem, systems, neighborhood_search)
     (; boundary_model, boundary_zones) = system
 
     # Store index of the fluid system. This is necessary for re-linking
@@ -710,3 +734,19 @@ function check_configuration(system::OpenBoundarySystem, systems,
                             "See the PointNeighbors.jl documentation for more details."))
     end
 end
+
+function check_boundary_zone_widths(boundary_zones, fluid_system)
+    support = compact_support(system_smoothing_kernel(fluid_system),
+                              initial_smoothing_length(fluid_system))
+    minimum_width = 2 * support
+
+    for (i, zone) in enumerate(boundary_zones)
+        if zone.zone_width < minimum_width
+            throw(ArgumentError("boundary zone $i has width $(zone.zone_width), " *
+                                "but must be at least two compact supports " *
+                                "($(minimum_width))"))
+        end
+    end
+
+    return nothing
+end

test/schemes/boundary/open_boundary/open_boundary.jl

@@ -3,6 +3,56 @@ include("mirroring.jl")
 include("boundary_zone.jl")
 include("pressure_model.jl")
 
+@testset verbose=true "Free-Surface Shifting Ramp" begin
+    particle_spacing = 0.1
+    smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+    smoothing_length = 0.5
+    compact_support = TrixiParticles.compact_support(smoothing_kernel, smoothing_length)
+
+    zone_width = 3 * compact_support
+    distances_from_free_surface = compact_support .* [0.5, 1.0, 1.5, 2.0]
+    coordinates = stack(SVector(-zone_width + distance, 0.5)
+                        for distance in distances_from_free_surface)
+    initial_condition = InitialCondition(; coordinates, density=1.0,
+                                         particle_spacing)
+
+    fluid = RectangularShape(particle_spacing, (2, 2), (0.0, 0.0), density=1.0)
+    fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                               smoothing_length, sound_speed=1.0,
+                                               shifting_technique=ParticleShiftingTechnique())
+
+    boundary_zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
+                                 face_normal=[1.0, 0.0], boundary_type=InFlow(),
+                                 open_boundary_layers=round(Int,
+                                                            zone_width / particle_spacing),
+                                 initial_condition, density=1.0, particle_spacing)
+    open_boundary = OpenBoundarySystem(boundary_zone; fluid_system, buffer_size=0,
+                                       boundary_model=BoundaryModelDynamicalPressureZhang())
+    open_boundary.boundary_zone_indices .= 1
+
+    unmodified_delta_v = [1.0, -2.0]
+    for particle in axes(open_boundary.cache.delta_v, 2)
+        open_boundary.cache.delta_v[:, particle] .= unmodified_delta_v
+    end
+
+    TrixiParticles.modify_shifting_at_free_surfaces!(open_boundary,
+                                                     initial_condition.coordinates,
+                                                     DummySemidiscretization())
+
+    kernel_max = TrixiParticles.smoothing_kernel(open_boundary, 0, 1)
+
+    @test iszero(open_boundary.cache.delta_v[:, 1])
+    @test iszero(open_boundary.cache.delta_v[:, 2])
+
+    dist_from_cutoff = 2 * compact_support - distances_from_free_surface[3]
+    shifting_weight = TrixiParticles.smoothing_kernel(open_boundary, dist_from_cutoff,
+                                                      3) / kernel_max
+    @test isapprox(open_boundary.cache.delta_v[:, 3],
+                   unmodified_delta_v * shifting_weight)
+
+    @test open_boundary.cache.delta_v[:, 4] == unmodified_delta_v
+end
+
 @testset verbose=true "Calculate Flow Rate" begin
     particle_spacing = 0.01
 

test/systems/open_boundary_system.jl

@@ -97,4 +97,41 @@
 
         @test repr("text/plain", system) == show_box
     end
+
+    @testset "boundary zone width" begin
+        particle_spacing = 0.05
+        smoothing_kernel = WendlandC2Kernel{2}()
+        smoothing_length = 0.1
+
+        fluid = RectangularShape(particle_spacing, (2, 2), (0.0, 0.0), density=1.0)
+        fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                   smoothing_length, sound_speed=1.0)
+
+        boundary_zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
+                                     particle_spacing, face_normal=(1.0, 0.0),
+                                     density=1.0, open_boundary_layers=7,
+                                     boundary_type=InFlow())
+
+        compact_support = TrixiParticles.compact_support(smoothing_kernel, smoothing_length)
+        error_str = "boundary zone 1 has width $(boundary_zone.zone_width), " *
+                    "but must be at least two compact supports ($(2 * compact_support))"
+
+        @test_nowarn OpenBoundarySystem(boundary_zone; buffer_size=0,
+                                        boundary_model=BoundaryModelCharacteristicsLastiwka(),
+                                        fluid_system)
+
+        @test_nowarn OpenBoundarySystem(boundary_zone; buffer_size=0,
+                                        boundary_model=BoundaryModelDynamicalPressureZhang(),
+                                        fluid_system)
+
+        fluid_system_with_shifting = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                                 smoothing_length,
+                                                                 sound_speed=1.0,
+                                                                 shifting_technique=ParticleShiftingTechnique())
+
+        @test_throws ArgumentError(error_str) OpenBoundarySystem(boundary_zone;
+                                                                 buffer_size=0,
+                                                                 boundary_model=BoundaryModelDynamicalPressureZhang(),
+                                                                 fluid_system=fluid_system_with_shifting)
+    end
 end