diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl index dd324a5bca..5e57c7c0bd 100644 --- a/examples/fluid/poiseuille_flow_2d.jl +++ b/examples/fluid/poiseuille_flow_2d.jl @@ -136,6 +136,7 @@ outflow = BoundaryZone(; boundary_face=face_out, face_normal=(.-(flow_direction) open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, boundary_model=open_boundary_model, + calculate_flow_rate=true, buffer_size=n_buffer_particles) # ========================================================================================== diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 01397a7732..6a58f7c3f9 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -9,3 +9,29 @@ function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system, foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search; points, parallelization_backend) end + +deactivate_out_of_bounds_particles!(system, ::Nothing, nhs, u, semi) = system + +function deactivate_out_of_bounds_particles!(system, ::SystemBuffer, + cell_list::FullGridCellList, u, semi) + (; min_corner, max_corner) = cell_list + (; cell_size) = get_neighborhood_search(system, semi) + + # Remove the padding layer (see comment in PointNeighbors: full_grid.jl line 60) + min_corner_ = min_corner .+ 1001 // 1000 .* cell_size + max_corner_ = max_corner .- 1001 // 1000 .* cell_size + + @threaded semi for particle in each_integrated_particle(system) + particle_position = current_coords(u, system, particle) + + if !all(min_corner_ .<= particle_position .<= max_corner_) + deactivate_particle!(system, particle, u) + end + end + + if count(system.buffer.active_particle) != system.buffer.active_particle_count[] + update_system_buffer!(system.buffer, semi) + end + + return system +end diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl index 6206bd25f3..884588033b 100644 --- a/src/io/read_vtk.jl +++ b/src/io/read_vtk.jl @@ -31,9 +31,11 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu")) │ eltype: …………………………………………………………… Float64 │ │ coordinate eltype: ……………………………… Float64 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘ -```` +``` """ -function vtk2trixi(file) +function vtk2trixi(file; element_type=Float64, coordinates_eltype=Float64) + ELTYPE = element_type + cELTYPE = coordinates_eltype vtk_file = ReadVTK.VTKFile(file) # Retrieve data fields (e.g., pressure, velocity, ...) @@ -42,7 +44,7 @@ function vtk2trixi(file) # Retrieve fields ndims = first(ReadVTK.get_data(field_data["ndims"])) - coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :] + coordinates = convert.(cELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :]) fields = ["velocity", "density", "pressure", "mass", "particle_spacing"] results = Dict{String, Array{Float64}}() @@ -52,11 +54,11 @@ function vtk2trixi(file) all_keys = keys(point_data) idx = findfirst(k -> occursin(field, k), all_keys) if idx !== nothing - results[field] = ReadVTK.get_data(point_data[all_keys[idx]]) + results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]])) else # Use zeros as default values when a field is missing results[field] = field in ["mass"] ? - zeros(size(coordinates, 2)) : zero(coordinates) + zeros(ELTYPE, size(coordinates, 2)) : zero(coordinates) @info "No '$field' field found in VTK file. Will be set to zero." end end @@ -65,7 +67,8 @@ function vtk2trixi(file) first(results["particle_spacing"]) : results["particle_spacing"] - return InitialCondition(; coordinates, particle_spacing=particle_spacing, + return InitialCondition(; coordinates, + particle_spacing=convert(ELTYPE, particle_spacing), velocity=results["velocity"], mass=results["mass"], density=results["density"], diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 07fa559d2e..ff2f2f85eb 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -409,6 +409,16 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem) vtk["pressure"] = [current_pressure(v, system, particle) for particle in eachparticle(system)] + if system.cache.calculate_flow_rate + Q_total = zero(eltype(system)) + for i in eachindex(system.cache.boundary_zones_flow_rate) + vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i][] + Q_total += system.cache.boundary_zones_flow_rate[i][] + end + + vtk["Q_total"] = Q_total + end + return vtk end diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 7088fea58e..55b86be737 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -81,6 +81,10 @@ There are three ways to specify the actual shape of the boundary zone: Per default it is set to zero (assuming a gauge pressure system). - For `EntropicallyDampedSPHSystem`: Use the initial pressure from the `InitialCondition` - For `WeaklyCompressibleSPHSystem`: Use the background pressure from the equation of state +- `sample_points=:default`: Either `:default` to automatically generate sample points on the boundary face (default), + or a matrix of dimensions `(ndims, n_points)` containing sample points + on the boundary face used to compute the volumetric flow rate. + Set to `nothing` to skip sampling. !!! note "Note" The reference values (`reference_velocity`, `reference_pressure`, `reference_density`) @@ -147,7 +151,7 @@ bidirectional_flow = BoundaryZone(; boundary_face=face_vertices, face_normal, !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R} +struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R, C} initial_condition :: IC spanning_set :: S zone_origin :: ZO @@ -156,6 +160,7 @@ struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R} face_normal :: FN rest_pressure :: ELTYPE # Only required for `BoundaryModelDynamicalPressureZhang` reference_values :: R + cache :: C # Note that the following can't be static type parameters, as all boundary zones in a system # must have the same type, so that we can loop over them in a type-stable way. average_inflow_velocity :: Bool @@ -167,7 +172,7 @@ end function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, open_boundary_layers::Integer, average_inflow_velocity=true, - boundary_type=BidirectionalFlow(), + boundary_type=BidirectionalFlow(), sample_points=:default, rest_pressure=zero(eltype(density)), reference_density=nothing, reference_pressure=nothing, reference_velocity=nothing) @@ -262,10 +267,13 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, ic.velocity .= stack(velocity_ref.(coordinates_svector, 0)) end + cache = (; + create_cache_boundary_zone(ic, boundary_face, face_normal_, sample_points)...) + return BoundaryZone(ic, spanning_set_, zone_origin, zone_width, flow_direction, face_normal_, rest_pressure, reference_values, - average_inflow_velocity, prescribed_density, prescribed_pressure, - prescribed_velocity) + cache, average_inflow_velocity, prescribed_density, + prescribed_pressure, prescribed_velocity) end function boundary_type_name(boundary_zone::BoundaryZone) @@ -297,10 +305,57 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::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=6)) + if hasproperty(boundary_zone.cache, :cross_sectional_area) + summary_line(io, "cross sectional area", + boundary_zone.cache.cross_sectional_area) + end summary_footer(io) end end +function create_cache_boundary_zone(initial_condition, boundary_face, face_normal, + sample_points::Nothing) + return (; sample_points) +end + +function create_cache_boundary_zone(initial_condition, boundary_face, face_normal, + sample_points) + (; particle_spacing) = initial_condition + area_increment = particle_spacing^(ndims(initial_condition) - 1) + if sample_points === :default + sample_points_ = extrude_geometry(boundary_face; particle_spacing, density=Inf, + direction=(-face_normal), n_extrude=1).coordinates + else + if !(sample_points isa Matrix && size(sample_points, 1) == ndims(initial_condition)) + throw(ArgumentError("`sample_points` must be a matrix with " * + "`ndims(initial_condition)` rows")) + end + + # TODO: Check if regular grid? + + sample_points_ = sample_points + end + + discrete_face_area = area_increment * size(sample_points_, 2) + + if ndims(initial_condition) == 3 + v1, v2, v3 = boundary_face + face_area = norm(cross(v2 - v1, v3 - v1)) + elseif ndims(initial_condition) == 2 + v1, v2 = boundary_face + face_area = norm(v2 - v1) + end + + if discrete_face_area > (face_area + eps(face_area)) + @warn "The sampled area of the boundary face " * + "($(discrete_face_area)) is larger than the actual face area ($(face_area)). " + end + + return (; sample_points=sample_points_, sample_velocity=copy(sample_points_), + shepard_coefficient=zeros(eltype(initial_condition), size(sample_points_, 2)), + area_increment, cross_sectional_area=discrete_face_area) +end + function set_up_boundary_zone(boundary_face, face_normal, density, particle_spacing, initial_condition, extrude_geometry, open_boundary_layers, boundary_type) diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl index 81ce8f9ff4..7e180396a1 100644 --- a/src/schemes/boundary/open_boundary/pressure_model.jl +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -67,51 +67,32 @@ function update_pressure_model!(system, v, u, semi, dt) if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values) @trixi_timeit timer() "update pressure model" begin - calculate_flow_rate_and_pressure!(system, v, u, dt) + calculate_pressure!(system, dt) end end return system end -function calculate_flow_rate_and_pressure!(system, v, u, dt) - (; pressure_reference_values) = system.cache - foreach_enumerate(pressure_reference_values) do (zone_id, pressure_model) - boundary_zone = system.boundary_zones[zone_id] - calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, v, u, dt) +function calculate_pressure!(system, dt) + (; pressure_reference_values, boundary_zones_flow_rate) = system.cache + + foreach_noalloc(pressure_reference_values, + boundary_zones_flow_rate) do (pressure_model, flow_rate) + calculate_pressure!(pressure_model, system, flow_rate[], dt) end return system end -function calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, v, u, dt) +function calculate_pressure!(pressure_model, system, current_flow_rate, dt) return pressure_model end -function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system, - boundary_zone, v, u, dt) - (; particle_spacing) = system.initial_condition +function calculate_pressure!(pressure_model::RCRWindkesselModel, system, + current_flow_rate, dt) (; characteristic_resistance, peripheral_resistance, compliance, flow_rate, pressure) = pressure_model - (; face_normal) = boundary_zone - - # Find particles within the current boundary zone - candidates = findall(particle -> boundary_zone == - current_boundary_zone(system, particle), - each_integrated_particle(system)) - - # Assuming negligible transverse velocity gradients within the boundary zone, - # the full area of the zone is taken as the representative cross-sectional - # area for volumetric flow-rate estimation. - cross_sectional_area = length(candidates) * particle_spacing^(ndims(system) - 1) - - # Division inside the `sum` closure to maintain GPU compatibility - velocity_avg = sum(candidates) do particle - return dot(current_velocity(v, system, particle), -face_normal) / length(candidates) - end - - # Compute volumetric flow rate: Q = v * A - current_flow_rate = velocity_avg * cross_sectional_area previous_pressure = pressure[] previous_flow_rate = flow_rate[] diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9396c3f1fb..bd05cf16f7 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,7 +1,7 @@ @doc raw""" OpenBoundarySystem(boundary_zone::BoundaryZone; fluid_system::AbstractFluidSystem, buffer_size::Integer, - boundary_model) + boundary_model, calculate_flow_rate=false) Open boundary system for in- and outflow particles. @@ -11,6 +11,8 @@ Open boundary system for in- and outflow particles. # Keywords - `fluid_system`: The corresponding fluid system - `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models)) +- `calculate_flow_rate=false`: Set to `true` to calculate the volumetric flow rate through each boundary zone. + Default is `false`. - `buffer_size`: Number of buffer particles. - `pressure_acceleration`: Pressure acceleration formulation for the system. Required only when using [`BoundaryModelDynamicalPressureZhang`](@ref). @@ -63,7 +65,8 @@ end function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; fluid_system::AbstractFluidSystem, buffer_size::Integer, - boundary_model, + boundary_model, calculate_flow_rate=false, + deactivate_isolated_particles=false, pressure_acceleration=boundary_model isa BoundaryModelDynamicalPressureZhang ? fluid_system.pressure_acceleration_formulation : @@ -82,10 +85,12 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; mass = copy(initial_conditions.mass) volume = similar(initial_conditions.density) - cache = (; + neighbor_counter = deactivate_isolated_particles ? + zeros(Int, nparticles(initial_conditions)) : nothing + cache = (; neighbor_counter=neighbor_counter, create_cache_shifting(initial_conditions, shifting_technique)..., create_cache_open_boundary(boundary_model, fluid_system, initial_conditions, - boundary_zones_)...) + calculate_flow_rate, boundary_zones_)...) fluid_system_index = Ref(0) @@ -110,6 +115,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; zone.face_normal, zone.rest_pressure, nothing, + zone.cache, zone.average_inflow_velocity, zone.prescribed_density, zone.prescribed_pressure, @@ -132,7 +138,7 @@ function initialize!(system::OpenBoundarySystem, semi) end function create_cache_open_boundary(boundary_model, fluid_system, initial_condition, - boundary_zones) + calculate_flow_rate, boundary_zones) reference_values = map(bz -> bz.reference_values, boundary_zones) ELTYPE = eltype(initial_condition) @@ -141,6 +147,26 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit density_reference_values = map(ref -> ref.reference_density, reference_values) velocity_reference_values = map(ref -> ref.reference_velocity, reference_values) + if any(pr -> isa(pr, AbstractPressureModel), pressure_reference_values) + calculate_flow_rate = true + end + + cache = (; pressure_reference_values=pressure_reference_values, + density_reference_values=density_reference_values, + velocity_reference_values=velocity_reference_values, calculate_flow_rate) + + if calculate_flow_rate + if any(zone -> isnothing(zone.cache.sample_points), boundary_zones) + throw(ArgumentError("`sample_points` must be specified for all boundary zones when " * + "`calculate_flow_rate` is true.\n" * + "Use `sample_points=:default` to automatically generate sample points.")) + end + + boundary_zones_flow_rate = ntuple(i -> Ref(zero(ELTYPE)), + Val(length(boundary_zones))) + cache = (; boundary_zones_flow_rate, cache...) + end + if boundary_model isa BoundaryModelCharacteristicsLastiwka characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) @@ -148,10 +174,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit return (; characteristics=characteristics, previous_characteristics=previous_characteristics, pressure=copy(initial_condition.pressure), - density=copy(initial_condition.density), - pressure_reference_values=pressure_reference_values, - density_reference_values=density_reference_values, - velocity_reference_values=velocity_reference_values) + density=copy(initial_condition.density), cache...) + elseif boundary_model isa BoundaryModelDynamicalPressureZhang # A separate array for the boundary pressure is required, # since it is specified independently from the computed pressure for the momentum equation. @@ -171,10 +195,7 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit cache = (; density_calculator=ContinuityDensity(), density_diffusion=density_diffusion_, pressure_boundary=pressure_boundary, - density_rest=density_rest, - pressure_reference_values=pressure_reference_values, - density_reference_values=density_reference_values, - velocity_reference_values=velocity_reference_values) + density_rest=density_rest, cache...) if fluid_system isa EntropicallyDampedSPHSystem # Density and pressure is stored in `v` @@ -186,10 +207,7 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit else return (; pressure=copy(initial_condition.pressure), - density=copy(initial_condition.density), - pressure_reference_values=pressure_reference_values, - density_reference_values=density_reference_values, - velocity_reference_values=velocity_reference_values) + density=copy(initial_condition.density), cache...) end end @@ -290,6 +308,12 @@ end return du end +function update_positions!(system::OpenBoundarySystem, v, u, v_ode, u_ode, semi, t) + cell_list = get_neighborhood_search(system, semi).cell_list + deactivate_out_of_bounds_particles!(system, buffer(system), cell_list, u, semi) +end + + function update_boundary_interpolation!(system::OpenBoundarySystem, v, u, v_ode, u_ode, semi, t) update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) @@ -305,6 +329,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) + calculate_neighbor_count!(system, shifting_technique(system), u, semi) + + # This must be called before `update_pressure_model!` and `check_domain!` + # to ensure quantities remain consistent with the current simulation state. + calculate_flow_rate!(system, v_ode, u_ode, semi) + @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) update_pressure_model!(system, v, u, semi, integrator.dt) @@ -323,6 +353,35 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system +function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, + v_ode, u_ode, semi) where {ELTYPE, NDIMS} + system.cache.calculate_flow_rate || return system + + @trixi_timeit timer() "flow rate calculation" begin + (; boundary_zones) = system + (; boundary_zones_flow_rate) = system.cache + + foreach_noalloc(boundary_zones, + boundary_zones_flow_rate) do (boundary_zone, flow_rate) + (; face_normal) = boundary_zone + (; sample_velocity, area_increment) = boundary_zone.cache + + interpolate_velocity!(system, boundary_zone, v_ode, u_ode, semi) + + # Compute volumetric flow rate: Q = ∫ v ⋅ n dA + velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity) + current_flow_rate = sum(velocities) do velocity + vn = dot(velocity, -face_normal) + return vn * area_increment + end + + flow_rate[] = current_flow_rate + end + end + + return system +end + function check_domain!(system, v, u, v_ode, u_ode, semi) (; boundary_zones, boundary_candidates, fluid_candidates, fluid_system) = system @@ -418,6 +477,12 @@ end return system end + if is_isolated(system, shifting_technique(system), particle) + deactivate_particle!(system, particle, u) + + return system + end + # Activate a new particle in simulation domain transfer_particle!(fluid_system, system, particle, particle_new, v_fluid, u_fluid, v, u) @@ -570,3 +635,87 @@ end return system end + +function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, + v_ode, u_ode, semi) + (; sample_points, sample_velocity, shepard_coefficient) = boundary_zone.cache + smoothing_length = initial_smoothing_length(system) + smoothing_kernel = system_smoothing_kernel(system) + + set_zero!(shepard_coefficient) + set_zero!(sample_velocity) + + foreach_system(semi) do neighbor_system + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, sample_points, neighbor_coords, + semi, + points=axes(sample_points, 2)) do point, neighbor, + pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + volume_b = m_b / current_density(v_neighbor, neighbor_system, neighbor) + W_ab = kernel(smoothing_kernel, distance, smoothing_length) + shepard_coefficient[point] += volume_b * W_ab + + velocity_neighbor = viscous_velocity(v_neighbor, neighbor_system, neighbor) + for i in axes(velocity_neighbor, 1) + sample_velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab + end + end + end + + @threaded semi for point in axes(sample_points, 2) + if shepard_coefficient[point] > eps() + for i in axes(sample_velocity, 1) + sample_velocity[i, point] /= shepard_coefficient[point] + end + end + end + + return system +end + +@inline calculate_neighbor_count!(system, ::Nothing, u, semi) = system + +@inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi) + (; neighbor_counter) = system.cache + + # Skip when `deactivate_isolated_particles` is not activated + isnothing(neighbor_counter) && return system + + set_zero!(neighbor_counter) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + system_coords = current_coordinates(u, system) + foreach_point_neighbor(system, system, system_coords, system_coords, semi; + points=each_integrated_particle(system)) do particle, + neighbor, + pos_diff, + distance + (particle == neighbor) && return + neighbor_counter[particle] += 1 + end + + return system +end + +@inline is_isolated(system, ::Nothing, particle) = false + +@inline function is_isolated(system, ::AbstractShiftingTechnique, particle) + (; particle_spacing) = system.initial_condition + (; neighbor_counter) = system.cache + + # Skip when `deactivate_isolated_particles` is not activated + isnothing(neighbor_counter) && return false + + min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, + compact_support(system, system)) / 10 + + if neighbor_counter[particle] < min_neighbors + return true + end + + return false +end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 35d14b5e04..67c40c8d3a 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -115,6 +115,11 @@ end @inline acceleration_source(system::AbstractFluidSystem) = system.acceleration +function update_positions!(system::AbstractFluidSystem, v, u, v_ode, u_ode, semi, t) + cell_list = get_neighborhood_search(system, semi).cell_list + deactivate_out_of_bounds_particles!(system, buffer(system), cell_list, u, semi) +end + function compute_density!(system, u, u_ode, semi, ::ContinuityDensity) # No density update with `ContinuityDensity` return system diff --git a/src/util.jl b/src/util.jl index 18389cd44b..305ebc953a 100644 --- a/src/util.jl +++ b/src/util.jl @@ -11,21 +11,20 @@ end @inline foreach_noalloc(func, collection::Tuple{}) = nothing -# Same as `foreach(enumerate(something))`, but without allocations. -# Note that compile times may increase if this is used with big tuples. -@inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1) -@inline foreach_enumerate(func, collection::Tuple{}, index) = nothing +@inline function foreach_noalloc(func, collection1, collection2) + element1 = first(collection1) + remaining_collection1 = Base.tail(collection1) + element2 = first(collection2) + remaining_collection2 = Base.tail(collection2) -@inline function foreach_enumerate(func, collection, index) - element = first(collection) - remaining_collection = Base.tail(collection) - - @inline func((index, element)) + func((element1, element2)) # Process remaining collection - foreach_enumerate(func, remaining_collection, index + 1) + foreach_noalloc(func, remaining_collection1, remaining_collection2) end +@inline foreach_noalloc(func, collection1::Tuple{}, collection2::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 diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index 619c048bbe..3e1fa13118 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -1,7 +1,7 @@ @testset verbose=true "Boundary Zone" begin @testset "`show`" begin inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), - particle_spacing=0.05, + particle_spacing=0.05, sample_points=nothing, face_normal=(1.0, 0.0), density=1.0, reference_density=0.0, reference_pressure=0.0, @@ -22,7 +22,7 @@ @test repr("text/plain", inflow) == show_box outflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), - particle_spacing=0.05, + particle_spacing=0.05, sample_points=nothing, reference_density=0.0, reference_pressure=0.0, reference_velocity=[0.0, 0.0], @@ -41,6 +41,48 @@ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", outflow) == show_box + + bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), + particle_spacing=0.05, sample_points=nothing, + reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + face_normal=(1.0, 0.0), density=1.0, + open_boundary_layers=4) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(bidirectional) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… bidirectional_flow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", bidirectional) == show_box + + zone = BoundaryZone(; boundary_face=([0.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)]), + particle_spacing=0.05, reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + face_normal=normalize([-1.0, 1.0]), density=1.0, + open_boundary_layers=4) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(zone) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… bidirectional_flow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + │ cross sectional area: ……………………… 1.0 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", zone) == show_box end @testset verbose=true "Illegal Inputs" begin @@ -127,9 +169,7 @@ outflow ] - @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in - boundary_zones - + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? @@ -192,9 +232,7 @@ outflow ] - @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in - boundary_zones - + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? @@ -232,9 +270,7 @@ outflow ] - @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in - boundary_zones - + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ? sqrt(eps()) : -sqrt(eps()) @@ -280,9 +316,7 @@ outflow ] - @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in - boundary_zones - + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones 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/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl index f35c2e867f..0fe47f0a59 100644 --- a/test/schemes/boundary/open_boundary/pressure_model.jl +++ b/test/schemes/boundary/open_boundary/pressure_model.jl @@ -1,8 +1,7 @@ @testset verbose=true "`RCRWindkesselModel`" begin @testset verbose=true "Show" begin pressure_model = RCRWindkesselModel(; peripheral_resistance=1.2, - compliance=0.5, - characteristic_resistance=2.3) + compliance=0.5, characteristic_resistance=2.3) show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -75,21 +74,19 @@ particle_spacing=0.1, face_normal=(-1.0, 0.0), density=1000.0, reference_velocity, reference_pressure=pressure_model, - open_boundary_layers=1, rest_pressure=p_0) + open_boundary_layers=17, rest_pressure=p_0) system = OpenBoundarySystem(boundary_zone; buffer_size=0, boundary_model=nothing, fluid_system=FluidSystemMockRCR(nothing, nothing)) system.boundary_zone_indices .= 1 - - u = system.initial_condition.coordinates v = system.initial_condition.velocity times = collect(tspan[1]:dt:tspan[2]) p_calculated = empty(times) for t in times - v[1, :] .= func(t) - TrixiParticles.calculate_flow_rate_and_pressure!(system, v, u, dt) + system.cache.boundary_zones_flow_rate[1][] = func(t) + TrixiParticles.calculate_pressure!(system, dt) # Store only values after the seventh cycle if t >= 7T