Skip to content

Fix open boundaries with periodic boundary zone #1237

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
efaulhaber wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Fix open boundaries with periodic boundary zone #1237

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
cd558a9
Fix open boundaries with periodic inflow or outflow
efaulhaber Jun 23, 2026
2c47280
Remove @info
efaulhaber Jun 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 7 additions & 3 deletions src/schemes/boundary/open_boundary/boundary_zones.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -487,8 +487,9 @@ end
span_dim = spanning_set[dim]
# Checks whether the projection of the particle position
# falls within the range of the zone
if !(0 <= dot(particle_position, span_dim) <= dot(span_dim, span_dim))

dot_span_dim = dot(span_dim, span_dim)
almostzero = 10 * eps(dot_span_dim)
if !(-almostzero <= dot(particle_position, span_dim) <= dot_span_dim + almostzero)
# Particle is not in boundary zone
return false
end
Expand All @@ -499,10 +500,13 @@ end
end

function update_boundary_zone_indices!(system, u, boundary_zones, semi)
periodic_box = get_neighborhood_search(system.fluid_system, system, semi).periodic_box

set_zero!(system.boundary_zone_indices)

@threaded semi for particle in each_integrated_particle(system)
particle_coords = current_coords(u, system, particle)
particle_coords = PointNeighbors.periodic_coords(current_coords(u, system, particle),
periodic_box)

for (zone_id, boundary_zone) in enumerate(boundary_zones)
# Check if boundary particle is in the boundary zone
Expand Down
33 changes: 22 additions & 11 deletions src/schemes/boundary/open_boundary/system.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -392,12 +392,14 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)

u_fluid = wrap_u(u_ode, fluid_system, semi)
v_fluid = wrap_v(v_ode, fluid_system, semi)
periodic_box = get_neighborhood_search(fluid_system, system, semi).periodic_box

boundary_candidates .= false

# Check the boundary particles whether they're leaving the boundary zone
@threaded semi for particle in each_integrated_particle(system)
particle_coords = current_coords(u, system, particle)
particle_coords = PointNeighbors.periodic_coords(current_coords(u, system, particle),
periodic_box)

# Check if boundary particle is outside the boundary zone
boundary_zone = current_boundary_zone(system, particle)
Expand All @@ -418,7 +420,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)

boundary_zone = current_boundary_zone(system, particle)
convert_particle!(system, fluid_system, boundary_zone, particle, particle_new,
v, u, v_fluid, u_fluid)
v, u, v_fluid, u_fluid, periodic_box)
end

update_system_buffer!(system.buffer)
Expand All @@ -428,7 +430,9 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)

# Check the fluid particles whether they're entering the boundary zone
@threaded semi for fluid_particle in each_integrated_particle(fluid_system)
fluid_coords = current_coords(u_fluid, fluid_system, fluid_particle)
fluid_coords = PointNeighbors.periodic_coords(current_coords(u_fluid, fluid_system,
fluid_particle),
periodic_box)

# Check if fluid particle is in any boundary zone
for boundary_zone in boundary_zones
Expand All @@ -449,7 +453,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)
particle_new = available_boundary_particles[i]

convert_particle!(fluid_system, system, particle, particle_new,
v, u, v_fluid, u_fluid)
v, u, v_fluid, u_fluid, periodic_box)
end

update_system_buffer!(system.buffer)
Expand All @@ -466,7 +470,7 @@ end
# Buffer particle is outside the boundary zone
@inline function convert_particle!(system::OpenBoundarySystem, fluid_system,
boundary_zone, particle, particle_new,
v, u, v_fluid, u_fluid)
v, u, v_fluid, u_fluid, periodic_box)
# Position relative to the origin of the transition face
relative_position = current_coords(u, system, particle) - boundary_zone.zone_origin

Expand All @@ -483,7 +487,8 @@ end
end

# Activate a new particle in simulation domain
transfer_particle!(fluid_system, system, particle, particle_new, v_fluid, u_fluid, v, u)
transfer_particle!(fluid_system, system, particle, particle_new,
v_fluid, u_fluid, v, u, periodic_box)

# Reset position of boundary particle back to the beginning of the boundary zone.
# If we translated it by exactly `zone_width` along `-face_normal`, rounding
Expand All @@ -497,7 +502,8 @@ end
end

# Verify the particle remains inside the boundary zone after the reset; deactivate it if not.
particle_coords = current_coords(u, system, particle)
particle_coords = PointNeighbors.periodic_coords(current_coords(u, system, particle),
periodic_box)
if !is_in_boundary_zone(boundary_zone, particle_coords)
deactivate_particle!(system, particle, v, u)

Expand All @@ -513,9 +519,11 @@ end

# Fluid particle is in boundary zone
@inline function convert_particle!(fluid_system::AbstractFluidSystem, system,
particle, particle_new, v, u, v_fluid, u_fluid)
particle, particle_new, v, u, v_fluid, u_fluid,
periodic_box)
# Activate particle in boundary zone
transfer_particle!(system, fluid_system, particle, particle_new, v, u, v_fluid, u_fluid)
transfer_particle!(system, fluid_system, particle, particle_new,
v, u, v_fluid, u_fluid, periodic_box)

# Deactivate particle in interior domain
deactivate_particle!(fluid_system, particle, v_fluid, u_fluid)
Expand All @@ -524,7 +532,7 @@ end
end

@inline function transfer_particle!(system_new, system_old, particle_old, particle_new,
v_new, u_new, v_old, u_old)
v_new, u_new, v_old, u_old, periodic_box)
# Activate new particle
system_new.buffer.active_particle[particle_new] = true

Expand All @@ -536,9 +544,12 @@ end
pressure = current_pressure(v_old, system_old, particle_old)
set_particle_pressure!(v_new, system_new, particle_new, pressure)

particle_coords = current_coords(u_old, system_old, particle_old)
particle_coords = PointNeighbors.periodic_coords(particle_coords, periodic_box)

# Exchange position and velocity
for dim in 1:ndims(system_new)
u_new[dim, particle_new] = u_old[dim, particle_old]
u_new[dim, particle_new] = particle_coords[dim]
v_new[dim, particle_new] = v_old[dim, particle_old]
end

Expand Down
55 changes: 55 additions & 0 deletions test/systems/open_boundary_system.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,4 +134,59 @@
boundary_model=BoundaryModelDynamicalPressureZhang(),
fluid_system=fluid_system_with_shifting)
end

@testset "periodic boundary conversion" begin
particle_spacing = 0.1
density = 1000.0

fluid = RectangularShape(particle_spacing, (20, 10), (0.0, 0.0);
density)
smoothing_kernel = WendlandC2Kernel{2}()
fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
smoothing_length=1.5 *
particle_spacing,
sound_speed=10.0, buffer_size=1)

outflow = BoundaryZone(; boundary_face=([2.0, 0.0], [2.0, 1.0]),
face_normal=(-1.0, 0.0), particle_spacing,
density, open_boundary_layers=1,
boundary_type=OutFlow())
open_boundary = OpenBoundarySystem(outflow; fluid_system,
boundary_model=BoundaryModelMirroringTafuni(),
buffer_size=1)

min_corner = [0.0, 0.0]
max_corner = [2.1, 1.0]
periodic_box = PeriodicBox(; min_corner, max_corner)
neighborhood_search = GridNeighborhoodSearch{2}(;
cell_list=FullGridCellList(;
min_corner,
max_corner),
update_strategy=ParallelUpdate(),
periodic_box)
semi = Semidiscretization(fluid_system, open_boundary;
neighborhood_search)
ode = semidiscretize(semi, (0.0, 1.0))
v_ode, u_ode = ode.u0.x
TrixiParticles.initialize!(open_boundary, semi)

u_fluid = TrixiParticles.wrap_u(u_ode, fluid_system, semi)
v_open_boundary = TrixiParticles.wrap_v(v_ode, open_boundary, semi)
u_open_boundary = TrixiParticles.wrap_u(u_ode, open_boundary, semi)

particle = findfirst(i -> u_fluid[1, i] > 1.9 && u_fluid[2, i] < 0.1,
axes(u_fluid, 2))
particle_new = findfirst(==(false), open_boundary.buffer.active_particle)

u_fluid[:, particle] .= [2.05, -eps(Float64)]

TrixiParticles.check_domain!(open_boundary, v_open_boundary, u_open_boundary,
v_ode, u_ode, semi)

@test !fluid_system.buffer.active_particle[particle]
@test open_boundary.buffer.active_particle[particle_new]
@test TrixiParticles.is_in_boundary_zone(outflow,
SVector(u_open_boundary[:,
particle_new]...))
end
end
Loading