Revise boundary zones and allow bidirectional flow (trixi-framework#623)

* fuse `InFlow` and `OutFlow` to `BoundaryZone`

* add docs

* adapt tests

* fix tests

* apply formatter

* add argument error

* quick fix

* implement suggestions

* switch from symbol to types

* extend words

* fix 3D pipe flow

* minor changes in docs

---------

Co-authored-by: LasNikas <niklas.nehe@web.de>
Co-authored-by: Sven Berger <berger.sven@gmail.com>
Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>

Author: 4 people

.typos.toml

@@ -7,4 +7,5 @@ Shepard = "Shepard"
 shepard = "shepard"
 Strack = "Strack"
 Lok = "Lok"
+PN = "PN"
 Typ = "Typ"

examples/fluid/pipe_flow_2d.jl

@@ -86,8 +86,9 @@ function velocity_function2d(pos, t)
 end
 
 plane_in = ([0.0, 0.0], [0.0, domain_size[2]])
-inflow = InFlow(; plane=plane_in, flow_direction,
-                open_boundary_layers, density=fluid_density, particle_spacing)
+inflow = BoundaryZone(; plane=plane_in,
+                      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(),
@@ -97,9 +98,9 @@ open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
                                          reference_velocity=velocity_function2d)
 
 plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]])
-outflow = OutFlow(; plane=plane_out,
-                  flow_direction, open_boundary_layers, density=fluid_density,
-                  particle_spacing)
+outflow = BoundaryZone(; plane=plane_out,
+                       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

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

src/general/semidiscretization.jl

@@ -855,3 +855,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 &&
+       boundary_zone isa BoundaryZone{BidirectionalFlow}
+        throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " *
+                            "Please specify inflow and outflow."))
+    end
+end

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -3,16 +3,18 @@
 
 Boundary model for `OpenBoundarySPHSystem`.
 This model uses the characteristic variables to propagate the appropriate values
-to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information
-about the method see [description below](@ref method_of_characteristics).
+to the outlet or inlet and have been proposed by Lastiwka et al. (2009).
+It requires a specific flow direction to be passed to the [`BoundaryZone`](@ref).
+For more information about the method see [description below](@ref method_of_characteristics).
 """
 struct BoundaryModelLastiwka end
 
 # Called from update callback via `update_open_boundary_eachstep!`
 @inline function update_boundary_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)
 
@@ -125,8 +127,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)
@@ -175,15 +178,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 transmitted 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/method_of_characteristics.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,
@@ -116,8 +115,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,
@@ -126,7 +123,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
@@ -143,12 +140,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
 
@@ -163,8 +161,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))
@@ -258,7 +255,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)
 
@@ -267,7 +264,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
 
@@ -282,6 +279,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)
+    relative_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 into the fluid domain.
+    if signbit(dot(relative_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/schemes/boundary/open_boundary/system.jl

@@ -345,9 +345,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)

src/visualization/write2vtk.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(),