Revise boundary zones and allow bidirectional flow

examples/fluid/pipe_flow_2d.jl

@@ -84,8 +84,9 @@ function velocity_function(pos, t)
  return SVector(prescribed_velocity, 0.0)
 end
 
-inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
- open_boundary_layers, density=fluid_density, particle_spacing)
+inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, domain_size[2]]),
+ plane_normal=flow_direction, open_boundary_layers,
+ density=fluid_density, particle_spacing, boundary_type=:inflow)
 
 open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
  boundary_model=BoundaryModelLastiwka(),
@@ -94,9 +95,9 @@ open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
  reference_pressure=pressure,
  reference_velocity=velocity_function)
 
-outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
- flow_direction, open_boundary_layers, density=fluid_density,
- particle_spacing)
+outflow = BoundaryZone(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
+  plane_normal=-flow_direction, open_boundary_layers,
+  density=fluid_density, particle_spacing, boundary_type=:outflow)
 
 open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
  boundary_model=BoundaryModelLastiwka(),

src/TrixiParticles.jl

@@ -54,8 +54,7 @@ include("visualization/recipes_plots.jl")
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
- BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow,
- OutFlow
+ BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, BoundaryZone
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
  PostprocessCallback, StepsizeCallback, UpdateCallback
 export ContinuityDensity, SummationDensity

src/general/semidiscretization.jl

@@ -854,3 +854,13 @@ function check_configuration(system::TotalLagrangianSPHSystem, systems)
  "`ContinuityDensity` is not yet supported for a `TotalLagrangianSPHSystem`"))
  end
 end
+
+function check_configuration(system::OpenBoundarySPHSystem, systems)
+ (; boundary_model, boundary_zone) = system
+
+ if boundary_model isa BoundaryModelLastiwka &&
+ first(typeof(boundary_zone).parameters) === BidirectionalFlow
+ throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " *
+ "Please specify inflow and outflow."))
+ end
+end

src/schemes/boundary/open_boundary/method_of_characteristics.jl

@@ -10,8 +10,9 @@ struct BoundaryModelLastiwka end
 
 @inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka,
  v, u, v_ode, u_ode, semi, t)
- (; density, pressure, cache, flow_direction,
+ (; density, pressure, cache, boundary_zone,
  reference_velocity, reference_pressure, reference_density) = system
+ (; flow_direction) = boundary_zone
 
  sound_speed = system_sound_speed(system.fluid_system)
 
@@ -114,8 +115,9 @@ end
 
 function evaluate_characteristics!(system, neighbor_system::FluidSystem,
  v, u, v_ode, u_ode, semi, t)
- (; volume, cache, flow_direction, density, pressure,
+ (; volume, cache, boundary_zone, density, pressure,
  reference_velocity, reference_pressure, reference_density) = system
+ (; flow_direction) = boundary_zone
  (; characteristics) = cache
 
  v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)
@@ -162,15 +164,15 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem,
  return system
 end
 
-@inline function prescribe_conditions!(characteristics, particle, ::OutFlow)
+@inline function prescribe_conditions!(characteristics, particle, ::BoundaryZone{OutFlow})
  # J3 is prescribed (i.e. determined from the exterior of the domain).
  # J1 and J2 is transimtted from the domain interior.
  characteristics[3, particle] = zero(eltype(characteristics))
 
  return characteristics
 end
 
-@inline function prescribe_conditions!(characteristics, particle, ::InFlow)
+@inline function prescribe_conditions!(characteristics, particle, ::BoundaryZone{InFlow})
  # Allow only J3 to propagate upstream to the boundary
  characteristics[1, particle] = zero(eltype(characteristics))
  characteristics[2, particle] = zero(eltype(characteristics))

src/schemes/boundary/open_boundary/system.jl

@@ -1,5 +1,5 @@
 @doc raw"""
- OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow};
+ OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
  fluid_system::FluidSystem, buffer_size::Integer,
  boundary_model,
  reference_velocity=nothing,
@@ -9,7 +9,7 @@
 Open boundary system for in- and outflow particles.
 
 # Arguments
-- `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary.
+- `boundary_zone`: See [`BoundaryZone`](@ref).
 
 # Keywords
 - `fluid_system`: The corresponding fluid system
@@ -39,15 +39,14 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  volume :: ARRAY1D # Array{ELTYPE, 1}: [particle]
  pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle]
  boundary_zone :: BZ
- flow_direction :: SVector{NDIMS, ELTYPE}
  reference_velocity :: RV
  reference_pressure :: RP
  reference_density :: RD
  buffer :: B
  update_callback_used :: Ref{Bool}
  cache :: C
 
- function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow};
+ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
  fluid_system::FluidSystem,
  buffer_size::Integer, boundary_model,
  reference_velocity=nothing,
@@ -98,8 +97,6 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  reference_density_ = wrap_reference_function(reference_density, Val(NDIMS))
  end
 
- flow_direction_ = boundary_zone.flow_direction
-
  cache = create_cache_open_boundary(boundary_model, initial_condition)
 
  return new{typeof(boundary_model), typeof(boundary_zone), NDIMS, ELTYPE,
@@ -108,7 +105,7 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  typeof(reference_density_), typeof(buffer),
  typeof(cache)}(initial_condition, fluid_system, boundary_model, mass,
  density, volume, pressure, boundary_zone,
- flow_direction_, reference_velocity_, reference_pressure_,
+ reference_velocity_, reference_pressure_,
  reference_density_, buffer, false, cache)
  end
 end
@@ -125,12 +122,13 @@ end
 
 timer_name(::OpenBoundarySPHSystem) = "open_boundary"
 vtkname(system::OpenBoundarySPHSystem) = "open_boundary"
+boundary_type_name(::BoundaryZone{ZT}) where {ZT} = string(nameof(ZT))
 
 function Base.show(io::IO, system::OpenBoundarySPHSystem)
  @nospecialize system # reduce precompilation time
 
  print(io, "OpenBoundarySPHSystem{", ndims(system), "}(")
- print(io, type2string(system.boundary_zone))
+ print(io, boundary_type_name(system.boundary_zone))
  print(io, ") with ", nparticles(system), " particles")
 end
 
@@ -145,8 +143,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem)
  summary_line(io, "#buffer_particles", system.buffer.buffer_size)
  summary_line(io, "fluid system", type2string(system.fluid_system))
  summary_line(io, "boundary model", type2string(system.boundary_model))
- summary_line(io, "boundary", type2string(system.boundary_zone))
- summary_line(io, "flow direction", system.flow_direction)
+ summary_line(io, "boundary type", boundary_type_name(system.boundary_zone))
  summary_line(io, "prescribed velocity", type2string(system.reference_velocity))
  summary_line(io, "prescribed pressure", type2string(system.reference_pressure))
  summary_line(io, "prescribed density", type2string(system.reference_density))
@@ -238,7 +235,7 @@ end
 
 # Outflow particle is outside the boundary zone
 @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
- boundary_zone::OutFlow, particle, v, u,
+ boundary_zone::BoundaryZone{OutFlow}, particle, v, u,
  v_fluid, u_fluid)
  deactivate_particle!(system, particle, u)
 
@@ -247,7 +244,7 @@ end
 
 # Inflow particle is outside the boundary zone
 @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
- boundary_zone::InFlow, particle, v, u,
+ boundary_zone::BoundaryZone{InFlow}, particle, v, u,
  v_fluid, u_fluid)
  (; spanning_set) = boundary_zone
 
@@ -262,6 +259,31 @@ end
  return system
 end
 
+# Buffer particle is outside the boundary zone
+@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
+ boundary_zone::BoundaryZone{BidirectionalFlow}, particle,
+ v, u, v_fluid, u_fluid)
+ particle_position = current_coords(u, system, particle) - boundary_zone.zone_origin
+
+ # Check if particle is in- or outside the fluid domain.
+ # `plane_normal` is always pointing in the fluid domain.
+ if signbit(dot(particle_position, boundary_zone.plane_normal))
+ deactivate_particle!(system, particle, u)
+
+ return system
+ end
+
+ # Activate a new particle in simulation domain
+ transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u)
+
+ # Reset position of boundary particle
+ for dim in 1:ndims(system)
+ u[dim, particle] = boundary_zone.plane_normal[dim]
+ end
+
+ return system
+end
+
 # Fluid particle is in boundary zone
 @inline function convert_particle!(fluid_system::FluidSystem, system,
  boundary_zone, particle, v, u, v_fluid, u_fluid)

src/visualization/write2vtk.jl

@@ -329,9 +329,8 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data
  for particle in active_particles(system)]
 
  if write_meta_data
- vtk["boundary_zone"] = type2string(system.boundary_zone)
+ vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters))
  vtk["width"] = round(system.boundary_zone.zone_width, digits=3)
- vtk["flow_direction"] = system.flow_direction
  vtk["velocity_function"] = type2string(system.reference_velocity)
  vtk["pressure_function"] = type2string(system.reference_pressure)
  vtk["density_function"] = type2string(system.reference_density)

test/general/buffer.jl

@@ -2,13 +2,14 @@
  # Mock fluid system
  struct FluidSystemMock3 <: TrixiParticles.FluidSystem{2, Nothing} end
 
- inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2,
- open_boundary_layers=2, density=1.0, flow_direction=[1.0, 0.0])
- system = OpenBoundarySPHSystem(inflow; fluid_system=FluidSystemMock3(),
+ zone = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2,
+ open_boundary_layers=2, density=1.0, plane_normal=[1.0, 0.0],
+ boundary_type=:inflow)
+ system = OpenBoundarySPHSystem(zone; fluid_system=FluidSystemMock3(),
  reference_density=0.0, reference_pressure=0.0,
  reference_velocity=[0, 0],
  boundary_model=BoundaryModelLastiwka(), buffer_size=0)
- system_buffer = OpenBoundarySPHSystem(inflow; buffer_size=5,
+ system_buffer = OpenBoundarySPHSystem(zone; buffer_size=5,
  reference_density=0.0, reference_pressure=0.0,
  reference_velocity=[0, 0],
  boundary_model=BoundaryModelLastiwka(),