From 700ebcda525d66908f0c7cc048126d7bd411006c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 May 2026 00:54:34 +0200 Subject: [PATCH] Fix open-boundary characteristic zone handling --- .../method_of_characteristics.jl | 3 + src/schemes/boundary/open_boundary/system.jl | 6 +- .../open_boundary/characteristic_variables.jl | 63 +++++++++++++++++++ 3 files changed, 69 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index ce04478a0c..5023b9023f 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -169,6 +169,7 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) # Particle is outside of the influence of fluid particles. # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1. if volume[particle] * smoothing_length^ndims(system) < eps(eltype(smoothing_length)) + zone_id = system.boundary_zone_indices[particle] # Using the average of the values at the previous time step for particles which # are outside of the influence of fluid particles. @@ -178,6 +179,8 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) counter = 0 for neighbor in each_integrated_particle(system) + system.boundary_zone_indices[neighbor] == zone_id || continue + # Make sure that only neighbors in the influence of # the fluid particles are used. # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1. diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 3e4754d4ad..bfc2d9b462 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -699,9 +699,9 @@ function check_configuration(system::OpenBoundarySystem, systems, system.fluid_system_index[] = fluid_system_index if boundary_model isa BoundaryModelCharacteristicsLastiwka && - any(zone -> isnothing(zone.flow_direction), boundary_zones) - throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " * - "Please specify `InFlow()` and `OutFlow()`.")) + any(zone -> zone.is_bidirectional, boundary_zones) + throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a directed boundary zone. " * + "Please specify `InFlow()` or `OutFlow()` instead of `BidirectionalFlow()`.")) end if first(PointNeighbors.requires_update(neighborhood_search)) diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 9d2e9fa03e..512453a792 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -18,6 +18,69 @@ # Add small offset to avoid "ArgumentError: density must be positive and larger than `eps()`" reference_density = (pos, t) -> 1000.0 * (t + sqrt(eps())) + @testset "Reject bidirectional flow" begin + initial_condition = rectangular_patch(particle_spacing, (2, 2)) + fluid_system = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=ContinuityDensity(), + state_equation=nothing) + + bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), + face_normal=[1.0, 0.0], open_boundary_layers, + density, particle_spacing) + boundary_system = OpenBoundarySystem(bidirectional; fluid_system, buffer_size=0, + boundary_model=BoundaryModelCharacteristicsLastiwka()) + + @test_throws ArgumentError Semidiscretization(fluid_system, boundary_system) + end + + @testset "Fallback is zone-local" begin + face_vertices = ([0.0, 0.0], [0.0, 0.5]) + face_vertices_far = ([10.0, 0.0], [10.0, 0.5]) + flow_direction = SVector(1.0, 0.0) + + inflow = BoundaryZone(; boundary_face=face_vertices, face_normal=flow_direction, + open_boundary_layers, boundary_type=InFlow(), + reference_velocity, reference_pressure, reference_density, + density, particle_spacing) + inflow_far = BoundaryZone(; boundary_face=face_vertices_far, + face_normal=flow_direction, + open_boundary_layers, boundary_type=InFlow(), + reference_velocity, reference_pressure, reference_density, + density, particle_spacing) + fluid = extrude_geometry(face_vertices; particle_spacing, n_extrude=4, + density, pressure, direction=flow_direction) + fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel, + smoothing_length, + sound_speed, + buffer_size=0, + density_calculator=ContinuityDensity()) + boundary_system = OpenBoundarySystem(inflow, inflow_far; + fluid_system, buffer_size=0, + boundary_model=BoundaryModelCharacteristicsLastiwka()) + semi = Semidiscretization(fluid_system, boundary_system) + ode = semidiscretize(semi, (0.0, 5.0)) + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, boundary_system, semi) + u = TrixiParticles.wrap_u(u0_ode, boundary_system, semi) + + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, 2.0) + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, 3.0) + + zone_1_particles = findall(==(1), boundary_system.boundary_zone_indices) + zone_2_particles = findall(==(2), boundary_system.boundary_zone_indices) + + @test any(!isapprox(characteristic, 0.0) + for characteristic in boundary_system.cache.characteristics[:, + zone_1_particles]) + @test all(isapprox(characteristic, 0.0) + for characteristic in boundary_system.cache.characteristics[:, + zone_2_particles]) + end + # Face vertices of open boundary face_vertices_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] face_vertices_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]]