Revise open boundaries

examples/fluid/pipe_flow_2d.jl

@@ -58,7 +58,7 @@ outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_siz
 
 NDIMS = ndims(pipe.fluid)
 
-n_buffer_particles = 5 * pipe.n_particles_per_dimension[2]^(NDIMS - 1)
+n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^(NDIMS - 1)
 
 # ==========================================================================================
 # ==== Fluid
@@ -73,12 +73,6 @@ kinematic_viscosity = maximum(prescribed_velocity) * domain_size[2] / reynolds_n
 
 viscosity = ViscosityAdami(nu=kinematic_viscosity)
 
-fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length,
-                                           sound_speed, viscosity=viscosity,
-                                           density_calculator=fluid_density_calculator,
-                                           shifting_technique=ParticleShiftingTechnique(),
-                                           buffer_size=n_buffer_particles)
-
 # Alternatively the WCSPH scheme can be used
 if wcsph
     state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
@@ -113,7 +107,7 @@ function velocity_function2d(pos, t)
     return SVector(prescribed_velocity)
 end
 
-open_boundary_model = BoundaryModelTafuniMirroring(; mirror_method=ZerothOrderMirroring())
+open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring())
 
 reference_velocity_in = velocity_function2d
 reference_pressure_in = nothing
@@ -131,7 +125,7 @@ reference_velocity_out = nothing
 reference_pressure_out = nothing
 reference_density_out = nothing
 boundary_type_out = OutFlow()
-plane_out = ([pipe.fluid_size[1], 0.0], [pipe.fluid_size[1], domain_size[2]])
+plane_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]])
 outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction),
                        open_boundary_layers, density=fluid_density, particle_spacing,
                        reference_density=reference_density_out,

src/TrixiParticles.jl

@@ -79,8 +79,8 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
        DensityDiffusionAntuono
 export tensile_instability_control
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
-       PressureMirroring, PressureZeroing, BoundaryModelLastiwkaCharacteristics,
-       BoundaryModelTafuniMirroring,
+       PressureMirroring, PressureZeroing, BoundaryModelCharacteristicsLastiwka,
+       BoundaryModelMirroringTafuni,
        BernoulliPressureExtrapolation
 export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring
 export HertzContactModel, LinearContactModel

src/general/semidiscretization.jl

@@ -979,9 +979,9 @@ function check_configuration(system::OpenBoundarySPHSystem, systems,
     fluid_system_index = findfirst(==(system.fluid_system), systems)
     system.fluid_system_index[] = fluid_system_index
 
-    if boundary_model isa BoundaryModelLastiwkaCharacteristics &&
+    if boundary_model isa BoundaryModelCharacteristicsLastiwka &&
        any(zone -> isnothing(zone.flow_direction), boundary_zones)
-        throw(ArgumentError("`BoundaryModelLastiwkaCharacteristics` needs a specific flow direction. " *
+        throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " *
                             "Please specify `InFlow()` and `OutFlow()`."))
     end
 

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -79,24 +79,35 @@ There are three ways to specify the actual shape of the boundary zone:
 !!! note "Note"
     The reference values (`reference_velocity`, `reference_pressure`, `reference_density`)
     can also be set to `nothing`.
-    In this case, they will either be extrapolated from the fluid domain ([BoundaryModelTafuniMirroring](@ref BoundaryModelTafuniMirroring))
-    or evolved using the characteristic flow variables ([BoundaryModelLastiwkaCharacteristics](@ref BoundaryModelLastiwkaCharacteristics)).
+    In this case, they will either be extrapolated from the fluid domain ([BoundaryModelMirroringTafuni](@ref BoundaryModelMirroringTafuni))
+    or evolved using the characteristic flow variables ([BoundaryModelCharacteristicsLastiwka](@ref BoundaryModelCharacteristicsLastiwka)).
 
 # Examples
 ```julia
 # 2D
 plane_points = ([0.0, 0.0], [0.0, 1.0])
 plane_normal=[1.0, 0.0]
 
+# Reference velocity function for the inflow: parabolic velocity profile
+velocity_in = (pos, t) -> SVector(4.0 * pos[2] * (1.0 - pos[2]), 0.0)
+
+# Reference pressure function for the inflow: y-dependent profile, sinusoidal in time
+pressure_in = (pos, t) -> pos[2] * sin(2pi * t)
+
 inflow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0,
-                      open_boundary_layers=4, boundary_type=InFlow())
+                      open_boundary_layers=4, boundary_type=InFlow(),
+                      reference_velocity=velocity_in, reference_pressure=pressure_in)
 
 # 3D
 plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
 plane_normal=[0.0, 0.0, 1.0]
 
+# Reference velocity function for the outflow: constant in flow direction
+velocity_out = [0.0, 0.0, -1.0]
+
 outflow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0,
-                       open_boundary_layers=4, boundary_type=OutFlow())
+                       open_boundary_layers=4, boundary_type=OutFlow(),
+                       reference_velocity=velocity_out)
 
 # 3D particles sampled as cylinder
 circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere())
@@ -109,13 +120,15 @@ bidirectional_flow = BoundaryZone(; plane=plane_points, plane_normal, particle_s
     This is an experimental feature and may change in any future releases.
 """
 struct BoundaryZone{IC, S, ZO, ZW, FD, PN, R}
-    initial_condition       :: IC
-    spanning_set            :: S
-    zone_origin             :: ZO
-    zone_width              :: ZW
-    flow_direction          :: FD
-    plane_normal            :: PN
-    reference_values        :: R
+    initial_condition :: IC
+    spanning_set      :: S
+    zone_origin       :: ZO
+    zone_width        :: ZW
+    flow_direction    :: FD
+    plane_normal      :: PN
+    reference_values  :: R
+    # Note that the following can't be static type parameters, as all boundary zones in a system
+    # must have the same type, so that we can loop over them in a type-stable way.
     average_inflow_velocity :: Bool
     prescribed_density      :: Bool
     prescribed_pressure     :: Bool
@@ -135,23 +148,10 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing,
     # `plane_normal` always points in fluid domain
     plane_normal_ = normalize(SVector(plane_normal...))
 
-    if boundary_type isa BidirectionalFlow
-        flow_direction = nothing
-
-    elseif boundary_type isa InFlow
-        # Unit vector pointing in downstream direction
-        flow_direction = plane_normal_
-
-    elseif boundary_type isa OutFlow
-        # Unit vector pointing in downstream direction
-        flow_direction = -plane_normal_
-    end
-
-    ic, spanning_set_, zone_origin,
-    zone_width = set_up_boundary_zone(plane, plane_normal_, flow_direction, density,
-                                      particle_spacing, initial_condition,
-                                      extrude_geometry, open_boundary_layers;
-                                      boundary_type=boundary_type)
+    ic, flow_direction, spanning_set_, zone_origin,
+    zone_width = set_up_boundary_zone(plane, plane_normal_, density, particle_spacing,
+                                      initial_condition, extrude_geometry,
+                                      open_boundary_layers, boundary_type)
 
     NDIMS = ndims(ic)
     ELTYPE = eltype(ic)
@@ -169,7 +169,7 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing,
             end
         end
         # We need this dummy for type stability reasons
-        velocity_dummy = SVector(ntuple(dim -> ELTYPE(Inf), NDIMS))
+        velocity_dummy = SVector(ntuple(dim -> convert(ELTYPE, Inf), NDIMS))
         velocity_ref = wrap_reference_function(reference_velocity, velocity_dummy)
     end
 
@@ -186,7 +186,7 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing,
             end
         end
         # We need this dummy for type stability reasons
-        pressure_dummy = ELTYPE(Inf)
+        pressure_dummy = convert(ELTYPE, Inf)
         pressure_ref = wrap_reference_function(reference_pressure, pressure_dummy)
     end
 
@@ -203,7 +203,7 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing,
             end
         end
         # We need this dummy for type stability reasons
-        density_dummy = ELTYPE(Inf)
+        density_dummy = convert(ELTYPE, Inf)
         density_ref = wrap_reference_function(reference_density, density_dummy)
     end
 
@@ -236,14 +236,13 @@ end
 function boundary_type_name(boundary_zone::BoundaryZone)
     (; flow_direction, plane_normal) = boundary_zone
 
-    # Bidirectional flow
-    isnothing(flow_direction) && return "bidirectional_flow"
-
-    # Outflow
-    signbit(dot(flow_direction, plane_normal)) && return "outflow"
-
-    # Inflow
-    return "inflow"
+    if isnothing(flow_direction)
+        return "bidirectional_flow"
+    elseif signbit(dot(flow_direction, plane_normal))
+        return "outflow"
+    else
+        return "inflow"
+    end
 end
 
 function Base.show(io::IO, boundary_zone::BoundaryZone)
@@ -267,29 +266,30 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone)
     end
 end
 
-function set_up_boundary_zone(plane, plane_normal, flow_direction, density,
-                              particle_spacing, initial_condition, extrude_geometry,
-                              open_boundary_layers; boundary_type)
+function set_up_boundary_zone(plane, plane_normal, density, particle_spacing,
+                              initial_condition, extrude_geometry, open_boundary_layers,
+                              boundary_type)
     if boundary_type isa InFlow
-        extrude_direction = -flow_direction
+        # Unit vector pointing in downstream direction
+        flow_direction = plane_normal
     elseif boundary_type isa OutFlow
-        extrude_direction = flow_direction
+        # Unit vector pointing in downstream direction
+        flow_direction = -plane_normal
     elseif boundary_type isa BidirectionalFlow
-        # `plane_normal` is always pointing in the fluid domain
-        extrude_direction = -plane_normal
+        flow_direction = nothing
     end
 
     # Sample particles in boundary zone
     if isnothing(initial_condition) && isnothing(extrude_geometry)
         initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
                                                             density,
-                                                            direction=extrude_direction,
+                                                            direction=(-plane_normal),
                                                             n_extrude=open_boundary_layers)
     elseif !isnothing(extrude_geometry)
         initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
                                                             particle_spacing,
                                                             density,
-                                                            direction=extrude_direction,
+                                                            direction=(-plane_normal),
                                                             n_extrude=open_boundary_layers)
     else
         initial_condition = initial_condition
@@ -310,23 +310,23 @@ function set_up_boundary_zone(plane, plane_normal, flow_direction, density,
         throw(ArgumentError("`plane_normal` is not normal to the boundary plane"))
     end
 
-    if boundary_type == InFlow()
+    if boundary_type isa InFlow
         # First vector of `spanning_vectors` is normal to the boundary plane
         dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction)
 
         # The vector must point in upstream direction for an inflow boundary.
         # Flip the normal vector to point in the opposite direction of `flow_direction`.
         spanning_set[:, 1] .*= -sign(dot_flow)
 
-    elseif boundary_type == OutFlow()
+    elseif boundary_type isa OutFlow
         # First vector of `spanning_vectors` is normal to the boundary plane
         dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction)
 
         # The vector must point in downstream direction for an outflow boundary.
         # Flip the normal vector to point in `flow_direction`.
         spanning_set[:, 1] .*= sign(dot_flow)
 
-    elseif boundary_type == BidirectionalFlow()
+    elseif boundary_type isa BidirectionalFlow
         # Flip the normal vector to point opposite to fluid domain
         spanning_set[:, 1] .*= -sign(dot_plane_normal)
     end
@@ -337,7 +337,7 @@ function set_up_boundary_zone(plane, plane_normal, flow_direction, density,
     # This check is only necessary when `initial_condition` or `extrude_geometry` are passed.
     ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin)
 
-    return ic, spanning_set_, zone_origin, zone_width
+    return ic, flow_direction, spanning_set_, zone_origin, zone_width
 end
 
 function calculate_spanning_vectors(plane, zone_width)
@@ -440,35 +440,53 @@ function wrap_reference_function(function_::Function, ref_dummy)
     return function_
 end
 
-function wrap_reference_function(constant_scalar_::Number, ref_dummy)
-    return @inline((coords, t)->constant_scalar_)
+function wrap_reference_function(constant_scalar::Number, ref_dummy)
+    return @inline((coords, t)->constant_scalar)
 end
 
-function wrap_reference_function(constant_vector_::AbstractVector,
+function wrap_reference_function(constant_vector::AbstractVector,
                                  ref_dummy::SVector{NDIMS, ELTYPE}) where {NDIMS, ELTYPE}
-    return @inline((coords, t)->SVector{NDIMS, ELTYPE}(constant_vector_))
+    return @inline((coords, t)->SVector{NDIMS, ELTYPE}(constant_vector))
 end
 
-function apply_reference_pressure(system, particle, pos, t)
-    (; pressure_references) = system.cache
+function reference_pressure(boundary_zone, v, system, particle, pos, t)
+    (; prescribed_pressure) = boundary_zone
+    (; pressure_reference_values) = system.cache
 
-    zone_id = system.boundary_zone_indices[particle]
+    if prescribed_pressure
+        zone_id = system.boundary_zone_indices[particle]
 
-    return apply_ith_function(pressure_references, zone_id, pos, t)
+        # `pressure_reference_values[zone_id](pos, t)`, but in a type-stable way
+        return apply_ith_function(pressure_reference_values, zone_id, pos, t)
+    else
+        return current_pressure(v, system, particle)
+    end
 end
 
-function apply_reference_density(system, particle, pos, t)
-    (; density_references) = system.cache
+function reference_density(boundary_zone, v, system, particle, pos, t)
+    (; prescribed_density) = boundary_zone
+    (; density_reference_values) = system.cache
 
-    zone_id = system.boundary_zone_indices[particle]
+    if prescribed_density
+        zone_id = system.boundary_zone_indices[particle]
 
-    return apply_ith_function(density_references, zone_id, pos, t)
+        # `density_reference_values[zone_id](pos, t)`, but in a type-stable way
+        return apply_ith_function(density_reference_values, zone_id, pos, t)
+    else
+        return current_density(v, system, particle)
+    end
 end
 
-function apply_reference_velocity(system, particle, pos, t)
-    (; velocity_references) = system.cache
+function reference_velocity(boundary_zone, v, system, particle, pos, t)
+    (; prescribed_velocity) = boundary_zone
+    (; velocity_reference_values) = system.cache
 
-    zone_id = system.boundary_zone_indices[particle]
+    if prescribed_velocity
+        zone_id = system.boundary_zone_indices[particle]
 
-    return apply_ith_function(velocity_references, zone_id, pos, t)
+        # `velocity_reference_values[zone_id](pos, t)`, but in a type-stable way
+        return apply_ith_function(velocity_reference_values, zone_id, pos, t)
+    else
+        return current_velocity(v, system, particle)
+    end
 end