Merge branch 'main' into more_semi_cleanup

Author: svchb

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -162,18 +162,19 @@ bidirectional_flow = BoundaryZone(; boundary_face=face_vertices, face_normal,
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 """
-struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R, C}
+struct BoundaryZone{NDIMS, ELTYPE, IC, S, R, C}
     initial_condition :: IC
     spanning_set      :: S
-    zone_origin       :: ZO
-    zone_width        :: ZW
-    flow_direction    :: FD
-    face_normal       :: FN
+    zone_origin       :: SVector{NDIMS, ELTYPE}
+    zone_width        :: ELTYPE
+    flow_direction    :: SVector{NDIMS, ELTYPE}
+    face_normal       :: SVector{NDIMS, ELTYPE}
     rest_pressure     :: ELTYPE # Only required for `BoundaryModelDynamicalPressureZhang`
     reference_values  :: R
     cache             :: C
     # 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.
+    is_bidirectional        :: Bool
     average_inflow_velocity :: Bool
     prescribed_density      :: Bool
     prescribed_pressure     :: Bool
@@ -202,7 +203,7 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
     NDIMS = ndims(ic)
     ELTYPE = eltype(ic)
     if !(reference_velocity isa Function || isnothing(reference_velocity) ||
-         (reference_velocity isa Vector && length(reference_velocity) == NDIMS))
+         (reference_velocity isa AbstractVector && length(reference_velocity) == NDIMS))
         throw(ArgumentError("`reference_velocity` must be either a function mapping " *
                             "each particle's coordinates and time to its velocity, " *
                             "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity"))
@@ -281,16 +282,18 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
     cache = (;
              create_cache_boundary_zone(ic, boundary_face, face_normal_, sample_points)...)
 
+    is_bidirectional = boundary_type isa BidirectionalFlow
+
     return BoundaryZone(ic, spanning_set_, zone_origin, zone_width,
                         flow_direction, face_normal_, rest_pressure, reference_values,
-                        cache, average_inflow_velocity, prescribed_density,
-                        prescribed_pressure, prescribed_velocity)
+                        cache, is_bidirectional, average_inflow_velocity,
+                        prescribed_density, prescribed_pressure, prescribed_velocity)
 end
 
 function boundary_type_name(boundary_zone::BoundaryZone)
-    (; flow_direction, face_normal) = boundary_zone
+    (; flow_direction, face_normal, is_bidirectional) = boundary_zone
 
-    if isnothing(flow_direction)
+    if is_bidirectional
         return "bidirectional_flow"
     elseif signbit(dot(flow_direction, face_normal))
         return "outflow"
@@ -378,7 +381,7 @@ function set_up_boundary_zone(boundary_face, face_normal, density, particle_spac
         # Unit vector pointing in downstream direction
         flow_direction = -face_normal
     elseif boundary_type isa BidirectionalFlow
-        flow_direction = nothing
+        flow_direction = zero(face_normal)
     end
 
     # Sample particles in boundary zone
@@ -565,7 +568,7 @@ function wrap_reference_function(constant_vector::AbstractVector,
     return @inline((coords, t)->SVector{NDIMS, ELTYPE}(constant_vector))
 end
 
-function reference_pressure(boundary_zone, v, system, particle, pos, t)
+@inline function reference_pressure(boundary_zone, v, system, particle, pos, t)
     (; prescribed_pressure) = boundary_zone
     (; pressure_reference_values) = system.cache
 
@@ -579,7 +582,7 @@ function reference_pressure(boundary_zone, v, system, particle, pos, t)
     end
 end
 
-function reference_density(boundary_zone, v, system, particle, pos, t)
+@inline function reference_density(boundary_zone, v, system, particle, pos, t)
     (; prescribed_density) = boundary_zone
     (; density_reference_values) = system.cache
 
@@ -593,7 +596,7 @@ function reference_density(boundary_zone, v, system, particle, pos, t)
     end
 end
 
-function reference_velocity(boundary_zone, v, system, particle, pos, t)
+@inline function reference_velocity(boundary_zone, v, system, particle, pos, t)
     (; prescribed_velocity) = boundary_zone
     (; velocity_reference_values) = system.cache
 

src/schemes/boundary/open_boundary/dynamical_pressure.jl

@@ -98,9 +98,9 @@ function write_v0!(v0, system::OpenBoundarySystem, ::BoundaryModelDynamicalPress
     return v0
 end
 
-function reference_pressure(boundary_zone, v,
-                            system::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang},
-                            particle, pos, t)
+@inline function reference_pressure(boundary_zone, v,
+                                    system::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang},
+                                    particle, pos, t)
     (; prescribed_pressure, rest_pressure) = boundary_zone
     (; pressure_reference_values) = system.cache
 

src/schemes/boundary/open_boundary/mirroring.jl

@@ -486,10 +486,11 @@ end
 
 # Only for inflow boundary zones
 function average_velocity!(v, u, system, boundary_zone, semi)
-    (; face_normal, zone_origin, initial_condition, flow_direction) = boundary_zone
+    (; face_normal, zone_origin, initial_condition, flow_direction,
+     is_bidirectional) = boundary_zone
 
     # Bidirectional flow
-    isnothing(flow_direction) && return v
+    is_bidirectional && return v
 
     # Outflow
     signbit(dot(flow_direction, face_normal)) && return v
@@ -526,10 +527,10 @@ end
 # Only for inflow boundary zones
 function project_velocity_on_face_normal!(v, system, particle, boundary_zone,
                                           boundary_model)
-    (; face_normal, flow_direction) = boundary_zone
+    (; face_normal, flow_direction, is_bidirectional) = boundary_zone
 
     # Bidirectional flow
-    isnothing(flow_direction) && return v
+    is_bidirectional && return v
 
     # Outflow
     signbit(dot(flow_direction, face_normal)) && return v

src/schemes/boundary/open_boundary/system.jl

@@ -114,20 +114,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
     # in the `BoundaryZone`, but they are not used in the actual simulation.
     # The reference values are extracted above in the "create cache" function
     # and then stored in `system.cache` as a `Tuple`.
-    boundary_zones_new = map(zone -> BoundaryZone(zone.initial_condition,
-                                                  zone.spanning_set,
-                                                  zone.zone_origin,
-                                                  zone.zone_width,
-                                                  zone.flow_direction,
-                                                  zone.face_normal,
-                                                  zone.rest_pressure,
-                                                  nothing,
-                                                  zone.cache,
-                                                  zone.average_inflow_velocity,
-                                                  zone.prescribed_density,
-                                                  zone.prescribed_pressure,
-                                                  zone.prescribed_velocity),
-                             boundary_zones_)
+    boundary_zones_new = map(zone -> @set(zone.reference_values = nothing), boundary_zones_)
 
     return OpenBoundarySystem(boundary_model, initial_conditions, fluid_system,
                               fluid_system_index, smoothing_kernel, smoothing_length, mass,
@@ -297,6 +284,11 @@ end
     return v
 end
 
+function calculate_dt(v_ode, u_ode, cfl_number, system::OpenBoundarySystem, semi)
+    # Open boundaries don't affect the timestep calculation
+    return Inf
+end
+
 @inline function add_velocity!(du, v, u, particle, system::OpenBoundarySystem, t)
     boundary_zone = current_boundary_zone(system, particle)