Skip to content

Fix open-boundary characteristic zone handling #1203

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

Open
svchb wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix open-boundary characteristic zone handling #1203

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
3 changes: 3 additions & 0 deletions src/schemes/boundary/open_boundary/method_of_characteristics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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.
Expand Down
6 changes: 3 additions & 3 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 @@ -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))
Expand Down
63 changes: 63 additions & 0 deletions test/schemes/boundary/open_boundary/characteristic_variables.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]]
Expand Down
Loading