Skip to content

Commit

Permalink
remove redundancy
Browse files Browse the repository at this point in the history
  • Loading branch information
@svchb
svchb committed Oct 10, 2024
1 parent 8c72576 commit 782f9de
Showing 1 changed file with 106 additions and 122 deletions.
228 changes: 106 additions & 122 deletions src/schemes/boundary/open_boundary/boundary_zones.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,76 @@
abstract type FlowBoundary{NDIMS} end

struct InFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS}
initial_condition :: IC
spanning_set :: S
zone_origin :: ZO
zone_width :: ZW
flow_direction :: FD
end

struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS}
initial_condition :: IC
spanning_set :: S
zone_origin :: ZO
zone_width :: ZW
flow_direction :: FD
end

function _create_flow_boundary(; plane, flow_direction, density, particle_spacing,
initial_condition, extrude_geometry,
open_boundary_layers::Integer, is_inflow::Bool)
if open_boundary_layers <= 0
throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero"))
end

# Unit vector pointing in downstream direction
flow_direction_ = normalize(SVector{length(flow_direction)}(flow_direction))

# Determine extrusion direction
extrusion_direction = is_inflow ? -flow_direction_ : flow_direction_

# Sample particles in boundary zone
if isnothing(initial_condition) && isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
density,
direction=extrusion_direction,
n_extrude=open_boundary_layers)
elseif !isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
particle_spacing,
density,
direction=extrusion_direction,
n_extrude=open_boundary_layers)
end

NDIMS = ndims(initial_condition)
ELTYPE = eltype(initial_condition)

zone_width = open_boundary_layers * initial_condition.particle_spacing

# Vectors spanning the boundary zone/box
spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width)

# First vector of `spanning_set` is normal to the plane.
# The normal vector must point in the correct direction.
normal_vector = normalize(spanning_set[:, 1])
dot_ = dot(normal_vector, flow_direction_)

if !isapprox(abs(dot_), 1.0, atol=1e-7)
throw(ArgumentError("`flow_direction` must be normal to the plane defined by `plane` points."))
end

# Flip the normal vector to point in the correct direction
spanning_set[:, 1] .*= is_inflow ? -sign(dot_) : sign(dot_)

spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set)

# Remove particles outside the boundary zone
ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin)

return NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_
end

@doc raw"""
InFlow(; plane, flow_direction, density, particle_spacing,
initial_condition=nothing, extrude_geometry=nothing,
Expand Down Expand Up @@ -71,67 +144,22 @@ inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, dens
!!! warning "Experimental Implementation"
This is an experimental feature and may change in any future releases.
"""
struct InFlow{NDIMS, IC, S, ZO, ZW, FD}
initial_condition :: IC
spanning_set :: S
zone_origin :: ZO
zone_width :: ZW
flow_direction :: FD

function InFlow(; plane, flow_direction, density, particle_spacing,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
if open_boundary_layers <= 0
throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero"))
end

# Unit vector pointing in downstream direction
flow_direction_ = normalize(SVector(flow_direction...))

# Sample particles in boundary zone
if isnothing(initial_condition) && isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
density,
direction=-flow_direction_,
n_extrude=open_boundary_layers)
elseif !isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
particle_spacing,
density,
direction=-flow_direction_,
n_extrude=open_boundary_layers)
end

NDIMS = ndims(initial_condition)
ELTYPE = eltype(initial_condition)

zone_width = open_boundary_layers * initial_condition.particle_spacing

# Vectors spanning the boundary zone/box
spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width)

# First vector of `spanning_vectors` is normal to the inflow plane.
# The normal vector must point in upstream direction for an inflow boundary.
dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_)

if !isapprox(abs(dot_), 1.0, atol=1e-7)
throw(ArgumentError("`flow_direction` is not normal to inflow plane"))
else
# Flip the normal vector to point in the opposite direction of `flow_direction`
spanning_set[:, 1] .*= -sign(dot_)
end

spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set)

# Remove particles outside the boundary zone.
# This check is only necessary when `initial_condition` or `extrude_geometry` are passed.
ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin)

return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin),
typeof(zone_width),
typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width,
flow_direction_)
end
function InFlow(; plane, flow_direction, density, particle_spacing,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(;
plane,
flow_direction,
density,
particle_spacing,
initial_condition,
extrude_geometry,
open_boundary_layers,
is_inflow=true)
return InFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin),
typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_,
zone_origin, zone_width,
flow_direction_)
end

@doc raw"""
Expand Down Expand Up @@ -207,66 +235,22 @@ outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, de
!!! warning "Experimental Implementation"
This is an experimental feature and may change in any future releases.
"""
struct OutFlow{NDIMS, IC, S, ZO, ZW, FD}
initial_condition :: IC
spanning_set :: S
zone_origin :: ZO
zone_width :: ZW
flow_direction :: FD

function OutFlow(; plane, flow_direction, density, particle_spacing,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
if open_boundary_layers <= 0
throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero"))
end

# Unit vector pointing in downstream direction
flow_direction_ = normalize(SVector(flow_direction...))

# Sample particles in boundary zone
if isnothing(initial_condition) && isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
density,
direction=flow_direction_,
n_extrude=open_boundary_layers)
elseif !isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
particle_spacing, density,
direction=flow_direction_,
n_extrude=open_boundary_layers)
end

NDIMS = ndims(initial_condition)
ELTYPE = eltype(initial_condition)

zone_width = open_boundary_layers * initial_condition.particle_spacing

# Vectors spanning the boundary zone/box
spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width)

# First vector of `spanning_vectors` is normal to the outflow plane.
# The normal vector must point in downstream direction for an outflow boundary.
dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_)

if !isapprox(abs(dot_), 1.0, atol=1e-7)
throw(ArgumentError("`flow_direction` is not normal to outflow plane"))
else
# Flip the normal vector to point in `flow_direction`
spanning_set[:, 1] .*= sign(dot_)
end

spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set)

# Remove particles outside the boundary zone.
# This check is only necessary when `initial_condition` or `extrude_geometry` are passed.
ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin)

return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin),
typeof(zone_width),
typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width,
flow_direction_)
end
function OutFlow(; plane, flow_direction, density, particle_spacing,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(;
plane,
flow_direction,
density,
particle_spacing,
initial_condition,
extrude_geometry,
open_boundary_layers,
is_inflow=false)
return OutFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin),
typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_,
zone_origin, zone_width,
flow_direction_)
end

@inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS
Expand Down Expand Up @@ -300,7 +284,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width)
return hcat(c, edge1, edge2)
end

@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords)
@inline function is_in_boundary_zone(boundary_zone::FlowBoundary, particle_coords)
(; zone_origin, spanning_set) = boundary_zone
particle_position = particle_coords - zone_origin

Expand Down

0 comments on commit 782f9de

Please sign in to comment.