diff --git a/NEWS.md b/NEWS.md index d4e2ccdd7d..e88701f1a6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,9 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### API Changes +- API for `OpenBoundarySPHSystem` and `BoundaryZone` changed. + It is now possible to pass multiple `BoundaryZone`s to a single `OpenBoundarySPHSystem`. + Reference values are now assigned individually to each `BoundaryZone`. (#866) - Combined transport velocity formulation (TVF) and particle shifting technique (PST) into one unified framework. The keyword argument `transport_velocity` now changed to `shifting_technique`. diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 85e3461ccc..c206f76484 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -12,7 +12,7 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -particle_spacing = 0.05 +particle_spacing = 0.02 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled boundary_layers = 4 @@ -21,7 +21,7 @@ boundary_layers = 4 # fully sampled. # Note: Due to the dynamics at the inlets and outlets of open boundaries, # it is recommended to use `open_boundary_layers > boundary_layers` -open_boundary_layers = 8 +open_boundary_layers = 6 # ========================================================================================== # ==== Experiment Setup @@ -30,65 +30,71 @@ tspan = (0.0, 2.0) # Boundary geometry and initial fluid particle positions domain_size = (1.0, 0.4) -flow_direction = [1.0, 0.0] reynolds_number = 100 -const prescribed_velocity = 2.0 +const prescribed_velocity = (1.0, 0.0) +flow_direction = [1.0, 0.0] -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2]) +open_boundary_size = (particle_spacing * open_boundary_layers, domain_size[2]) fluid_density = 1000.0 -# For this particular example, it is necessary to have a background pressure. -# Otherwise the suction at the outflow is to big and the simulation becomes unstable. -pressure = 1000.0 - -sound_speed = 20 * prescribed_velocity - -state_equation = nothing +sound_speed = 10 * maximum(abs.(prescribed_velocity)) -pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, - pressure=pressure, n_layers=boundary_layers, +pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density, + n_layers=boundary_layers, velocity=prescribed_velocity, faces=(false, false, true, true)) -# Shift pipe walls in negative x-direction for the inflow -pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0) +inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true)) + +min_coords_outlet = (pipe.fluid_size[1], 0.0) +outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true)) 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 -wcsph = false +wcsph = true smoothing_length = 1.5 * particle_spacing smoothing_kernel = WendlandC2Kernel{NDIMS}() fluid_density_calculator = ContinuityDensity() -kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number +kinematic_viscosity = maximum(prescribed_velocity) * domain_size[2] / reynolds_number 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, - exponent=1, background_pressure=pressure) - alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) - viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + exponent=1) + density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, state_equation, smoothing_kernel, + density_diffusion=density_diffusion, smoothing_length, viscosity=viscosity, shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) +else + # Alternatively the EDAC scheme can be used + state_equation = nothing + + 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) end # ========================================================================================== @@ -98,63 +104,61 @@ function velocity_function2d(pos, t) # Use this for a time-dependent inflow velocity # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) - return SVector(prescribed_velocity, 0.0) + return SVector(prescribed_velocity) end -open_boundary_model = BoundaryModelLastiwka() +open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()) +reference_velocity_in = velocity_function2d +reference_pressure_in = nothing +reference_density_in = nothing boundary_type_in = InFlow() plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers, density=fluid_density, particle_spacing, - boundary_type=boundary_type_in) - -reference_velocity_in = velocity_function2d -reference_pressure_in = pressure -reference_density_in = fluid_density -open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_in, - reference_pressure=reference_pressure_in, - reference_velocity=reference_velocity_in) - + reference_density=reference_density_in, + reference_pressure=reference_pressure_in, + reference_velocity=reference_velocity_in, + initial_condition=inlet.fluid, boundary_type=boundary_type_in) + +reference_velocity_out = nothing +reference_pressure_out = nothing +reference_density_out = nothing boundary_type_out = OutFlow() -plane_out = ([domain_size[1], 0.0], [domain_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, - boundary_type=boundary_type_out) - -reference_velocity_out = velocity_function2d -reference_pressure_out = pressure -reference_density_out = fluid_density -open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_out, - reference_pressure=reference_pressure_out, - reference_velocity=reference_velocity_out) + reference_density=reference_density_out, + reference_pressure=reference_pressure_out, + reference_velocity=reference_velocity_out, + initial_condition=outlet.fluid, boundary_type=boundary_type_out) + +open_boundary = OpenBoundarySPHSystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles) + # ========================================================================================== # ==== Boundary -viscosity_boundary = ViscosityAdami(nu=1e-4) -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, +wall = union(pipe.boundary, inlet.boundary, outlet.boundary) +viscosity_boundary = viscosity +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, AdamiPressureExtrapolation(), state_equation=state_equation, viscosity=viscosity_boundary, smoothing_kernel, smoothing_length) -boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) +boundary_system = BoundarySPHSystem(wall, boundary_model) # ========================================================================================== # ==== Simulation -min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2) -max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2) +min_corner = minimum(wall.coordinates .- particle_spacing, dims=2) +max_corner = maximum(wall.coordinates .+ particle_spacing, dims=2) nhs = GridNeighborhoodSearch{NDIMS}(; cell_list=FullGridCellList(; min_corner, max_corner), update_strategy=ParallelUpdate()) -semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, - boundary_system, neighborhood_search=nhs, +semi = Semidiscretization(fluid_system, open_boundary, boundary_system, + neighborhood_search=nhs, parallelization_backend=PolyesterBackend()) ode = semidiscretize(semi, tspan) diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl index 10035f4308..595859affd 100644 --- a/examples/fluid/pipe_flow_3d.jl +++ b/examples/fluid/pipe_flow_3d.jl @@ -22,28 +22,24 @@ open_boundary_layers = 6 # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 2.0) - -function velocity_function3d(pos, t) - # Use this for a time-dependent inflow velocity - # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) - - return SVector(prescribed_velocity, 0.0, 0.0) -end +tspan = (0.0, 0.5) domain_size = (1.0, 0.4, 0.4) - -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2], domain_size[3]) - +const prescribed_velocity = (1.0, 0.0, 0.0) flow_direction = [1.0, 0.0, 0.0] +open_boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0) +min_coords_outlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0) + # setup simulation trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - domain_size=domain_size, boundary_size=boundary_size, + domain_size=domain_size, open_boundary_size=open_boundary_size, flow_direction=flow_direction, faces=(false, false, true, true, true, true), - tspan=tspan, reference_velocity=velocity_function3d, - open_boundary_layers=open_boundary_layers, + tspan=tspan, prescribed_velocity=prescribed_velocity, + open_boundary_layers=open_boundary_layers, min_coords_inlet=min_coords_inlet, + min_coords_outlet=min_coords_outlet, plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0], [0.0, 0.0, domain_size[3]]), plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0], diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index b15c8873aa..0aff00a73a 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -79,7 +79,8 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr DensityDiffusionAntuono export tensile_instability_control export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni, + PressureMirroring, PressureZeroing, BoundaryModelCharacteristicsLastiwka, + BoundaryModelMirroringTafuni, BernoulliPressureExtrapolation export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring export HertzContactModel, LinearContactModel diff --git a/src/general/buffer.jl b/src/general/buffer.jl index 3e1c28e8ca..214d9765bc 100644 --- a/src/general/buffer.jl +++ b/src/general/buffer.jl @@ -36,6 +36,10 @@ function allocate_buffer(initial_condition, buffer::SystemBuffer) return union(initial_condition, buffer_ic) end +# By default, there is no buffer. +# Dispatch by system type to handle systems that provide a buffer. +@inline buffer(system) = nothing + @inline update_system_buffer!(buffer::Nothing, semi) = buffer # TODO `resize` allocates. Find a non-allocating version diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 9532d812ab..f6cc8015c9 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -972,17 +972,17 @@ end function check_configuration(system::OpenBoundarySPHSystem, systems, neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch) - (; boundary_model, boundary_zone) = system + (; boundary_model, boundary_zones) = system # Store index of the fluid system. This is necessary for re-linking # in case we use Adapt.jl to create a new semidiscretization. fluid_system_index = findfirst(==(system.fluid_system), systems) system.fluid_system_index[] = fluid_system_index - if boundary_model isa BoundaryModelLastiwka && - boundary_zone isa BoundaryZone{BidirectionalFlow} - throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " * - "Please specify inflow and outflow.")) + if boundary_model isa BoundaryModelCharacteristicsLastiwka && + any(zone -> isnothing(zone.flow_direction), boundary_zones) + throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " * + "Please specify `InFlow()` and `OutFlow()`.")) end if first(PointNeighbors.requires_update(neighborhood_search)) @@ -1011,10 +1011,8 @@ function set_system_links(system::OpenBoundarySPHSystem, semi) system.pressure, system.boundary_candidates, system.fluid_candidates, - system.boundary_zone, - system.reference_velocity, - system.reference_pressure, - system.reference_density, + system.boundary_zone_indices, + system.boundary_zones, system.buffer, system.cache) end diff --git a/src/general/system.jl b/src/general/system.jl index 5b4d5aecb8..6c4aed9b0f 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -37,17 +37,18 @@ initialize!(system, semi) = system # Number of particles in the system whose positions are to be integrated (corresponds to the size of u and du) @inline n_moving_particles(system) = nparticles(system) -@inline eachparticle(system) = Base.OneTo(nparticles(system)) +@inline eachparticle(system::System) = active_particles(system) +@inline eachparticle(initial_condition) = Base.OneTo(nparticles(initial_condition)) # Wrapper for systems with `SystemBuffer` -@inline each_moving_particle(system) = each_moving_particle(system, system.buffer) +@inline each_moving_particle(system) = each_moving_particle(system, buffer(system)) @inline each_moving_particle(system, ::Nothing) = Base.OneTo(n_moving_particles(system)) -@inline active_coordinates(u, system) = active_coordinates(u, system, system.buffer) +@inline active_coordinates(u, system) = active_coordinates(u, system, buffer(system)) @inline active_coordinates(u, system, ::Nothing) = current_coordinates(u, system) -@inline active_particles(system) = active_particles(system, system.buffer) -@inline active_particles(system, ::Nothing) = eachparticle(system) +@inline active_particles(system) = active_particles(system, buffer(system)) +@inline active_particles(system, ::Nothing) = Base.OneTo(nparticles(system)) # This should not be dispatched by system type. We always expect to get a column of `A`. @propagate_inbounds function extract_svector(A, system, i) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 2d8b220f05..c17410d698 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -137,12 +137,12 @@ function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box; write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) # Store particle index - vtk["index"] = active_particles(system) + vtk["index"] = eachparticle(system) vtk["time"] = t vtk["ndims"] = ndims(system) vtk["particle_spacing"] = [particle_spacing(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if write_meta_data vtk["solver_version"] = git_hash @@ -294,20 +294,20 @@ end function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true) vtk["velocity"] = view(v, 1:ndims(system), :) vtk["mass"] = [hydrodynamic_mass(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["radius"] = [particle_radius(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] return vtk end function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["density"] = [current_density(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] # Indexing the pressure is a workaround for slicing issue (see https://github.com/JuliaSIMD/StrideArrays.jl/issues/88) vtk["pressure"] = [current_pressure(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if system.surface_normal_method !== nothing vtk["surf_normal"] = [surface_normal(system, particle) @@ -406,7 +406,7 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d n_fixed_particles = nparticles(system) - n_moving_particles(system) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["jacobian"] = [det(deformation_gradient(system, particle)) for particle in eachparticle(system)] @@ -438,19 +438,11 @@ end function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["density"] = [current_density(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["pressure"] = [current_pressure(v, system, particle) - for particle in active_particles(system)] - - if write_meta_data - vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters)) - vtk["width"] = round(system.boundary_zone.zone_width, digits=3) - vtk["velocity_function"] = type2string(system.reference_velocity) - vtk["pressure_function"] = type2string(system.reference_pressure) - vtk["density_function"] = type2string(system.reference_density) - end + for particle in eachparticle(system)] return vtk end diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index 6f73fb2c68..100d0f4963 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -218,7 +218,7 @@ end function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true) vtk["velocity"] = [advection_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if write_meta_data vtk["signed_distances"] = system.signed_distances end diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 30e589e1c9..6738552f4b 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -7,8 +7,10 @@ struct OutFlow end @doc raw""" BoundaryZone(; plane, plane_normal, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer, boundary_type=BidirectionalFlow(), - average_inflow_velocity=true) + open_boundary_layers::Integer, average_inflow_velocity=true, + boundary_type=BidirectionalFlow(), + reference_density=nothing, reference_pressure=nothing, + reference_velocity=nothing) Boundary zone for [`OpenBoundarySPHSystem`](@ref). @@ -63,48 +65,100 @@ There are three ways to specify the actual shape of the boundary zone: anisotropic buffer-particles distribution, resulting in a potential numerical instability. Averaging mitigates these effects. +- `reference_velocity`: Reference velocity is either a function mapping each particle's coordinates + and time to its velocity, or, for a constant fluid velocity, + a vector holding this velocity. +- `reference_pressure`: Reference pressure is either a function mapping each particle's coordinates + and time to its pressure, or a scalar for a constant pressure over all particles. +- `reference_density`: Reference density is either a function mapping each particle's coordinates + and time to its density, or a scalar for a constant density over all particles. + +!!! 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 ([BoundaryModelMirroringTafuni](@ref BoundaryModelMirroringTafuni)) + or evolved using the characteristic flow variables ([BoundaryModelCharacteristicsLastiwka](@ref BoundaryModelCharacteristicsLastiwka)). # Examples -```julia +```jldoctest; output = false # 2D plane_points = ([0.0, 0.0], [0.0, 1.0]) -plane_normal=[1.0, 0.0] +plane_normal = [1.0, 0.0] + +# Constant reference velocity: +velocity_const = [1.0, 0.0] + +inflow_1 = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, + density=1.0, open_boundary_layers=4, boundary_type=InFlow(), + reference_velocity=velocity_const) + +# Reference velocity as a function (parabolic velocity profile): +velocity_func = (pos, t) -> SVector(4.0 * pos[2] * (1.0 - pos[2]), 0.0) -inflow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0, - open_boundary_layers=4, boundary_type=InFlow()) +inflow_2 = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, + density=1.0, open_boundary_layers=4, boundary_type=InFlow(), + reference_velocity=velocity_func) # 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] +plane_normal = [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()) +# Constant reference pressure: +pressure_const = 0.0 + +outflow_1 = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0, + open_boundary_layers=4, boundary_type=OutFlow(), + reference_pressure=pressure_const) + +# Reference pressure as a function (y-dependent profile, sinusoidal in time): +pressure_func = (pos, t) -> pos[2] * sin(2pi * t) + +outflow_2 = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0, + open_boundary_layers=4, boundary_type=OutFlow(), + reference_pressure=pressure_func) # 3D particles sampled as cylinder circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) bidirectional_flow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, - density=1.0, extrude_geometry=circle, open_boundary_layers=4) + density=1.0, boundary_type=BidirectionalFlow(), + extrude_geometry=circle, open_boundary_layers=4) + +# output +┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ +│ BoundaryZone │ +│ ════════════ │ +│ boundary type: ………………………………………… bidirectional_flow │ +│ #particles: ………………………………………………… 234 │ +│ width: ……………………………………………………………… 0.4 │ +└──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct BoundaryZone{BT, IC, S, ZO, ZW, FD, PN} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD - plane_normal :: PN - boundary_type :: BT +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 + # 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 + prescribed_velocity :: Bool end function BoundaryZone(; plane, plane_normal, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer, boundary_type=BidirectionalFlow(), - average_inflow_velocity=true) + open_boundary_layers::Integer, average_inflow_velocity=true, + boundary_type=BidirectionalFlow(), + reference_density=nothing, reference_pressure=nothing, + reference_velocity=nothing) if open_boundary_layers <= 0 throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) end @@ -112,52 +166,147 @@ 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 + 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) + if !(reference_velocity isa Function || isnothing(reference_velocity) || + (reference_velocity isa Vector && 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")) + else + if reference_velocity isa Function + test_result = reference_velocity(zeros(NDIMS), 0.0) + if length(test_result) != NDIMS + throw(ArgumentError("`velocity` function must be of dimension $NDIMS")) + end + end + # We need this dummy for type stability reasons + velocity_dummy = SVector(ntuple(dim -> convert(ELTYPE, Inf), NDIMS)) + velocity_ref = wrap_reference_function(reference_velocity, velocity_dummy) + end - elseif boundary_type isa InFlow - # Unit vector pointing in downstream direction - flow_direction = plane_normal_ + if !(reference_pressure isa Function || reference_pressure isa Real || + isnothing(reference_pressure)) + throw(ArgumentError("`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "or a scalar")) + else + if reference_pressure isa Function + test_result = reference_pressure(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_pressure` function must be a scalar function")) + end + end + # We need this dummy for type stability reasons + pressure_dummy = convert(ELTYPE, Inf) + pressure_ref = wrap_reference_function(reference_pressure, pressure_dummy) + end - elseif boundary_type isa OutFlow - # Unit vector pointing in downstream direction - flow_direction = -plane_normal_ + if !(reference_density isa Function || reference_density isa Real || + isnothing(reference_density)) + throw(ArgumentError("`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "or a scalar")) + else + if reference_density isa Function + test_result = reference_density(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_density` function must be a scalar function")) + end + end + # We need this dummy for type stability reasons + density_dummy = convert(ELTYPE, Inf) + density_ref = wrap_reference_function(reference_density, density_dummy) 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) + prescribed_pressure = isnothing(reference_pressure) ? false : true + prescribed_density = isnothing(reference_density) ? false : true + prescribed_velocity = isnothing(reference_velocity) ? false : true + + reference_values = (reference_velocity=velocity_ref, reference_pressure=pressure_ref, + reference_density=density_ref) + + coordinates_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, ic.coordinates) + + if prescribed_pressure + ic.pressure .= pressure_ref.(coordinates_svector, 0) + end + if prescribed_density + ic.density .= density_ref.(coordinates_svector, 0) + ic.mass .= ic.density * ic.particle_spacing^NDIMS + end + if prescribed_velocity + ic.velocity .= stack(velocity_ref.(coordinates_svector, 0)) + end return BoundaryZone(ic, spanning_set_, zone_origin, zone_width, - flow_direction, plane_normal_, boundary_type, - average_inflow_velocity) + flow_direction, plane_normal_, reference_values, + average_inflow_velocity, prescribed_density, prescribed_pressure, + prescribed_velocity) +end + +function boundary_type_name(boundary_zone::BoundaryZone) + (; flow_direction, plane_normal) = boundary_zone + + 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) + @nospecialize boundary_zone # reduce precompilation time + + print(io, "BoundaryZone(") + print(io, ") with ", nparticles(boundary_zone.initial_condition), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone) + @nospecialize boundary_zone # reduce precompilation time + + if get(io, :compact, false) + show(io, boundary_zone) + else + summary_header(io, "BoundaryZone") + summary_line(io, "boundary type", boundary_type_name(boundary_zone)) + summary_line(io, "#particles", nparticles(boundary_zone.initial_condition)) + summary_line(io, "width", round(boundary_zone.zone_width, digits=3)) + summary_footer(io) + 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 @@ -205,7 +354,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) @@ -237,7 +386,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) return hcat(c, edge1, edge2) end -@inline function is_in_boundary_zone(boundary_zone::BoundaryZone, particle_coords) +@inline function is_in_boundary_zone(boundary_zone, particle_coords) (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin @@ -261,6 +410,27 @@ end return true end +function update_boundary_zone_indices!(system, u, boundary_zones, semi) + set_zero!(system.boundary_zone_indices) + + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + for (zone_id, boundary_zone) in enumerate(boundary_zones) + # Check if boundary particle is in the boundary zone + if is_in_boundary_zone(boundary_zone, particle_coords) + system.boundary_zone_indices[particle] = zone_id + end + end + end + + return system +end + +function current_boundary_zone(system, particle) + return system.boundary_zones[system.boundary_zone_indices[particle]] +end + function remove_outside_particles(initial_condition, spanning_set, zone_origin) (; coordinates, density, particle_spacing) = initial_condition @@ -276,3 +446,64 @@ function remove_outside_particles(initial_condition, spanning_set, zone_origin) return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), particle_spacing) end + +function wrap_reference_function(function_::Nothing, ref_dummy) + # Return a dummy value for type stability + return @inline((coords, t)->ref_dummy) +end + +function wrap_reference_function(function_::Function, ref_dummy) + # Already a function + return function_ +end + +function wrap_reference_function(constant_scalar::Number, ref_dummy) + return @inline((coords, t)->constant_scalar) +end + +function wrap_reference_function(constant_vector::AbstractVector, + ref_dummy::SVector{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} + return @inline((coords, t)->SVector{NDIMS, ELTYPE}(constant_vector)) +end + +function reference_pressure(boundary_zone, v, system, particle, pos, t) + (; prescribed_pressure) = boundary_zone + (; pressure_reference_values) = system.cache + + if prescribed_pressure + zone_id = system.boundary_zone_indices[particle] + + # `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 reference_density(boundary_zone, v, system, particle, pos, t) + (; prescribed_density) = boundary_zone + (; density_reference_values) = system.cache + + if prescribed_density + zone_id = system.boundary_zone_indices[particle] + + # `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 reference_velocity(boundary_zone, v, system, particle, pos, t) + (; prescribed_velocity) = boundary_zone + (; velocity_reference_values) = system.cache + + if prescribed_velocity + zone_id = system.boundary_zone_indices[particle] + + # `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 diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index a37a9cd605..4e25500462 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -1,5 +1,5 @@ @doc raw""" - BoundaryModelLastiwka(; extrapolate_reference_values=nothing) + BoundaryModelCharacteristicsLastiwka(; extrapolate_reference_values=nothing) Boundary model for [`OpenBoundarySPHSystem`](@ref). This model uses the characteristic variables to propagate the appropriate values @@ -19,58 +19,55 @@ For more information about the method see [description below](@ref method_of_cha Note that even without this extrapolation feature, the reference values don't need to be prescribed - they're computed from the characteristics. """ -struct BoundaryModelLastiwka{T} +struct BoundaryModelCharacteristicsLastiwka{T} extrapolate_reference_values::T - function BoundaryModelLastiwka(; extrapolate_reference_values=nothing) + function BoundaryModelCharacteristicsLastiwka(; extrapolate_reference_values=nothing) return new{typeof(extrapolate_reference_values)}(extrapolate_reference_values) end end # Called from update callback via `update_open_boundary_eachstep!` -@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, +@inline function update_boundary_quantities!(system, + boundary_model::BoundaryModelCharacteristicsLastiwka, v, u, v_ode, u_ode, semi, t) - (; density, pressure, cache, boundary_zone, - reference_velocity, reference_pressure, reference_density) = system - (; flow_direction) = boundary_zone - - fluid_system = corresponding_fluid_system(system, semi) + (; density, pressure, cache, boundary_zones, fluid_system) = system sound_speed = system_sound_speed(fluid_system) if !isnothing(boundary_model.extrapolate_reference_values) - (; prescribed_pressure, prescribed_velocity, prescribed_density) = cache v_fluid = wrap_v(v_ode, fluid_system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) @trixi_timeit timer() "extrapolate and correct values" begin extrapolate_values!(system, boundary_model.extrapolate_reference_values, - v, v_fluid, u, u_fluid, semi, t; - prescribed_pressure, prescribed_velocity, - prescribed_density) + v, v_fluid, u, u_fluid, semi) end end # Update quantities based on the characteristic variables @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; flow_direction) = boundary_zone + particle_position = current_coords(u, system, particle) J1 = cache.characteristics[1, particle] J2 = cache.characteristics[2, particle] J3 = cache.characteristics[3, particle] - rho_ref = reference_value(reference_density, density[particle], - particle_position, t) + rho_ref = reference_density(boundary_zone, v, system, particle, + particle_position, t) + density[particle] = rho_ref + ((-J1 + (J2 + J3) / 2) / sound_speed^2) - p_ref = reference_value(reference_pressure, pressure[particle], - particle_position, t) + p_ref = reference_pressure(boundary_zone, v, system, particle, particle_position, t) + pressure[particle] = p_ref + (J2 + J3) / 2 - v_current = current_velocity(v, system, particle) - v_ref = reference_value(reference_velocity, v_current, - particle_position, t) - rho = density[particle] + v_ref = reference_velocity(boundary_zone, v, system, particle, particle_position, t) + + rho = current_density(v, system, particle) v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction for dim in 1:ndims(system) @@ -78,19 +75,21 @@ end end end - if boundary_zone.average_inflow_velocity - # Even if the velocity is prescribed, this boundary model computes the velocity for each particle individually. - # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, - # resulting in a potential numerical instability. Averaging mitigates these effects. - average_velocity!(v, u, system, boundary_model, boundary_zone, semi) + for boundary_zone in boundary_zones + if boundary_zone.average_inflow_velocity + # Even if the velocity is prescribed, this boundary model computes the velocity for each particle individually. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_model, boundary_zone, semi) + end end return system end # Called from semidiscretization -function update_boundary_model!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, - semi, t) +function update_boundary_model!(system, ::BoundaryModelCharacteristicsLastiwka, + v, u, v_ode, u_ode, semi, t) @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end @@ -103,9 +102,8 @@ end # J2: Propagates downstream to the local flow # J3: Propagates upstream to the local flow function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) - (; volume, cache, boundary_zone) = system + (; volume, cache, fluid_system, density, pressure) = system (; characteristics, previous_characteristics) = cache - fluid_system = corresponding_fluid_system(system, semi) @threaded semi for particle in eachparticle(system) previous_characteristics[1, particle] = characteristics[1, particle] @@ -117,14 +115,56 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) set_zero!(volume) # Evaluate the characteristic variables with the fluid system - evaluate_characteristics!(system, fluid_system, v, u, v_ode, u_ode, semi, t) + v_fluid = wrap_v(v_ode, fluid_system, semi) + u_fluid = wrap_u(u_ode, fluid_system, semi) + + system_coords = current_coordinates(u, system) + fluid_coords = current_coordinates(u_fluid, fluid_system) + sound_speed = system_sound_speed(fluid_system) + + # Loop over all fluid neighbors within the kernel cutoff + foreach_point_neighbor(system, fluid_system, system_coords, fluid_coords, semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + boundary_zone = current_boundary_zone(system, particle) + (; flow_direction) = boundary_zone + + neighbor_position = current_coords(u_fluid, fluid_system, neighbor) + + # Determine current and prescribed quantities + rho_b = current_density(v_fluid, fluid_system, neighbor) + + rho_ref = reference_density(boundary_zone, v, system, particle, + neighbor_position, t) + + p_b = current_pressure(v_fluid, fluid_system, neighbor) + + p_ref = reference_pressure(boundary_zone, v, system, particle, neighbor_position, t) + + v_b = current_velocity(v_fluid, fluid_system, neighbor) + + v_neighbor_ref = reference_velocity(boundary_zone, v, system, particle, + neighbor_position, t) + + # Determine characteristic variables + density_term = -sound_speed^2 * (rho_b - rho_ref) + pressure_term = p_b - p_ref + velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) + + kernel_ = smoothing_kernel(fluid_system, distance, particle) + + characteristics[1, particle] += (density_term + pressure_term) * kernel_ + characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ + characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ + + volume[particle] += kernel_ + end # Only some of the in-/outlet particles are in the influence of the fluid particles. # Thus, we compute the characteristics for the particles that are outside the influence # of fluid particles by using the average of the values of the previous time step. # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119 @threaded semi for particle in each_moving_particle(system) - # Particle is outside of the influence of fluid particles if isapprox(volume[particle], 0) @@ -161,93 +201,43 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) characteristics[2, particle] /= volume[particle] characteristics[3, particle] /= volume[particle] end - prescribe_conditions!(characteristics, particle, boundary_zone) - end - - return system -end - -function evaluate_characteristics!(system, neighbor_system::FluidSystem, - v, u, v_ode, u_ode, semi, t) - (; 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) - u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - sound_speed = system_sound_speed(neighbor_system) - - # Loop over all fluid neighbors within the kernel cutoff - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - neighbor_position = current_coords(u_neighbor_system, neighbor_system, neighbor) - - # Determine current and prescribed quantities - rho_b = current_density(v_neighbor_system, neighbor_system, neighbor) - rho_ref = reference_value(reference_density, density[particle], - neighbor_position, t) - - p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor) - p_ref = reference_value(reference_pressure, pressure[particle], - neighbor_position, t) - - v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - v_particle = current_velocity(v, system, particle) - v_neighbor_ref = reference_value(reference_velocity, v_particle, - neighbor_position, t) - - # Determine characteristic variables - density_term = -sound_speed^2 * (rho_b - rho_ref) - pressure_term = p_b - p_ref - velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) - kernel_ = smoothing_kernel(neighbor_system, distance, particle) + boundary_zone = current_boundary_zone(system, particle) + (; flow_direction, plane_normal) = boundary_zone - characteristics[1, particle] += (density_term + pressure_term) * kernel_ - characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ - characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ + # Outflow + if signbit(dot(flow_direction, plane_normal)) + # 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)) - volume[particle] += kernel_ + else # Inflow + # Allow only J3 to propagate upstream to the boundary + characteristics[1, particle] = zero(eltype(characteristics)) + characteristics[2, particle] = zero(eltype(characteristics)) + end end return system end -@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)) +function average_velocity!(v, u, system, ::BoundaryModelCharacteristicsLastiwka, + boundary_zone, semi) + (; flow_direction, plane_normal) = boundary_zone - return characteristics -end - -@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)) - - return characteristics -end + # This is an outflow. Only apply averaging at the inflow. + signbit(dot(flow_direction, plane_normal)) && return v -function average_velocity!(v, u, system, ::BoundaryModelLastiwka, boundary_zone, semi) - # Only apply averaging at the inflow - return v -end - -function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone{InFlow}, - semi) + particles_in_zone = findall(particle -> boundary_zone == + current_boundary_zone(system, particle), + each_moving_particle(system)) # Division inside the `sum` closure to maintain GPU compatibility - avg_velocity = sum(each_moving_particle(system)) do particle - return current_velocity(v, system, particle) / system.buffer.active_particle_count[] + avg_velocity = sum(particles_in_zone) do particle + return current_velocity(v, system, particle) / length(particles_in_zone) end - @threaded semi for particle in each_moving_particle(system) + @threaded semi for particle in particles_in_zone # Set the velocity of the ghost node to the average velocity of the fluid domain for dim in eachindex(avg_velocity) @inbounds v[dim, particle] = avg_velocity[dim] diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 45a98ad7c2..e1f33fed7e 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -47,7 +47,7 @@ The interpolated values at the ghost nodes are then assigned to the correspondin struct ZerothOrderMirroring end @doc raw""" - BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) + BoundaryModelMirroringTafuni(; mirror_method=FirstOrderMirroring()) Boundary model for the `OpenBoundarySPHSystem`. This model implements the method of [Tafuni et al. (2018)](@cite Tafuni2018) to extrapolate the properties from the fluid domain @@ -61,62 +61,56 @@ We provide three different mirroring methods: - [`FirstOrderMirroring`](@ref): Uses a first order correction based on the gradient of the interpolated values . - [`SimpleMirroring`](@ref): Similar to the first order mirroring, but does not use the gradient of the interpolated values. """ -struct BoundaryModelTafuni{MM} +struct BoundaryModelMirroringTafuni{MM} mirror_method::MM end -function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) - return BoundaryModelTafuni(mirror_method) +function BoundaryModelMirroringTafuni(; mirror_method=FirstOrderMirroring()) + return BoundaryModelMirroringTafuni(mirror_method) end -function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, +function update_boundary_quantities!(system, boundary_model::BoundaryModelMirroringTafuni, v, u, v_ode, u_ode, semi, t) - (; reference_pressure, reference_density, reference_velocity, boundary_zone, - pressure, density, cache) = system - (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache + (; boundary_zones, pressure, density, fluid_system, cache) = system @trixi_timeit timer() "extrapolate and correct values" begin - fluid_system = corresponding_fluid_system(system, semi) - v_fluid = wrap_v(v_ode, fluid_system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) - extrapolate_values!(system, boundary_model.mirror_method, v, v_fluid, - u, u_fluid, semi, t; prescribed_pressure, - prescribed_density, prescribed_velocity) - end - - if !prescribed_velocity && boundary_zone.average_inflow_velocity - # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. - # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, - # resulting in a potential numerical instability. Averaging mitigates these effects. - average_velocity!(v, u, system, boundary_zone, semi) + extrapolate_values!(system, boundary_model.mirror_method, + v, v_fluid, u, u_fluid, semi) end - if prescribed_pressure - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) + for boundary_zone in boundary_zones + (; average_inflow_velocity, prescribed_velocity) = boundary_zone - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) + if !prescribed_velocity && average_inflow_velocity + # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_zone, semi) end end - if prescribed_density - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) + @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone + + particle_coords = current_coords(u, system, particle) - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) + if prescribed_pressure + pressure[particle] = reference_pressure(boundary_zone, v, system, particle, + particle_coords, t) end - end - if prescribed_velocity - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) - v_particle = current_velocity(v, system, particle) + if prescribed_density + density[particle] = reference_density(boundary_zone, v, system, particle, + particle_coords, t) + end - v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) + if prescribed_velocity + v_ref = reference_velocity(boundary_zone, v, system, particle, + particle_coords, t) for dim in eachindex(v_ref) @inbounds v[dim, particle] = v_ref[dim] @@ -125,16 +119,15 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni end end -update_boundary_model!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system +function update_boundary_model!(system, ::BoundaryModelMirroringTafuni, v, u, v_ode, u_ode, + semi, t) + return system +end function extrapolate_values!(system, mirror_method::Union{FirstOrderMirroring, SimpleMirroring}, - v_open_boundary, v_fluid, u_open_boundary, u_fluid, - semi, t; prescribed_density=false, - prescribed_pressure=false, prescribed_velocity=false) - (; pressure, density, boundary_zone) = system - - fluid_system = corresponding_fluid_system(system, semi) + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi) + (; pressure, density, fluid_system) = system # Static indices to avoid allocations two_to_end = SVector{ndims(system)}(2:(ndims(system) + 1)) @@ -149,6 +142,9 @@ function extrapolate_values!(system, # We can do this because we require the neighborhood search to support querying neighbors # of arbitrary positions (see `PointNeighbors.requires_update`). @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone + particle_coords = current_coords(u_open_boundary, system, particle) ghost_node_position = mirror_position(particle_coords, boundary_zone) @@ -280,12 +276,8 @@ function extrapolate_values!(system, end function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, - v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi, t; - prescribed_density=false, prescribed_pressure=false, - prescribed_velocity=false) - (; pressure, density, boundary_zone) = system - - fluid_system = corresponding_fluid_system(system, semi) + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi) + (; pressure, density, fluid_system) = system # Use the fluid-fluid nhs, since the boundary particles are mirrored into the fluid domain nhs = get_neighborhood_search(fluid_system, fluid_system, semi) @@ -297,6 +289,9 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, # We can do this because we require the neighborhood search to support querying neighbors # of arbitrary positions (see `PointNeighbors.requires_update`). @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_pressure, prescribed_density, prescribed_velocity) = boundary_zone + particle_coords = current_coords(u_open_boundary, system, particle) ghost_node_position = mirror_position(particle_coords, boundary_zone) @@ -486,10 +481,15 @@ function mirror_position(particle_coords, boundary_zone) return particle_coords - 2 * dist * boundary_zone.plane_normal end -average_velocity!(v, u, system, boundary_zone, semi) = v +# Only for inflow boundary zones +function average_velocity!(v, u, system, boundary_zone, semi) + (; plane_normal, zone_origin, initial_condition, flow_direction) = boundary_zone -function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, semi) - (; plane_normal, zone_origin, initial_condition) = boundary_zone + # Bidirectional flow + isnothing(flow_direction) && return v + + # Outflow + signbit(dot(flow_direction, plane_normal)) && return v # We only use the extrapolated velocity in the vicinity of the transition region. # Otherwise, if the boundary zone is too large, averaging would be excessively influenced @@ -497,15 +497,20 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se max_dist = initial_condition.particle_spacing * 110 / 100 candidates = findall(x -> dot(x - zone_origin, -plane_normal) <= max_dist, - reinterpret(reshape, SVector{ndims(system), eltype(u)}, - active_coordinates(u, system))) + reinterpret(reshape, SVector{ndims(system), eltype(u)}, u)) + + particles_in_zone = findall(particle -> boundary_zone == + current_boundary_zone(system, particle), + each_moving_particle(system)) + + intersect!(candidates, particles_in_zone) # Division inside the `sum` closure to maintain GPU compatibility avg_velocity = sum(candidates) do particle return current_velocity(v, system, particle) / length(candidates) end - @threaded semi for particle in each_moving_particle(system) + @threaded semi for particle in particles_in_zone # Set the velocity of the ghost node to the average velocity of the fluid domain for dim in eachindex(avg_velocity) @inbounds v[dim, particle] = avg_velocity[dim] @@ -515,10 +520,16 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se return v end -project_velocity_on_plane_normal!(v, system, particle, boundary_zone) = v +# Only for inflow boundary zones +function project_velocity_on_plane_normal!(v, system, particle, boundary_zone) + (; plane_normal, flow_direction) = boundary_zone + + # Bidirectional flow + isnothing(flow_direction) && return v + + # Outflow + signbit(dot(flow_direction, plane_normal)) && return v -function project_velocity_on_plane_normal!(v, system, particle, - boundary_zone::BoundaryZone{InFlow}) # Project `vel` on the normal direction of the boundary zone # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: # "Because flow from the inlet interface occurs perpendicular to the boundary, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4cadc62119..a7fe53308d 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,10 +1,7 @@ @doc raw""" OpenBoundarySPHSystem(boundary_zone::BoundaryZone; fluid_system::FluidSystem, buffer_size::Integer, - boundary_model, - reference_velocity=nothing, - reference_pressure=nothing, - reference_density=nothing) + boundary_model) Open boundary system for in- and outflow particles. @@ -15,184 +12,136 @@ Open boundary system for in- and outflow particles. - `fluid_system`: The corresponding fluid system - `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models)) - `buffer_size`: 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 - the velocity of particle ``i`` or, for a constant fluid velocity, - a vector holding this velocity. -- `reference_pressure`: Reference pressure is either a function mapping each particle's coordinates - and time to its pressure, a vector holding the pressure of each particle, - or a scalar for a constant pressure over all particles. -- `reference_density`: Reference density is either a function mapping each particle's coordinates - and time to its density, a vector holding the density of each particle, - or a scalar for a constant density over all particles. - -!!! 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 ([BoundaryModelTafuni](@ref BoundaryModelTafuni)) - or evolved using the characteristic flow variables ([BoundaryModelLastiwka](@ref BoundaryModelLastiwka)). !!! warning "Experimental Implementation" - This is an experimental feature and may change in future releases. - It is GPU-compatible (e.g., with CUDA.jl and AMDGPU.jl), but currently **not** supported with Metal.jl. + This is an experimental feature and may change in any future releases. """ -struct OpenBoundarySPHSystem{BM, ELTYPE, NDIMS, IC, FS, FSI, ARRAY1D, BC, FC, BZ, RV, - RP, RD, B, C} <: System{NDIMS} - boundary_model :: BM - initial_condition :: IC - fluid_system :: FS - fluid_system_index :: FSI - smoothing_length :: ELTYPE - mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] - density :: ARRAY1D # Array{ELTYPE, 1}: [particle] - volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] - pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] - boundary_candidates :: BC # Array{UInt32, 1}: [particle] - fluid_candidates :: FC # Array{UInt32, 1}: [particle] - boundary_zone :: BZ - reference_velocity :: RV - reference_pressure :: RP - reference_density :: RD - buffer :: B - cache :: C +struct OpenBoundarySPHSystem{BM, ELTYPE, NDIMS, IC, FS, FSI, ARRAY1D, BC, FC, BZI, BZ, + B, C} <: System{NDIMS} + boundary_model :: BM + initial_condition :: IC + fluid_system :: FS + fluid_system_index :: FSI + smoothing_length :: ELTYPE + mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] + density :: ARRAY1D # Array{ELTYPE, 1}: [particle] + volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] + pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] + boundary_candidates :: BC # Array{Bool, 1}: [particle] + fluid_candidates :: FC # Array{Bool, 1}: [particle] + boundary_zone_indices :: BZI # Array{UInt8, 1}: [particle] + boundary_zones :: BZ + buffer :: B + cache :: C end function OpenBoundarySPHSystem(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, fluid_candidates, - boundary_zone, reference_velocity, - reference_pressure, reference_density, buffer, cache) + boundary_zone_indices, boundary_zone, buffer, cache) OpenBoundarySPHSystem{typeof(boundary_model), eltype(mass), ndims(initial_condition), typeof(initial_condition), typeof(fluid_system), typeof(fluid_system_index), typeof(mass), typeof(boundary_candidates), typeof(fluid_candidates), - typeof(boundary_zone), typeof(reference_velocity), - typeof(reference_pressure), typeof(reference_density), + typeof(boundary_zone_indices), typeof(boundary_zone), typeof(buffer), typeof(cache)}(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, - fluid_candidates, boundary_zone, - reference_velocity, reference_pressure, - reference_density, buffer, cache) + fluid_candidates, boundary_zone_indices, + boundary_zone, buffer, cache) end -function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; - fluid_system::FluidSystem, - buffer_size::Integer, boundary_model, - reference_velocity=nothing, - reference_pressure=nothing, - reference_density=nothing) - (; initial_condition) = boundary_zone +function OpenBoundarySPHSystem(boundary_zones::Union{BoundaryZone, Nothing}...; + fluid_system::FluidSystem, buffer_size::Integer, + boundary_model) + boundary_zones_ = filter(bz -> !isnothing(bz), boundary_zones) + reference_values_ = map(bz -> bz.reference_values, boundary_zones_) - buffer = SystemBuffer(nparticles(initial_condition), buffer_size) + initial_conditions = union((bz.initial_condition for bz in boundary_zones)...) - initial_condition = allocate_buffer(initial_condition, buffer) + buffer = SystemBuffer(nparticles(initial_conditions), buffer_size) - NDIMS = ndims(initial_condition) + initial_conditions = allocate_buffer(initial_conditions, buffer) - pressure = copy(initial_condition.pressure) - mass = copy(initial_condition.mass) - density = copy(initial_condition.density) - volume = similar(initial_condition.density) + pressure = copy(initial_conditions.pressure) + mass = copy(initial_conditions.mass) + density = copy(initial_conditions.density) + volume = similar(initial_conditions.density) - if !(reference_velocity isa Function || isnothing(reference_velocity) || - (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) - throw(ArgumentError("`reference_velocity` must be either a function mapping " * - "each particle's coordinates and time to its velocity, " * - "an array where the ``i``-th column holds the velocity of particle ``i`` " * - "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) - else - if reference_velocity isa Function - test_result = reference_velocity(zeros(NDIMS), 0.0) - if length(test_result) != NDIMS - throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) - end - end - reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) - end - - if !(reference_pressure isa Function || reference_pressure isa Real || - isnothing(reference_pressure)) - throw(ArgumentError("`reference_pressure` must be either a function mapping " * - "each particle's coordinates and time to its pressure, " * - "a vector holding the pressure of each particle, or a scalar")) - else - if reference_pressure isa Function - test_result = reference_pressure(zeros(NDIMS), 0.0) - if length(test_result) != 1 - throw(ArgumentError("`reference_pressure` function must be a scalar function")) - end - end - reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) - end - - if !(reference_density isa Function || reference_density isa Real || - isnothing(reference_density)) - throw(ArgumentError("`reference_density` must be either a function mapping " * - "each particle's coordinates and time to its density, " * - "a vector holding the density of each particle, or a scalar")) - else - if reference_density isa Function - test_result = reference_density(zeros(NDIMS), 0.0) - if length(test_result) != 1 - throw(ArgumentError("`reference_density` function must be a scalar function")) - end - end - reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) - end - - cache = create_cache_open_boundary(boundary_model, initial_condition, - reference_density, reference_velocity, - reference_pressure) + cache = create_cache_open_boundary(boundary_model, initial_conditions, + reference_values_) fluid_system_index = Ref(0) smoothing_length = initial_smoothing_length(fluid_system) - boundary_candidates = fill(false, nparticles(initial_condition)) + boundary_candidates = fill(false, nparticles(initial_conditions)) fluid_candidates = fill(false, nparticles(fluid_system)) - return OpenBoundarySPHSystem(boundary_model, initial_condition, fluid_system, + boundary_zone_indices = zeros(Int, nparticles(initial_conditions)) + + # Create new `BoundaryZone`s with `reference_values` set to `nothing` for type stability. + # `reference_values` are only used as API feature to temporarily store the reference values + # in the `BoundaryZone`, but they are not used in the actual simulation. + boundary_zones_new = map(zone -> BoundaryZone(zone.initial_condition, + zone.spanning_set, + zone.zone_origin, + zone.zone_width, + zone.flow_direction, + zone.plane_normal, + nothing, + zone.average_inflow_velocity, + zone.prescribed_density, + zone.prescribed_pressure, + zone.prescribed_velocity), + boundary_zones) + + return OpenBoundarySPHSystem(boundary_model, initial_conditions, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, fluid_candidates, - boundary_zone, reference_velocity_, - reference_pressure_, reference_density_, buffer, cache) + boundary_zone_indices, boundary_zones_new, buffer, cache) end -function create_cache_open_boundary(boundary_model, initial_condition, - reference_density, reference_velocity, - reference_pressure) - ELTYPE = eltype(initial_condition) +function initialize!(system::OpenBoundarySPHSystem, semi) + (; boundary_zones) = system - prescribed_pressure = isnothing(reference_pressure) ? false : true - prescribed_velocity = isnothing(reference_velocity) ? false : true - prescribed_density = isnothing(reference_density) ? false : true + update_boundary_zone_indices!(system, initial_coordinates(system), boundary_zones, semi) - if boundary_model isa BoundaryModelTafuni - return (; prescribed_pressure=prescribed_pressure, - prescribed_density=prescribed_density, - prescribed_velocity=prescribed_velocity) - end + return system +end + +function create_cache_open_boundary(boundary_model, initial_condition, reference_values) + ELTYPE = eltype(initial_condition) + + # Separate `reference_values` into pressure, density and velocity reference values + pressure_reference_values = map(ref -> ref.reference_pressure, reference_values) + density_reference_values = map(ref -> ref.reference_density, reference_values) + velocity_reference_values = map(ref -> ref.reference_velocity, reference_values) - characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) - previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) + if boundary_model isa BoundaryModelCharacteristicsLastiwka + characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) + previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) - return (; characteristics=characteristics, - previous_characteristics=previous_characteristics, - prescribed_pressure=prescribed_pressure, - prescribed_density=prescribed_density, prescribed_velocity=prescribed_velocity) + return (; characteristics=characteristics, + previous_characteristics=previous_characteristics, + pressure_reference_values=pressure_reference_values, + density_reference_values=density_reference_values, + velocity_reference_values=velocity_reference_values) + else + return (; pressure_reference_values=pressure_reference_values, + density_reference_values=density_reference_values, + velocity_reference_values=velocity_reference_values) + end 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, boundary_type_name(system.boundary_zone)) print(io, ") with ", nparticles(system), " particles") end @@ -205,13 +154,9 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) summary_header(io, "OpenBoundarySPHSystem{$(ndims(system))}") summary_line(io, "#particles", nparticles(system)) summary_line(io, "#buffer_particles", system.buffer.buffer_size) + summary_line(io, "#boundary_zones", length(system.boundary_zones)) summary_line(io, "fluid system", type2string(system.fluid_system)) summary_line(io, "boundary model", type2string(system.boundary_model)) - 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)) - summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3)) summary_footer(io) end end @@ -220,13 +165,11 @@ end return ELTYPE end +@inline buffer(system::OpenBoundarySPHSystem) = system.buffer + # The `UpdateCallback` is required to update particle positions between time steps @inline requires_update_callback(system::OpenBoundarySPHSystem) = true -function corresponding_fluid_system(system::OpenBoundarySPHSystem, semi) - return system.fluid_system -end - function smoothing_length(system::OpenBoundarySPHSystem, particle) return system.smoothing_length end @@ -263,19 +206,17 @@ end # This function is called by the `UpdateCallback`, as the integrator array might be modified function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ode, semi, t) + (; boundary_model) = system + u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) - # Update density, pressure and velocity based on the characteristic variables. - # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 - @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, - system.boundary_model, - v, u, - v_ode, - u_ode, - semi, t) + # Update density, pressure and velocity based on the specific boundary model + @trixi_timeit timer() "update boundary quantities" begin + update_boundary_quantities!(system, boundary_model, v, u, v_ode, u_ode, semi, t) + end return system end @@ -283,8 +224,7 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) = system function check_domain!(system, v, u, v_ode, u_ode, semi) - (; boundary_zone, boundary_candidates, fluid_candidates) = system - fluid_system = corresponding_fluid_system(system, semi) + (; boundary_zones, boundary_candidates, fluid_candidates, fluid_system) = system u_fluid = wrap_u(u_ode, fluid_system, semi) v_fluid = wrap_v(v_ode, fluid_system, semi) @@ -296,6 +236,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle_coords = current_coords(u, system, particle) # Check if boundary particle is outside the boundary zone + boundary_zone = current_boundary_zone(system, particle) if !is_in_boundary_zone(boundary_zone, particle_coords) boundary_candidates[particle] = true end @@ -311,6 +252,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle = crossed_boundary_particles[i] particle_new = available_fluid_particles[i] + boundary_zone = current_boundary_zone(system, particle) convert_particle!(system, fluid_system, boundary_zone, particle, particle_new, v, u, v_fluid, u_fluid) end @@ -324,9 +266,11 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) @threaded semi for fluid_particle in each_moving_particle(fluid_system) fluid_coords = current_coords(u_fluid, fluid_system, fluid_particle) - # Check if fluid particle is in boundary zone - if is_in_boundary_zone(boundary_zone, fluid_coords) - fluid_candidates[fluid_particle] = true + # Check if fluid particle is in any boundary zone + for boundary_zone in boundary_zones + if is_in_boundary_zone(boundary_zone, fluid_coords) + fluid_candidates[fluid_particle] = true + end end end @@ -340,7 +284,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle = crossed_fluid_particles[i] particle_new = available_boundary_particles[i] - convert_particle!(fluid_system, system, boundary_zone, particle, particle_new, + convert_particle!(fluid_system, system, particle, particle_new, v, u, v_fluid, u_fluid) end @@ -350,39 +294,15 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) # Since particles have been transferred, the neighborhood searches must be updated update_nhs!(semi, u_ode) - return system -end - -# Outflow particle is outside the boundary zone -@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{OutFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) - deactivate_particle!(system, particle, u) - - return system -end - -# Inflow particle is outside the boundary zone -@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{InFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) - (; spanning_set) = boundary_zone - - # Activate a new particle in simulation domain - transfer_particle!(fluid_system, system, particle, particle_new, v_fluid, u_fluid, v, u) - - # Reset position of boundary particle - for dim in 1:ndims(system) - u[dim, particle] += spanning_set[1][dim] - end + update_boundary_zone_indices!(system, u, boundary_zones, semi) return system end # Buffer particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{BidirectionalFlow}, - particle, particle_new, v, u, v_fluid, u_fluid) + boundary_zone, particle, particle_new, + 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. @@ -406,8 +326,7 @@ end # Fluid particle is in boundary zone @inline function convert_particle!(fluid_system::FluidSystem, system, - boundary_zone, particle, particle_new, - v, u, v_fluid, u_fluid) + particle, particle_new, v, u, v_fluid, u_fluid) # Activate particle in boundary zone transfer_particle!(system, fluid_system, particle, particle_new, v, u, v_fluid, u_fluid) @@ -457,32 +376,6 @@ function write_u0!(u0, system::OpenBoundarySPHSystem) return u0 end -wrap_reference_function(::Nothing, ::Val) = nothing - -function wrap_reference_function(function_::Function, ::Val) - # Already a function - return function_ -end - -# Name the function so that the summary box does know which kind of function this is -function wrap_reference_function(constant_scalar_::Number, ::Val) - return constant_scalar(coords, t) = constant_scalar_ -end - -# For vectors and tuples -# Name the function so that the summary box does know which kind of function this is -function wrap_reference_function(constant_vector_, ::Val{NDIMS}) where {NDIMS} - return constant_vector(coords, t) = SVector{NDIMS}(constant_vector_) -end - -function reference_value(value::Function, quantity, position, t) - return value(position, t) -end - -# This method is used when extrapolating quantities from the domain -# instead of using the method of characteristics -reference_value(value::Nothing, quantity, position, t) = quantity - # To account for boundary effects in the viscosity term of the RHS, use the viscosity model # of the neighboring particle systems. @inline function viscosity_model(system::OpenBoundarySPHSystem, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 70dfd0cc72..a20ed70ac1 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -245,6 +245,8 @@ end return ndims(system) + 2 end +@inline buffer(system::EntropicallyDampedSPHSystem) = system.buffer + system_correction(system::EntropicallyDampedSPHSystem) = system.correction @inline function current_pressure(v, system::EntropicallyDampedSPHSystem, particle) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index dda2ead5f1..46a51cdb8c 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -238,6 +238,8 @@ end return ndims(system) + 1 end +@inline buffer(system::WeaklyCompressibleSPHSystem) = system.buffer + system_correction(system::WeaklyCompressibleSPHSystem) = system.correction @inline function current_velocity(v, system::WeaklyCompressibleSPHSystem) diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index bc93aa8c2f..c2b27ca75b 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -37,7 +37,6 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST, CM} <: SolidSystem{NDIM acceleration :: SVector{NDIMS, ELTYPE} source_terms :: ST contact_model :: CM - buffer :: Nothing function DEMSystem(initial_condition, contact_model; damping_coefficient=0.0001, acceleration=ntuple(_ -> 0.0, @@ -65,7 +64,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST, CM} <: SolidSystem{NDIM typeof(mass), typeof(source_terms), typeof(contact_model)}(initial_condition, mass, radius, damping_coefficient, acceleration_, source_terms, - contact_model, nothing) + contact_model) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index c0ef99b30a..60fc7a3990 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -78,7 +78,6 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, penalty_force :: PF viscosity :: V source_terms :: ST - buffer :: Nothing end function TotalLagrangianSPHSystem(initial_condition, @@ -123,7 +122,7 @@ function TotalLagrangianSPHSystem(initial_condition, n_moving_particles, young_modulus, poisson_ratio, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, - penalty_force, viscosity, source_terms, nothing) + penalty_force, viscosity, source_terms) end function Base.show(io::IO, system::TotalLagrangianSPHSystem) diff --git a/src/util.jl b/src/util.jl index 88be2096a2..48f5fe4be6 100644 --- a/src/util.jl +++ b/src/util.jl @@ -11,6 +11,17 @@ end @inline foreach_noalloc(func, collection::Tuple{}) = nothing +# Returns `functions[index](args...)`, but in a type-stable way for a heterogeneous tuple `functions` +@inline function apply_ith_function(functions, index, args...) + if index == 1 + # Found the function to apply, apply it and return + return first(functions)(args...) + end + + # Process remaining functions + apply_ith_function(Base.tail(functions), index - 1, args...) +end + # Print informative message at startup function print_startup_message() s = """ diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index a84a69317d..54085d2005 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -302,62 +302,52 @@ @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (WCSPH)" begin @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), joinpath(examples_dir(), "fluid", - "pipe_flow_2d.jl"), - wcsph=true) + "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (EDAC)" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (EDAC)" begin + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=false, + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (EDAC)" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (EDAC)" begin + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=false, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(), boundary_type_in=BidirectionalFlow(), - boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_velocity_out=nothing) + boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (WCSPH)" begin @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - wcsph=true, sound_speed=20.0, pressure=0.0, - open_boundary_model=BoundaryModelTafuni(), boundary_type_in=BidirectionalFlow(), - boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_pressure_out=nothing, - reference_velocity_out=nothing) + boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin - steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=10, + steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=5, reltol=1e-3) @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), extra_callback=steady_state_reached, tspan=(0.0, 1.5), viscosity_boundary=nothing) @@ -367,11 +357,12 @@ end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin - steady_state_reached = SteadyStateReachedCallback(; interval=1, interval_size=10, + steady_state_reached = SteadyStateReachedCallback(; interval=1, interval_size=5, reltol=1e-3) @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), extra_callback=steady_state_reached, dtmax=2e-3, tspan=(0.0, 1.5), viscosity_boundary=nothing) @@ -381,7 +372,7 @@ end @trixi_testset "fluid/pipe_flow_3d.jl" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_3d.jl")) @test sol.retcode == ReturnCode.Success diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 038e0a2b89..3c4db90621 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -305,7 +305,7 @@ end end # Test open boundaries and steady-state callback - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (WCSPH)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), @@ -318,7 +318,7 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (EDAC)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (EDAC)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), @@ -330,13 +330,13 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (EDAC)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (EDAC)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(), + open_boundary_model=BoundaryModelMirroringTafuni(), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, @@ -349,16 +349,15 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (WCSPH)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), wcsph=true, sound_speed=20.0f0, - pressure=0.0f0, - open_boundary_model=BoundaryModelTafuni(; - mirror_method=ZerothOrderMirroring()), + open_boundary_model=BoundaryModelMirroringTafuni(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, @@ -373,14 +372,14 @@ end end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin - steady_state_reached = SteadyStateReachedCallback(; dt=0.002f0, - interval_size=10, + steady_state_reached = SteadyStateReachedCallback(; dt=0.002f0, interval_size=5, reltol=1.0f-3) @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), extra_callback=steady_state_reached, tspan=(0.0f0, 1.5f0), parallelization_backend=Main.parallelization_backend, @@ -393,13 +392,14 @@ end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin steady_state_reached = SteadyStateReachedCallback(; interval=1, - interval_size=10, + interval_size=5, reltol=1.0f-3) @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), extra_callback=steady_state_reached, + open_boundary_model=BoundaryModelCharacteristicsLastiwka(), dtmax=2.0f-3, tspan=(0.0f0, 1.5f0), parallelization_backend=Main.parallelization_backend, diff --git a/test/general/buffer.jl b/test/general/buffer.jl index 18758590ea..b953d7cfa8 100644 --- a/test/general/buffer.jl +++ b/test/general/buffer.jl @@ -6,15 +6,13 @@ 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()) + reference_density=1.0, reference_pressure=0.0, + reference_velocity=[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) + boundary_model=BoundaryModelCharacteristicsLastiwka(), + buffer_size=0) system_buffer = OpenBoundarySPHSystem(zone; buffer_size=5, - reference_density=0.0, reference_pressure=0.0, - reference_velocity=[0, 0], - boundary_model=BoundaryModelLastiwka(), + boundary_model=BoundaryModelCharacteristicsLastiwka(), fluid_system=FluidSystemMock3()) n_particles = nparticles(system) diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index dea8c356ab..862f2ba37a 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -140,7 +140,7 @@ v_ode = vcat(vec(v1), v2) # Avoid `SystemBuffer` barrier - TrixiParticles.each_moving_particle(system::Union{System1, System2}) = TrixiParticles.eachparticle(system) + TrixiParticles.each_moving_particle(system::Union{System1, System2}) = Base.OneTo(nparticles(system)) TrixiParticles.add_source_terms!(dv_ode, v_ode, u_ode, semi, 0.0) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index b14906b630..ce80fced84 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -1,4 +1,94 @@ @testset verbose=true "Boundary Zone" begin + @testset "`show`" begin + inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + plane_normal=(1.0, 0.0), density=1.0, + reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + open_boundary_layers=4, boundary_type=InFlow()) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(inflow) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… inflow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", inflow) == show_box + + outflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + plane_normal=(1.0, 0.0), density=1.0, open_boundary_layers=4, + boundary_type=OutFlow()) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(outflow) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… outflow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", outflow) == show_box + end + + @testset verbose=true "Illegal Inputs" begin + plane = ([0.0, 0.0], [0.0, 1.0]) + flow_direction = (1.0, 0.0) + + error_str = "`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 2 for a 2D problem holding this velocity" + + reference_velocity = 1.0 + @test_throws ArgumentError(error_str) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density=0, + reference_pressure=0, + reference_velocity, + open_boundary_layers=2, + boundary_type=InFlow()) + + 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) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density=0, + reference_velocity=[1.0, + 1.0], reference_pressure, + open_boundary_layers=2, + boundary_type=InFlow()) + + 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) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density, + reference_velocity=[1.0, + 1.0], + reference_pressure=0, + open_boundary_layers=2, + boundary_type=InFlow()) + end @testset verbose=true "Boundary Zone 2D" begin particle_spacing = 0.2 open_boundary_layers = 4 @@ -35,8 +125,8 @@ zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (first(typeof(boundary_zone).parameters) === - TrixiParticles.InFlow) ? -1 : 1 + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? + -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -99,8 +189,8 @@ zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (first(typeof(boundary_zone).parameters) === - TrixiParticles.InFlow) ? -1 : 1 + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? + -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -137,7 +227,7 @@ @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones - perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ? sqrt(eps()) : -sqrt(eps()) @@ -183,7 +273,7 @@ @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones - perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ? eps() : -eps() point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 570818f8f9..1920521caa 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -15,7 +15,8 @@ # Prescribed quantities reference_velocity = (pos, t) -> SVector(t, 0.0) reference_pressure = (pos, t) -> 50_000.0 * t - reference_density = (pos, t) -> 1000.0 * t + # Add small offset to avoid "ArgumentError: density must be positive and larger than `eps()`" + reference_density = (pos, t) -> 1000.0 * (t + sqrt(eps())) # Plane points of open boundary plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] @@ -36,10 +37,12 @@ flow_direction = flow_directions[j] inflow = BoundaryZone(; plane=plane_points, particle_spacing, density, plane_normal=flow_direction, open_boundary_layers, - boundary_type=InFlow()) + boundary_type=InFlow(), reference_velocity, + reference_pressure, reference_density) outflow = BoundaryZone(; plane=plane_points, particle_spacing, density, plane_normal=(-flow_direction), open_boundary_layers, - boundary_type=OutFlow()) + boundary_type=OutFlow(), reference_velocity, + reference_pressure, reference_density) boundary_zones = [ inflow, @@ -49,7 +52,7 @@ @testset "`$(TrixiParticles.boundary_type_name(boundary_zone))`" for boundary_zone in boundary_zones - sign_ = (first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow) ? + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? 1 : -1 fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, density, pressure, @@ -62,10 +65,7 @@ boundary_system = OpenBoundarySPHSystem(boundary_zone; fluid_system, buffer_size=0, - boundary_model=BoundaryModelLastiwka(), - reference_velocity, - reference_pressure, - reference_density) + boundary_model=BoundaryModelCharacteristicsLastiwka()) semi = Semidiscretization(fluid_system, boundary_system) @@ -101,13 +101,13 @@ v, u, v0_ode, u0_ode, semi, t1) evaluated_vars1 = boundary_system.cache.characteristics - if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow + if TrixiParticles.boundary_type_name(boundary_zone) == "inflow" @test all(isapprox.(evaluated_vars1[1, :], 0.0)) @test all(isapprox.(evaluated_vars1[2, :], 0.0)) @test all(isapprox.(evaluated_vars1[3, 1:n_influenced], J3(t1))) @test all(isapprox.(evaluated_vars1[3, (n_influenced + 1):end], 0.0)) - elseif first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow + elseif TrixiParticles.boundary_type_name(boundary_zone) == "outflow" @test all(isapprox.(evaluated_vars1[1, 1:n_influenced], J1(t1))) @test all(isapprox.(evaluated_vars1[2, 1:n_influenced], J2(t1))) @test all(isapprox.(evaluated_vars1[1:2, (n_influenced + 1):end], 0.0)) @@ -121,13 +121,13 @@ v, u, v0_ode, u0_ode, semi, t2) evaluated_vars2 = boundary_system.cache.characteristics - if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow + if TrixiParticles.boundary_type_name(boundary_zone) == "inflow" @test all(isapprox.(evaluated_vars2[1, :], 0.0)) @test all(isapprox.(evaluated_vars2[2, :], 0.0)) @test all(isapprox.(evaluated_vars2[3, 1:n_influenced], J3(t2))) @test all(isapprox.(evaluated_vars2[3, (n_influenced + 1):end], J3(t1))) - elseif first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow + elseif TrixiParticles.boundary_type_name(boundary_zone) == "outflow" @test all(isapprox.(evaluated_vars2[1, 1:n_influenced], J1(t2))) @test all(isapprox.(evaluated_vars2[2, 1:n_influenced], J2(t2))) @test all(isapprox.(evaluated_vars2[1, (n_influenced + 1):end], J1(t1))) diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 6e3fb9cd6f..ae7afc26bf 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -61,11 +61,12 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelMirroringTafuni(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary, semi) v_open_boundary = zero(inflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -75,9 +76,7 @@ TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false, - prescribed_velocity=false) + domain_fluid.coordinates, semi) # Checked visually in ParaView: # trixi2vtk(fluid_system.initial_condition, filename="fluid", # v=domain_fluid.velocity, p=domain_fluid.pressure) @@ -158,11 +157,12 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelMirroringTafuni(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary, semi) v_open_boundary = zero(inflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -172,9 +172,7 @@ TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false, - prescribed_velocity=false) + domain_fluid.coordinates, semi) # Checked visually in ParaView: # trixi2vtk(fluid_system.initial_condition, filename="fluid", # v=domain_fluid.velocity, p=domain_fluid.pressure) @@ -231,11 +229,12 @@ open_boundary_layers=open_boundary_layers, density=1000.0, particle_spacing, average_inflow_velocity=true) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelMirroringTafuni(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary_in) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_in, semi) v_open_boundary = zero(inflow.initial_condition.velocity) u_open_boundary = inflow.initial_condition.coordinates @@ -246,10 +245,10 @@ TrixiParticles.extrapolate_values!(open_boundary_in, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0) + domain_fluid.coordinates, semi) TrixiParticles.average_velocity!(v_open_boundary, u_open_boundary, open_boundary_in, - inflow, semi) + first(open_boundary_in.boundary_zones), semi) # Since the velocity profile increases linearly in positive x-direction, # we can use the first velocity entry as a representative value. @@ -282,12 +281,13 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelMirroringTafuni(), buffer_size=0) # Temporary semidiscretization just to extrapolate the pressure into the outflow system semi = Semidiscretization(fluid_system, open_boundary_out) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_out, semi) v_open_boundary = zero(outflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -297,8 +297,7 @@ TrixiParticles.extrapolate_values!(open_boundary_out, mirror_method, v_open_boundary, v_fluid, outflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false) + domain_fluid.coordinates, semi) plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) @@ -306,12 +305,13 @@ plane_normal=[1.0, 0.0], open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelMirroringTafuni(), buffer_size=0) # Temporary semidiscretization just to extrapolate the pressure into the outflow system semi = Semidiscretization(fluid_system, open_boundary_in) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_in, semi) v_open_boundary = zero(inflow.initial_condition.velocity) @@ -320,8 +320,7 @@ TrixiParticles.extrapolate_values!(open_boundary_in, mirror_method, v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false) + domain_fluid.coordinates, semi) return fluid_system, open_boundary_in, open_boundary_out, v_fluid end @@ -335,7 +334,7 @@ v_fluid = mirror(pressure_func, mirror_method) p_fluid = [TrixiParticles.current_pressure(v_fluid, fluid_system, particle) - for particle in TrixiParticles.active_particles(fluid_system)] + for particle in TrixiParticles.eachparticle(fluid_system)] fluid_system.initial_condition.pressure .= p_fluid open_boundary_in.initial_condition.pressure .= open_boundary_in.pressure diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index 0fdf98f907..3ce14e7634 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -1,71 +1,19 @@ @testset verbose=true "`OpenBoundarySPHSystem`" begin - @testset verbose=true "Illegal Inputs" begin - plane = ([0.0, 0.0], [0.0, 1.0]) - flow_direction = (1.0, 0.0) + @testset "`show`" begin # Mock fluid system struct FluidSystemMock2 <: TrixiParticles.FluidSystem{2} end TrixiParticles.initial_smoothing_length(system::FluidSystemMock2) = 1.0 TrixiParticles.nparticles(system::FluidSystemMock2) = 1 - inflow = BoundaryZone(; plane, particle_spacing=0.1, - plane_normal=flow_direction, density=1.0, - open_boundary_layers=2, boundary_type=InFlow()) - - error_str = "`reference_velocity` must be either a function mapping " * - "each particle's coordinates and time to its velocity, " * - "an array where the ``i``-th column holds the velocity of particle ``i`` " * - "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" - - reference_velocity = 1.0 - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density=0, - reference_pressure=0, - reference_velocity) - - error_str = "`reference_pressure` must be either a function mapping " * - "each particle's coordinates and time to its pressure, " * - "a vector holding the pressure of each particle, or a scalar" - - reference_pressure = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density=0, - reference_velocity=[1.0, - 1.0], - reference_pressure) - - error_str = "`reference_density` must be either a function mapping " * - "each particle's coordinates and time to its density, " * - "a vector holding the density of each particle, or a scalar" - - reference_density = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density, - reference_velocity=[1.0, - 1.0], - reference_pressure=0) - end - @testset "`show`" begin inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, plane_normal=(1.0, 0.0), density=1.0, open_boundary_layers=4, boundary_type=InFlow()) system = OpenBoundarySPHSystem(inflow; buffer_size=0, - boundary_model=BoundaryModelLastiwka(), - reference_density=0.0, - reference_pressure=0.0, - reference_velocity=[0.0, 0.0], + boundary_model=BoundaryModelCharacteristicsLastiwka(), fluid_system=FluidSystemMock2()) - show_compact = "OpenBoundarySPHSystem{2}(InFlow) with 80 particles" + show_compact = "OpenBoundarySPHSystem{2}() with 80 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -73,28 +21,21 @@ │ ════════════════════════ │ │ #particles: ………………………………………………… 80 │ │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 1 │ │ fluid system: …………………………………………… FluidSystemMock2 │ - │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary type: ………………………………………… InFlow │ - │ prescribed velocity: ………………………… constant_vector │ - │ prescribed pressure: ………………………… constant_scalar │ - │ prescribed density: …………………………… constant_scalar │ - │ width: ……………………………………………………………… 0.2 │ + │ boundary model: ……………………………………… BoundaryModelCharacteristicsLastiwka │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box - outflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + outflow = BoundaryZone(; plane=([5.0, 0.0], [5.0, 1.0]), particle_spacing=0.05, plane_normal=(1.0, 0.0), density=1.0, open_boundary_layers=4, boundary_type=OutFlow()) system = OpenBoundarySPHSystem(outflow; buffer_size=0, - boundary_model=BoundaryModelLastiwka(), - reference_density=0.0, - reference_pressure=0.0, - reference_velocity=[0.0, 0.0], + boundary_model=BoundaryModelMirroringTafuni(), fluid_system=FluidSystemMock2()) - show_compact = "OpenBoundarySPHSystem{2}(OutFlow) with 80 particles" + show_compact = "OpenBoundarySPHSystem{2}() with 80 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -102,13 +43,28 @@ │ ════════════════════════ │ │ #particles: ………………………………………………… 80 │ │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 1 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary model: ……………………………………… BoundaryModelMirroringTafuni │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", system) == show_box + + system = OpenBoundarySPHSystem(outflow, inflow; buffer_size=0, + boundary_model=BoundaryModelMirroringTafuni(), + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}() with 160 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 160 │ + │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 2 │ │ fluid system: …………………………………………… FluidSystemMock2 │ - │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary type: ………………………………………… OutFlow │ - │ prescribed velocity: ………………………… constant_vector │ - │ prescribed pressure: ………………………… constant_scalar │ - │ prescribed density: …………………………… constant_scalar │ - │ width: ……………………………………………………………… 0.2 │ + │ boundary model: ……………………………………… BoundaryModelMirroringTafuni │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box