diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 422a0f2e26..473f0f27f2 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -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 @@ -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 diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 7a56b2cfcc..da5dc2ab33 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -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) @@ -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) @@ -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 @@ -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) @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 @@ -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 diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index 9332080505..66f81e3d20 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -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