diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 20ef2abee..3e744bcb3 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -70,18 +70,17 @@ function velocity_function(pos, t) return SVector(prescribed_velocity, 0.0) # SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) end -open_boundary_in = OpenBoundarySPHSystem(([0.0, 0.0], [0.0, domain_size[2]]), InFlow(), - sound_speed; particle_spacing, - flow_direction, open_boundary_layers, - density=fluid_density, buffer=n_buffer_particles, +inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing) + +open_boundary_in = OpenBoundarySPHSystem(inflow, sound_speed; buffer=n_buffer_particles, reference_pressure=pressure, reference_velocity=velocity_function) -open_boundary_out = OpenBoundarySPHSystem(([domain_size[1], 0.0], - [domain_size[1], domain_size[2]]), OutFlow(), - sound_speed; particle_spacing, - flow_direction, open_boundary_layers, - density=fluid_density, buffer=n_buffer_particles, +outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), + flow_direction, open_boundary_layers, density=fluid_density, + particle_spacing) +open_boundary_out = OpenBoundarySPHSystem(outflow, sound_speed; buffer=n_buffer_particles, reference_pressure=pressure, reference_velocity=velocity_function) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl new file mode 100644 index 000000000..cc8366471 --- /dev/null +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -0,0 +1,275 @@ +""" + InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer=0) + +Inflow boundary zone for [`OpenBoundarySPHSystem`](@ref) + +# Keywords +- `plane`: Points defining the boundary zones front plane. + The points must either span a rectangular plane in 3D or a line in 2D. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in upstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for boundary zone (optional). + Particles outside the boundary zone will be removed. +- `extrude_geometry=nothing`: Extrude a geometry inside the boundary zone (optional). + For further information see [`extrude_geometry`](@ref). + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction=[1.0, 0.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction=[0.0, 0.0, 1.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` +""" +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=0) + if open_boundary_layers < sqrt(eps()) + 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 and normal vector of " * + "inflow-plane do not correspond")) + else + # Flip the inflow vector correspondingly + spanning_set[:, 1] .*= -dot_ + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove paricles outside the boundary zone. + # This check is only necessary when custom `initial_condition` 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 +end + +""" + OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer=0) + +Outflow boundary zone for [`OpenBoundarySPHSystem`](@ref) + +# Keywords +- `plane`: Points defining the boundary zones front plane. + The points must either span a rectangular plane in 3D or a line in 2D. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in upstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for boundary zone (optional). + Particles outside the boundary zone will be removed. +- `extrude_geometry=nothing`: Extrude a geometry inside the boundary zone (optional). + For further information see [`extrude_geometry`](@ref). + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction = [1.0, 0.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction = [0.0, 0.0, 1.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` +""" +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=0) + if open_boundary_layers < sqrt(eps()) + 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 and normal vector of " * + "outflow-plane do not correspond")) + else + # Flip the inflow vector correspondingly + spanning_set[:, 1] .*= dot_ + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove paricles outside the boundary zone. + # This check is only necessary when custom `initial_condition` 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 +end + +@inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS + +# function calculate_spanning_vectors(plane::Shapes, zone_width) +# # TODO: Handle differently +# end + +function calculate_spanning_vectors(plane, zone_width) + return spanning_vectors(plane, zone_width), SVector(plane[1]...) +end + +function spanning_vectors(plane_points, zone_width) + + # Convert to tuple + return spanning_vectors(tuple(plane_points...), zone_width) +end + +function spanning_vectors(plane_points::NTuple{2}, zone_width) + plane_size = plane_points[2] - plane_points[1] + + # Calculate normal vector of plane + b = Vector(normalize([-plane_size[2]; plane_size[1]]) * zone_width) + + return hcat(b, plane_size) +end + +function spanning_vectors(plane_points::NTuple{3}, zone_width) + # Vectors spanning the plane + edge1 = plane_points[2] - plane_points[1] + edge2 = plane_points[3] - plane_points[1] + + if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) + throw(ArgumentError("the provided points do not span a rectangular plane")) + end + + # Calculate normal vector of plane + c = Vector(normalize(cross(edge2, edge1)) * zone_width) + + return hcat(c, edge1, edge2) +end + +function remove_outside_particles(initial_condition, spanning_set, zone_origin) + (; coordinates, density, particle_spacing) = initial_condition + + in_zone = trues(nparticles(initial_condition)) + + for particle in eachparticle(initial_condition) + current_position = current_coords(coordinates, initial_condition, particle) + particle_position = current_position - zone_origin + + for dim in 1:ndims(initial_condition) + span_dim = spanning_set[dim] + # Checks whether the projection of the particle position + # falls within the range of the zone. + if !(0 <= dot(particle_position, span_dim) <= dot(span_dim, span_dim)) + + # Particle is not in boundary zone. + in_zone[particle] = false + end + end + end + + return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), + particle_spacing) +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 24c8c058b..9bd0af42d 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,16 +1,4 @@ -""" - InFlow - -Inflow boundary zone for [`OpenBoundarySPHSystem`](@ref) -""" -struct InFlow end - -""" - OutFlow - -Outflow boundary zone for [`OpenBoundarySPHSystem`](@ref) -""" -struct OutFlow end +include("boundary_zones.jl") """ OpenBoundarySPHSystem(plane_points, boundary_zone::Union{InFlow, OutFlow}, @@ -25,22 +13,10 @@ to the outlet or inlet and has been proposed by Lastiwka et al (2009). For more about the method see [Open Boundary System](@ref open_boundary). # Arguments -- `plane_points`: Points defining the boundary zones front plane. - The points must either span a rectangular plane in 3D or a line in 2D. - See description above for more information. - `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary. - `sound_speed`: Speed of sound. # Keywords -- `sample_plane`: For customized particle sampling in the boundary zone, this can be either - points defining a 3D plane (2D line), particle coordinates defining a specific - shape or a specific [`InitialCondition`](@ref) type. - The geometry will be extruded in upstream direction with [`extrude_geometry`](@ref). - Default is `plane_points` which fully samples the boundary zone with particles. -- `particle_spacing`: The spacing between the particles in the boundary zone. -- `flow_direction`: Vector defining the flow direction. -- `open_boundary_layers`: Number of particle layers in upstream direction. -- `density`: Density of each particle to define the mass of each particle (see [`InitialCondition`](@ref)). - `buffer`: Number of buffer particles. - `reference_velocity`: Reference velocity is either a function mapping each particle's coordinates and time to its velocity, an array where the ``i``-th column holds @@ -54,32 +30,8 @@ about the method see [Open Boundary System](@ref open_boundary). and time to its density, a vector holding the density of each particle, or a scalar for a constant density over all particles. Density is constant zero by default. - -# Examples -```julia -# 2D inflow -plane_points = ([0.0, 0.0], [0.0, 1.0]) -flow_direction=[1.0, 0.0] - -system = OpenBoundarySPHSystem(plane_points, InFlow(), 10.0; particle_spacing=0.1, - open_boundary_layers=4, density=1.0, flow_direction) - -# 3D outflow -plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) -flow_direction=[0.0, 0.0, 1.0] - -system = OpenBoundarySPHSystem(plane_points, OutFlow(), 10.0; particle_spacing=0.1, - open_boundary_layers=4, density=1.0, flow_direction) - -# 3D particles sampled as cylinder -circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) - -system = OpenBoundarySPHSystem(plane_points, InFlow(), 10.0; particle_spacing=0.1, - sample_geometry=circle, - open_boundary_layers=4, density=1.0, flow_direction) -``` """ -struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, S, RV, RP, +struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, RV, RP, RD, B} <: System{NDIMS} initial_condition :: IC mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] @@ -91,35 +43,18 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, S, sound_speed :: ELTYPE boundary_zone :: BZ flow_direction :: SVector{NDIMS, ELTYPE} - zone_origin :: SVector{NDIMS, ELTYPE} - spanning_set :: S reference_velocity :: RV reference_pressure :: RP reference_density :: RD buffer :: B update_callback_used :: Ref{Bool} - function OpenBoundarySPHSystem(plane_points, boundary_zone::Union{InFlow, OutFlow}, - sound_speed; - sample_geometry=plane_points, particle_spacing, - flow_direction, open_boundary_layers::Integer=0, density, + function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}, sound_speed; buffer=nothing, - reference_velocity=zeros(length(flow_direction)), - reference_pressure=0.0, reference_density=density) - if open_boundary_layers < sqrt(eps()) - 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 boundary zone with particles in corresponding direction. - (boundary_zone isa OutFlow) && (direction = flow_direction_) - (boundary_zone isa InFlow) && (direction = -flow_direction_) - - # Sample particles in boundary zone. - initial_condition = extrude_geometry(sample_geometry; particle_spacing, direction, - n_extrude=open_boundary_layers, density) + reference_velocity=zeros(ndims(boundary_zone)), + reference_pressure=0.0, + reference_density=first(boundary_zone.initial_condition.density)) + (; initial_condition) = boundary_zone (buffer ≠ nothing) && (buffer = SystemBuffer(nparticles(initial_condition), buffer)) initial_condition = allocate_buffer(initial_condition, buffer) @@ -150,33 +85,6 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, S, reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end - zone_width = open_boundary_layers * initial_condition.particle_spacing - - # Vectors spanning the boundary zone/box. - spanning_set = spanning_vectors(plane_points, zone_width) - - # First vector of `spanning_vectors` is normal to the in-/outflow plane. - # The normal vector must point in downstream direction for an outflow boundary and - # for an inflow boundary the normal vector must point in upstream direction. - # Thus, rotate the normal vector correspondingly. - if isapprox(dot(normalize(spanning_set[:, 1]), flow_direction_), 1.0, atol=1e-7) - # Normal vector points in downstream direction. - # Flip the inflow vector in upstream direction - (boundary_zone isa InFlow) && (spanning_set[:, 1] .*= -1) - elseif isapprox(dot(normalize(spanning_set[:, 1]), flow_direction_), -1.0, - atol=1e-7) - # Normal vector points in upstream direction. - # Flip the outflow vector in downstream direction - (boundary_zone isa OutFlow) && (spanning_set[:, 1] .*= -1) - else - throw(ArgumentError("flow direction and normal vector of " * - "$(typeof(boundary_zone))-plane do not correspond")) - end - - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) - - zone_origin = SVector(plane_points[1]...) - mass = copy(initial_condition.mass) pressure = [reference_pressure_(initial_condition.coordinates[:, i], 0.0) for i in eachparticle(initial_condition)] @@ -186,15 +94,16 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, S, characteristics = zeros(ELTYPE, 3, length(mass)) previous_characteristics = zeros(ELTYPE, 3, length(mass)) + flow_direction_ = boundary_zone.flow_direction + return new{typeof(boundary_zone), NDIMS, ELTYPE, typeof(initial_condition), - typeof(mass), typeof(characteristics), typeof(spanning_set_), + typeof(mass), typeof(characteristics), typeof(reference_velocity_), typeof(reference_pressure_), typeof(reference_density_), typeof(buffer)}(initial_condition, mass, density, volume, pressure, characteristics, previous_characteristics, sound_speed, - boundary_zone, flow_direction_, zone_origin, - spanning_set_, reference_velocity_, reference_pressure_, - reference_density_, buffer, false) + boundary_zone, flow_direction_, reference_velocity_, + reference_pressure_, reference_density_, buffer, false) end end @@ -205,7 +114,7 @@ function Base.show(io::IO, system::OpenBoundarySPHSystem) @nospecialize system # reduce precompilation time print(io, "OpenBoundarySPHSystem{", ndims(system), "}(") - print(io, system.boundary_zone) + print(io, type2string(system.boundary_zone)) print(io, ") with ", nparticles(system), " particles") end @@ -222,12 +131,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) else summary_line(io, "#particles", nparticles(system)) end - summary_line(io, "boundary", system.boundary_zone) + summary_line(io, "boundary", type2string(system.boundary_zone)) summary_line(io, "flow direction", system.flow_direction) summary_line(io, "prescribed velocity", string(nameof(system.reference_velocity))) summary_line(io, "prescribed pressure", string(nameof(system.reference_pressure))) summary_line(io, "prescribed density", string(nameof(system.reference_density))) - summary_line(io, "width", round(norm(system.spanning_set[1]), digits=3)) + summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3)) summary_footer(io) end end @@ -246,41 +155,11 @@ end return system.pressure[particle] end -function spanning_vectors(plane_points, zone_width) - - # Convert to tuple - return spanning_vectors(tuple(plane_points...), zone_width) -end - -function spanning_vectors(plane_points::NTuple{2}, zone_width) - plane_size = plane_points[2] - plane_points[1] - - # Calculate normal vector of plane - b = Vector(normalize([-plane_size[2]; plane_size[1]]) * zone_width) - - return hcat(b, plane_size) -end - -function spanning_vectors(plane_points::NTuple{3}, zone_width) - # Vectors spanning the plane - edge1 = plane_points[2] - plane_points[1] - edge2 = plane_points[3] - plane_points[1] - - if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) - throw(ArgumentError("the provided points do not span a rectangular plane")) - end - - # Calculate normal vector of plane - c = Vector(normalize(cross(edge2, edge1)) * zone_width) - - return hcat(c, edge1, edge2) -end - -@inline function within_boundary_zone(particle_coords, system) - (; zone_origin, spanning_set) = system +@inline function (boundary_zone::Union{InFlow, OutFlow})(particle_coords) + (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin - for dim in 1:ndims(system) + for dim in 1:ndims(boundary_zone) span_dim = spanning_set[dim] # Checks whether the projection of the particle position # falls within the range of the zone. @@ -501,7 +380,7 @@ function check_fluid_domain!(system, fluid_system::FluidSystem, particle, neighborhood_search = get_neighborhood_search(system, fluid_system, semi) # Check if the particle position is outside the boundary zone. - if !within_boundary_zone(particle_coords, system) + if !boundary_zone(particle_coords) transform_particle!(system, fluid_system, boundary_zone, particle, v, u, v_fluid, u_fluid) end @@ -511,7 +390,7 @@ function check_fluid_domain!(system, fluid_system::FluidSystem, particle, fluid_coords = current_coords(u_fluid, fluid_system, neighbor) # Check if neighbor position is in boundary zone - if within_boundary_zone(fluid_coords, system) + if boundary_zone(fluid_coords) transform_particle!(fluid_system, system, boundary_zone, neighbor, v, u, v_fluid, u_fluid) end @@ -522,7 +401,8 @@ end # Outflow particle is outside the boundary zone @inline function transform_particle!(system::OpenBoundarySPHSystem, fluid_system, - ::OutFlow, particle, v, u, v_fluid, u_fluid) + boundary_zone::OutFlow, particle, v, u, + v_fluid, u_fluid) deactivate_particle!(system, particle, u) return system @@ -530,8 +410,9 @@ end # Inflow particle is outside the boundary zone @inline function transform_particle!(system::OpenBoundarySPHSystem, fluid_system, - ::InFlow, particle, v, u, v_fluid, u_fluid) - (; spanning_set) = system + boundary_zone::InFlow, particle, v, u, + v_fluid, u_fluid) + (; spanning_set) = boundary_zone # Activate a new particle in simulation domain transfer_quantities!(fluid_system, system, particle, v_fluid, u_fluid, v, u) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 3bdf98528..2acc4b0df 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -282,7 +282,7 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data if write_meta_data vtk["boundary_zone"] = type2string(system.boundary_zone) vtk["sound_speed"] = system.sound_speed - vtk["width"] = round(norm(system.spanning_set[1]), digits=3) + vtk["width"] = round(system.boundary_zone.zone_width, digits=3) vtk["flow_direction"] = system.flow_direction vtk["velocity_function"] = string(nameof(system.reference_velocity)) vtk["pressure_function"] = string(nameof(system.reference_pressure)) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl new file mode 100644 index 000000000..35b1076c3 --- /dev/null +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -0,0 +1,221 @@ +@testset verbose=true "Boundary Zone" begin + @testset verbose=true "Boundary Zone 2D" begin + particle_spacing = 0.2 + open_boundary_layers = 4 + + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset "Points $(i)" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + + plane_size = point_2 - point_1 + + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset "Flow Direction $(j)" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "$boundary_zone" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test sign_ * flow_directions[j] ≈ + normalize(boundary_zone.spanning_set[1]) + @test zone_width ≈ norm(boundary_zone.spanning_set[1]) + end + end + end + end + + @testset verbose=true "Boundary Zone 3D" begin + particle_spacing = 0.05 + open_boundary_layers = 4 + + plane_points_1 = [ + [0.0, 0.0, 0.0], + [0.3113730847835541, 0.19079485535621643, -0.440864622592926], + ] + plane_points_2 = [ + [1.0, 0.0, 0.0], + [-0.10468611121177673, 0.252103328704834, -0.44965094327926636], + ] + plane_points_3 = [ + [0.0, 1.0, 0.0], + [0.3113730847835541, 0.25057315826416016, -0.02374829351902008], + ] + + @testset "Points $(i)" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + point_3 = plane_points_3[i] + + edge1 = point_2 - point_1 + edge2 = point_3 - point_1 + + flow_directions = [ + normalize(cross(edge1, edge2)), + -normalize(cross(edge1, edge2)), + ] + + @testset "Flow Direction $(j)" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "$boundary_zone" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test plane_points_3[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[3] + @test sign_ * flow_directions[j] ≈ + normalize(boundary_zone.spanning_set[1]) + @test zone_width ≈ norm(boundary_zone.spanning_set[1]) + end + end + end + end + + @testset verbose=true "Particle In Boundary Zone 2D" begin + plane_points = [[-0.2, -0.5], [0.3, 0.6]] + plane_size = plane_points[2] - plane_points[1] + + flow_direction = normalize([-plane_size[2], plane_size[1]]) + + inflow = InFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? eps() : -eps() + + point_3 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin + + query_points = Dict( + "Behind" => ([-1.0, -1.0], false), + "Before" => ([2.0, 2.0], false), + "On Point 1" => (plane_points[1], true), + "On Point 2" => (plane_points[2], true), + "On Point 3" => (point_3, true), + "Closely On Point 1" => (plane_points[1] .+ perturb_ * flow_direction, + false)) + + TrixiParticles.@autoinfiltrate + + @testset "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == boundary_zone(particle_position) + end + end + end + + @testset verbose=true "Particle In Boundary Zone 3D" begin + point1 = [-0.2, -0.5, 0.0] + point2 = [0.3, 0.5, 0.0] + point3 = [0.13111173850909402, -0.665555869254547, 0.0] + + flow_direction = normalize(cross(point2 - point1, point3 - point1)) + + inflow = InFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? eps() : -eps() + query_points = Dict( + "Behind" => ([-1.0, -1.0, 1.2], false), + "Before" => ([2.0, 2.0, -1.2], false), + "On Point 1" => (point1, true), + "On Point 2" => (point2, true), + "On Point 3" => (point3, true), + "On Point 4" => (boundary_zone.spanning_set[1] + boundary_zone.zone_origin, + true), + "Closely On Point 1" => (point1 .+ perturb_ * flow_direction, false)) + + @testset "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == boundary_zone(particle_position) + end + end + end + + @testset verbose=true "Illegal Inputs" begin + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + flow_direction = [0.0, 0.0, 1.0] + + error_str = "the provided points do not span a rectangular plane" + + @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + @test_throws ArgumentError(error_str) OutFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + rectangular_plane = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + flow_direction = [0.0, 1.0, 0.0] + + error_str = "flow direction and normal vector of " * + "inflow-plane do not correspond" + + @test_throws ArgumentError(error_str) InFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + error_str = "flow direction and normal vector of " * + "outflow-plane do not correspond" + + @test_throws ArgumentError(error_str) OutFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + end +end diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index fafa0c72b..da535ce82 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -1,6 +1,5 @@ @testset verbose=true "Characteristic Variables" begin particle_spacing = 0.1 - boundary_zones = [InFlow(), OutFlow()] # Number of boundary particles in the influence of fluid particles influenced_particles = [20, 52, 26] @@ -19,29 +18,37 @@ reference_density = (pos, t) -> 1000.0 * t # Plane points of open boundary - point1s = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] - point2s = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] - - @testset "$boundary_zone" for boundary_zone in boundary_zones - @testset "Points $(i)" for i in eachindex(point1s) - n_influenced = influenced_particles[i] - - plane_points = [point1s[i], point2s[i]] - - plane_size = plane_points[2] - plane_points[1] - flow_directions = [ - normalize([-plane_size[2], plane_size[1]]), - -normalize([-plane_size[2], plane_size[1]]), + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset "Points $(i)" for i in eachindex(plane_points_1) + n_influenced = influenced_particles[i] + + plane_points = [plane_points_1[i], plane_points_2[i]] + + plane_size = plane_points[2] - plane_points[1] + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset "Flow Direction $(j)" for j in eachindex(flow_directions) + flow_direction = flow_directions[j] + inflow = InFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + outflow = OutFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, ] - @testset "Flow Direction $(j)" for j in eachindex(flow_directions) - flow_direction = flow_directions[j] - - inlet_system = OpenBoundarySPHSystem(plane_points, boundary_zone, - sound_speed; flow_direction, - particle_spacing, open_boundary_layers, - density, reference_velocity, - reference_pressure, reference_density) + @testset "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + boundary_system = OpenBoundarySPHSystem(boundary_zone, sound_speed; + reference_velocity, + reference_pressure, + reference_density) sign_ = (boundary_zone isa InFlow) ? 1 : -1 fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, @@ -52,13 +59,13 @@ density_calculator=ContinuityDensity(), smoothing_length, sound_speed) - semi = Semidiscretization(fluid_system, inlet_system) + semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, (0.0, 5.0)) v0_ode, u0_ode = ode.u0.x - v = TrixiParticles.wrap_v(v0_ode, inlet_system, semi) - u = TrixiParticles.wrap_u(u0_ode, inlet_system, semi) + v = TrixiParticles.wrap_v(v0_ode, boundary_system, semi) + u = TrixiParticles.wrap_u(u0_ode, boundary_system, semi) # ==== Characteristic Variables # `J1 = -sound_speed^2 * (rho - rho_ref) + (p - p_ref)` @@ -82,9 +89,9 @@ # First evaluation # Particles not influenced by the fluid have zero values t1 = 2.0 - TrixiParticles.evaluate_characteristics!(inlet_system, + TrixiParticles.evaluate_characteristics!(boundary_system, v, u, v0_ode, u0_ode, semi, t1) - evaluated_vars1 = inlet_system.characteristics + evaluated_vars1 = boundary_system.characteristics if boundary_zone isa InFlow @test all(isapprox.(evaluated_vars1[1, :], 0.0)) @@ -102,9 +109,9 @@ # Second evaluation # Particles not influenced by the fluid have previous values t2 = 3.0 - TrixiParticles.evaluate_characteristics!(inlet_system, + TrixiParticles.evaluate_characteristics!(boundary_system, v, u, v0_ode, u0_ode, semi, t2) - evaluated_vars2 = inlet_system.characteristics + evaluated_vars2 = boundary_system.characteristics if boundary_zone isa InFlow @test all(isapprox.(evaluated_vars2[1, :], 0.0)) diff --git a/test/schemes/boundary/open_boundary/open_boundary.jl b/test/schemes/boundary/open_boundary/open_boundary.jl new file mode 100644 index 000000000..5fe09c703 --- /dev/null +++ b/test/schemes/boundary/open_boundary/open_boundary.jl @@ -0,0 +1,2 @@ +include("characteristic_variables.jl") +include("boundary_zone.jl") diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index c25bbba21..5cf7cd5bd 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -1,5 +1,5 @@ include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") include("boundary/dummy_particles/dummy_particles.jl") include("boundary/monaghan_kajtar/monaghan_kajtar.jl") -include("boundary/open_boundary/characteristic_variables.jl") +include("boundary/open_boundary/open_boundary.jl") include("fluid/fluid.jl") diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index d8b48d620..d1abfdb40 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -1,212 +1,63 @@ @testset verbose=true "OpenBoundarySPHSystem" begin - @testset verbose=true "Boundary Zone 2D" begin - boundary_zones = [InFlow(), OutFlow()] - particle_spacing = 0.2 - open_boundary_layers = 4 - - point1s = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] - point2s = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] - - @testset "$boundary_zone" for boundary_zone in boundary_zones - @testset "Points $(i)" for i in eachindex(point1s) - plane_size = point2s[i] - point1s[i] - flow_directions = [ - normalize([-plane_size[2], plane_size[1]]), - -normalize([-plane_size[2], plane_size[1]]), - ] - - plane_points = [point1s[i], point2s[i]] - - @testset "Flow Direction $(j)" for j in eachindex(flow_directions) - system = OpenBoundarySPHSystem(plane_points, boundary_zone, 1.0; - flow_direction=flow_directions[j], - particle_spacing, - open_boundary_layers, - density=1.0) - zone_width = open_boundary_layers * - system.initial_condition.particle_spacing - sign_ = (boundary_zone isa InFlow) ? -1 : 1 - - @test point1s[i] == system.zone_origin - @test point2s[i] - system.zone_origin == system.spanning_set[2] - @test sign_ * flow_directions[j] ≈ normalize(system.spanning_set[1]) - @test zone_width ≈ norm(system.spanning_set[1]) - end - end - end - end - - @testset verbose=true "Boundary Zone 3D" begin - boundary_zones = [InFlow(), OutFlow()] - particle_spacing = 0.05 - open_boundary_layers = 4 - - point1s = [ - [0.0, 0.0, 0.0], - [0.3113730847835541, 0.19079485535621643, -0.440864622592926], - ] - point2s = [ - [1.0, 0.0, 0.0], - [-0.10468611121177673, 0.252103328704834, -0.44965094327926636], - ] - point3s = [ - [0.0, 1.0, 0.0], - [0.3113730847835541, 0.25057315826416016, -0.02374829351902008], - ] - - @testset "$boundary_zone" for boundary_zone in boundary_zones - @testset "Points $(i)" for i in eachindex(point1s) - edge1 = point2s[i] - point1s[i] - edge2 = point3s[i] - point1s[i] - - flow_directions = [ - normalize(cross(edge1, edge2)), - -normalize(cross(edge1, edge2)), - ] - - plane_points = [point1s[i], point2s[i], point3s[i]] - - @testset "Flow Direction $(j)" for j in eachindex(flow_directions) - system = OpenBoundarySPHSystem(plane_points, boundary_zone, 1.0; - flow_direction=flow_directions[j], - particle_spacing, - open_boundary_layers, - density=1.0) - zone_width = open_boundary_layers * - system.initial_condition.particle_spacing - sign_ = (boundary_zone isa InFlow) ? -1 : 1 - - @test point1s[i] == system.zone_origin - @test point2s[i] - system.zone_origin == system.spanning_set[2] - @test point3s[i] - system.zone_origin == system.spanning_set[3] - @test sign_ * flow_directions[j] ≈ normalize(system.spanning_set[1]) - @test zone_width ≈ norm(system.spanning_set[1]) - end - end - end - end - - @testset verbose=true "Particle In Boundary Zone 2D" begin - plane_points = [[-0.2, -0.5], [0.3, 0.5]] - plane_size = plane_points[2] - plane_points[1] - - flow_direction = normalize([-plane_size[2], plane_size[1]]) - system = OpenBoundarySPHSystem(plane_points, InFlow(), 1.0; flow_direction, - density=1.0, - particle_spacing=0.1, open_boundary_layers=4) - - query_points = Dict( - "Behind" => ([-1.0, -1.0], false), - "Before" => ([2.0, 2.0], false), - "On Point 1" => (plane_points[1], true), - "On Point 2" => (plane_points[2], true), - "On Point 3" => (system.spanning_set[1] + system.zone_origin, true), - "Closely On Point 1" => (plane_points[1] .+ eps() * flow_direction, false)) - - @testset "$key" for key in keys(query_points) - (particle_position, evaluation) = query_points[key] - - @test evaluation == - TrixiParticles.within_boundary_zone(particle_position, system) - end - end - - @testset verbose=true "Particle In Boundary Zone 3D" begin - point1 = [-0.2, -0.5, 0.0] - point2 = [0.3, 0.5, 0.0] - point3 = [0.13111173850909402, -0.665555869254547, 0.0] - - flow_direction = normalize(cross(point2 - point1, point3 - point1)) - system = OpenBoundarySPHSystem([point1, point2, point3], InFlow(), 1.0; - flow_direction, - density=1.0, particle_spacing=0.1, - open_boundary_layers=4) - - query_points = Dict( - "Behind" => ([-1.0, -1.0, 1.2], false), - "Before" => ([2.0, 2.0, -1.2], false), - "On Point 1" => (point1, true), - "On Point 2" => (point2, true), - "On Point 3" => (point3, true), - "On Point 4" => (system.spanning_set[1] + system.zone_origin, true), - "Closely On Point 1" => (point1 .+ eps() * flow_direction, false)) - - @testset "$key" for key in keys(query_points) - (particle_position, evaluation) = query_points[key] - - @test evaluation == - TrixiParticles.within_boundary_zone(particle_position, system) - end - end - @testset verbose=true "Illegal Inputs" begin - no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] - flow_direction = [0.0, 0.0, 1.0] - error_str = "the provided points do not span a rectangular plane" - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(no_rectangular_plane, - InFlow(), 1.0; - flow_direction, - particle_spacing=0.1, - open_boundary_layers=2, - density=1.0) - - rectangular_plane = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] - flow_direction = [0.0, 1.0, 0.0] - error_str = "flow direction and normal vector of " * - "InFlow-plane do not correspond" - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(rectangular_plane, - InFlow(), 1.0; - flow_direction, - particle_spacing=0.1, - open_boundary_layers=2, - density=1.0) - plane = ([0.0, 0.0], [0.0, 1.0]) flow_direction = (1.0, 0.0) + + inflow = InFlow(; plane, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=2) + error_str = "`reference_velocity` must be either a function mapping " * "each particle's coordinates and time to its velocity or a " * "vector of length 2 for a 2D problem" reference_velocity = 1.0 - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(plane, InFlow(), 1.0; - flow_direction, - particle_spacing=0.1, - reference_velocity, - open_boundary_layers=2, - density=1.0) + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow, 1.0; + reference_velocity) error_str = "`reference_pressure` must be either a function mapping " * "each particle's coordinates and time to its pressure or a scalar" reference_pressure = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(plane, InFlow(), 1.0; - flow_direction, - particle_spacing=0.1, - reference_pressure, - open_boundary_layers=2, - density=1.0) + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow, 1.0; + reference_pressure) error_str = "`reference_density` must be either a function mapping " * "each particle's coordinates and time to its density or a scalar" reference_density = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(plane, InFlow(), 1.0; - flow_direction, - particle_spacing=0.1, - reference_density, - open_boundary_layers=2, - density=1.0) + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow, 1.0; + reference_density) end @testset "show" begin - system = OpenBoundarySPHSystem(([0.0, 0.0], [0.0, 1.0]), InFlow(), 1.0; - flow_direction=(1.0, 0.0), density=1.0, - particle_spacing=0.05, open_boundary_layers=4) + inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(inflow, 1.0) + + show_compact = "OpenBoundarySPHSystem{2}(InFlow) with 80 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 80 │ + │ boundary: ……………………………………………………… InFlow │ + │ flow direction: ……………………………………… [1.0, 0.0] │ + │ prescribed velocity: ………………………… constant_vector │ + │ prescribed pressure: ………………………… constant_scalar │ + │ prescribed density: …………………………… constant_scalar │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + + outflow = OutFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(outflow, 1.0) - show_compact = "OpenBoundarySPHSystem{2}(InFlow()) with 80 particles" + show_compact = "OpenBoundarySPHSystem{2}(OutFlow) with 80 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ │ OpenBoundarySPHSystem{2} │ │ ════════════════════════ │ │ #particles: ………………………………………………… 80 │ - │ boundary: ……………………………………………………… InFlow() │ + │ boundary: ……………………………………………………… OutFlow │ │ flow direction: ……………………………………… [1.0, 0.0] │ │ prescribed velocity: ………………………… constant_vector │ │ prescribed pressure: ………………………… constant_scalar │