From a489785efb2e742e49a26988df57f796e3b218ae Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 19 Nov 2025 15:10:31 +0100 Subject: [PATCH 01/34] add cross-sectional area --- .../boundary/open_boundary/pressure_model.jl | 18 +++++++++--------- .../boundary/open_boundary/pressure_model.jl | 12 +++++++----- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl index 81ce8f9ff4..95ed7f762d 100644 --- a/src/schemes/boundary/open_boundary/pressure_model.jl +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -19,6 +19,7 @@ It is derived from an electrical circuit analogy and consists of three elements: to store and release volume; in other words, it models the "stretchiness" of the vessel walls. Analogous to a capacitor in an electrical circuit, it absorbs blood when pressure rises and releases it during diastole. The presence of ``C`` smooths pulsatile flow and produces a more uniform outflow profile. +- `cross_sectional_area`: Representative cross-sectional area through which the volumetric flow is evaluated. Lumped-parameter models for the vascular system are well described in the literature (e.g. [Westerhof2008](@cite)). A practical step-by-step procedure for identifying the corresponding model parameters is provided by [Gasser2021](@cite). @@ -32,17 +33,19 @@ struct RCRWindkesselModel{ELTYPE <: Real, P, FR} <: AbstractPressureModel characteristic_resistance :: ELTYPE peripheral_resistance :: ELTYPE compliance :: ELTYPE + cross_sectional_area :: ELTYPE pressure :: P flow_rate :: FR end # The default constructor needs to be accessible for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. -function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance) +function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance, + cross_sectional_area) pressure = Ref(zero(compliance)) flow_rate = Ref(zero(compliance)) return RCRWindkesselModel(characteristic_resistance, peripheral_resistance, compliance, - pressure, flow_rate) + cross_sectional_area, pressure, flow_rate) end function Base.show(io::IO, ::MIME"text/plain", pressure_model::RCRWindkesselModel) @@ -57,6 +60,7 @@ function Base.show(io::IO, ::MIME"text/plain", pressure_model::RCRWindkesselMode summary_line(io, "peripheral_resistance", pressure_model.peripheral_resistance) summary_line(io, "compliance", pressure_model.compliance) + summary_line(io, "cross sectional area", pressure_model.cross_sectional_area) summary_footer(io) end end @@ -90,9 +94,8 @@ end function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system, boundary_zone, v, u, dt) - (; particle_spacing) = system.initial_condition - (; characteristic_resistance, peripheral_resistance, compliance, - flow_rate, pressure) = pressure_model + (; characteristic_resistance, peripheral_resistance, compliance, cross_sectional_area, + flow_rate, pressure) = pressure_model (; face_normal) = boundary_zone # Find particles within the current boundary zone @@ -101,10 +104,7 @@ function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, s 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) - + # the full area of the zone is taken to compute the velocity average. # 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) diff --git a/test/schemes/boundary/open_boundary/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl index f35c2e867f..64f03f2578 100644 --- a/test/schemes/boundary/open_boundary/pressure_model.jl +++ b/test/schemes/boundary/open_boundary/pressure_model.jl @@ -1,7 +1,7 @@ @testset verbose=true "`RCRWindkesselModel`" begin @testset verbose=true "Show" begin pressure_model = RCRWindkesselModel(; peripheral_resistance=1.2, - compliance=0.5, + compliance=0.5, cross_sectional_area=1.2, characteristic_resistance=2.3) show_box = """ @@ -11,6 +11,7 @@ │ characteristic_resistance: ………… 2.3 │ │ peripheral_resistance: …………………… 1.2 │ │ compliance: ………………………………………………… 0.5 │ + │ cross sectional area: ……………………… 1.2 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", pressure_model) == show_box @@ -67,11 +68,12 @@ function simulate_rcr(R1, R2, C, tspan, dt, p_0, func) pressure_model = RCRWindkesselModel(; peripheral_resistance=R2, compliance=C, - characteristic_resistance=R1) - # Define a boundary zone with height=1.0 to ensure a unit volume, - # so velocity directly corresponds to flow rate. + characteristic_resistance=R1, + cross_sectional_area=1.0) reference_velocity(pos, t) = SVector(func(t), 0.0) - boundary_zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), + # Use an arbitrarily wide `boundary_face` to verify that the computed flow rate + # is invariant with respect to the boundary-zone front. + boundary_zone = BoundaryZone(; boundary_face=([0.0, -2.5], [0.0, 3.0]), particle_spacing=0.1, face_normal=(-1.0, 0.0), density=1000.0, reference_velocity, reference_pressure=pressure_model, From 7e51a814bfc6e2d7ea693f96d0724bfa8e0d73aa Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 22 Nov 2025 08:54:34 +0100 Subject: [PATCH 02/34] add ELTYPE --- src/io/read_vtk.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl index 90daa721a3..d87a586bf2 100644 --- a/src/io/read_vtk.jl +++ b/src/io/read_vtk.jl @@ -29,8 +29,10 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu")) │ #particles: ………………………………………………… 100 │ │ particle spacing: ………………………………… 0.1 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘ +``` """ -function vtk2trixi(file) +function vtk2trixi(file; element_type=Float64) + ELTYPE = element_type vtk_file = ReadVTK.VTKFile(file) # Retrieve data fields (e.g., pressure, velocity, ...) @@ -39,7 +41,7 @@ function vtk2trixi(file) # Retrieve fields ndims = first(ReadVTK.get_data(field_data["ndims"])) - coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :] + coordinates = convert.(ELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :]) fields = ["velocity", "density", "pressure", "mass", "particle_spacing"] results = Dict{String, Array{Float64}}() @@ -49,11 +51,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 @@ -62,7 +64,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"], From e0417cbc7c09d043c6c20276ddff7aebbca3a107 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 22 Nov 2025 09:28:07 +0100 Subject: [PATCH 03/34] apply formatter --- src/schemes/boundary/open_boundary/pressure_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl index 95ed7f762d..92a6b575f2 100644 --- a/src/schemes/boundary/open_boundary/pressure_model.jl +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -95,7 +95,7 @@ end function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system, boundary_zone, v, u, dt) (; characteristic_resistance, peripheral_resistance, compliance, cross_sectional_area, - flow_rate, pressure) = pressure_model + flow_rate, pressure) = pressure_model (; face_normal) = boundary_zone # Find particles within the current boundary zone From d0f5aba25a766899b63000eea4aa664d2e5cf5e1 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 26 Nov 2025 16:14:15 +0100 Subject: [PATCH 04/34] fix outside particle --- .../boundary/open_boundary/boundary_zones.jl | 18 +++++++++++++----- src/schemes/boundary/open_boundary/system.jl | 17 +++++++++++++++-- 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index d4031b926a..0b5c8b9acb 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -440,11 +440,19 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi) end end - # This only occurs if `face_normal` is not exactly normal to the `boundary_face`. - # In such cases, particles that are actually outside the simulation domain (outflow particles) - # may be incorrectly kept active as inflow particles and therefore cannot be assigned to any boundary zone. - # This issue should not occur and has been fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926 . - @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle" + # Assert that every active buffer particle is assigned to a boundary zone. + # This should always be true if the boundary zone geometry is set up correctly. + # However, rare edge cases during particle conversion (`convert_particle!`) + # may leave a particle unassigned. Potential causes for failure: + # - `face_normal` is not exactly normal to the `boundary_face`. + # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926.) + # - Large downstream expansion can shift an inflow particle to the zone edge; + # even after upstream adjustment it may remain outside + # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/XXX.) + # - Floating-point rounding when a particle lies almost exactly on the `boundary_face` + # during transition, causing a reset just outside the zone. + # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/XXX.) + @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle" end return system diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4f409716e6..303d4c0049 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -421,9 +421,22 @@ end # 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 + # Reset position of boundary particle back into the boundary zone. + # Particles lying almost exactly on the `boundary_face` might be reset + # slightly outside the boundary zone. + # Thus, we use a shortened vector to account for numerical precision issues. + reset_dist = boundary_zone.zone_width - sqrt(eps(boundary_zone.zone_width)) + reset_vector = -boundary_zone.face_normal * reset_dist for dim in 1:ndims(system) - u[dim, particle] += boundary_zone.spanning_set[1][dim] + u[dim, particle] += reset_vector[dim] + end + + # Verify the particle remains inside the boundary zone after the reset; deactivate it if not. + particle_coords = current_coords(u, system, particle) + if !is_in_boundary_zone(boundary_zone, particle_coords) + deactivate_particle!(system, particle, u) + + return system end impose_rest_density!(v, system, particle, system.boundary_model) From 148e1982caf28cf4e1b3536e6f0a8d1a89c84be2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 26 Nov 2025 21:51:31 +0100 Subject: [PATCH 05/34] add PR link --- src/schemes/boundary/open_boundary/boundary_zones.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 0b5c8b9acb..fd396617ef 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -444,14 +444,14 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi) # This should always be true if the boundary zone geometry is set up correctly. # However, rare edge cases during particle conversion (`convert_particle!`) # may leave a particle unassigned. Potential causes for failure: - # - `face_normal` is not exactly normal to the `boundary_face`. - # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926.) + # - `face_normal` is not exactly normal to the `boundary_face` + # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926). # - Large downstream expansion can shift an inflow particle to the zone edge; # even after upstream adjustment it may remain outside - # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/XXX.) + # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). # - Floating-point rounding when a particle lies almost exactly on the `boundary_face` - # during transition, causing a reset just outside the zone. - # (Fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/XXX.) + # during transition, causing a reset just outside the zone + # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle" end From 0558efb218497a26fcc866af00c16f1b340efa4c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 26 Nov 2025 17:29:10 +0100 Subject: [PATCH 06/34] deactivate lost particle --- src/schemes/boundary/open_boundary/system.jl | 36 +++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4f409716e6..9216f37970 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -82,7 +82,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; mass = copy(initial_conditions.mass) volume = similar(initial_conditions.density) - cache = (; + cache = (; neighbor_counter=zeros(Int, nparticles(initial_conditions)), create_cache_shifting(initial_conditions, shifting_technique)..., create_cache_open_boundary(boundary_model, fluid_system, initial_conditions, boundary_zones_)...) @@ -305,6 +305,8 @@ 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) + deactivate_lost_particles!(system, shifting_technique(system), u, 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) @@ -551,3 +553,35 @@ end return system end + +@inline deactivate_lost_particles!(system, ::Nothing, u, semi) = system + +@inline function deactivate_lost_particles!(system, ::AbstractShiftingTechnique, u, semi) + (; neighbor_counter) = system.cache + + 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 + + @threaded semi for particle in each_integrated_particle(system) + # Deactivate particles that have three or fewer neighbors. + # The threshold of 3 also accounts for possible clustered pairs. + if neighbor_counter[particle] <= 3 + @warn "deactivate particle $particle" + deactivate_particle!(system, particle, u) + end + end + + update_system_buffer!(system.buffer, semi) + + return system +end From 4e4f5a967150f1b87c2d0b0f8f7af19831f60201 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 27 Nov 2025 07:49:22 +0100 Subject: [PATCH 07/34] use ideal_neighbor_count --- src/schemes/boundary/open_boundary/system.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9216f37970..0b924e4a59 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -557,6 +557,7 @@ end @inline deactivate_lost_particles!(system, ::Nothing, u, semi) = system @inline function deactivate_lost_particles!(system, ::AbstractShiftingTechnique, u, semi) + (; particle_spacing) = system.initial_condition (; neighbor_counter) = system.cache set_zero!(neighbor_counter) @@ -572,10 +573,13 @@ end neighbor_counter[particle] += 1 end + # Set 20% of the ideal neighbor count as threshold + min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, + compact_support(system, system)) / 5 + @threaded semi for particle in each_integrated_particle(system) - # Deactivate particles that have three or fewer neighbors. - # The threshold of 3 also accounts for possible clustered pairs. - if neighbor_counter[particle] <= 3 + # Deactivate particles that have too few neighbors + if neighbor_counter[particle] < min_neighbors @warn "deactivate particle $particle" deactivate_particle!(system, particle, u) end From 0826bf784d2c7fd462eed19399dc468b443f7ff6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 27 Nov 2025 10:49:44 +0100 Subject: [PATCH 08/34] reduce threshold --- src/schemes/boundary/open_boundary/system.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 0b924e4a59..59800cce33 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -573,14 +573,13 @@ end neighbor_counter[particle] += 1 end - # Set 20% of the ideal neighbor count as threshold + # Set 10% of the ideal neighbor count as threshold min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, - compact_support(system, system)) / 5 + compact_support(system, system)) / 10 @threaded semi for particle in each_integrated_particle(system) # Deactivate particles that have too few neighbors if neighbor_counter[particle] < min_neighbors - @warn "deactivate particle $particle" deactivate_particle!(system, particle, u) end end From 23981de0c3d3631e1428599110e8d60798191b1b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 27 Nov 2025 14:25:27 +0100 Subject: [PATCH 09/34] make it better --- src/schemes/boundary/open_boundary/system.jl | 33 +++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 59800cce33..d0c793d60f 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -305,7 +305,7 @@ 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) - deactivate_lost_particles!(system, shifting_technique(system), u, semi) + calculate_neighbor_count!(system, shifting_technique(system), u, semi) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) @@ -420,6 +420,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) @@ -554,9 +560,9 @@ end return system end -@inline deactivate_lost_particles!(system, ::Nothing, u, semi) = system +@inline calculate_neighbor_count!(system, ::Nothing, u, semi) = system -@inline function deactivate_lost_particles!(system, ::AbstractShiftingTechnique, u, semi) +@inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi) (; particle_spacing) = system.initial_condition (; neighbor_counter) = system.cache @@ -573,18 +579,21 @@ end neighbor_counter[particle] += 1 end - # Set 10% of the ideal neighbor count as threshold + 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 + min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, compact_support(system, system)) / 10 - @threaded semi for particle in each_integrated_particle(system) - # Deactivate particles that have too few neighbors - if neighbor_counter[particle] < min_neighbors - deactivate_particle!(system, particle, u) - end + if neighbor_counter[particle] < min_neighbors + return true end - update_system_buffer!(system.buffer, semi) - - return system + return false end From e578b68144902de62c9bc6a78dd92efc79b022e2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 28 Nov 2025 09:23:32 +0100 Subject: [PATCH 10/34] deactivate out of bounds particle --- src/general/neighborhood_search.jl | 19 +++++++++++++++++++ src/schemes/boundary/open_boundary/system.jl | 6 ++++++ src/schemes/fluid/fluid.jl | 5 +++++ 3 files changed, 30 insertions(+) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 01397a7732..ea3046de6d 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -9,3 +9,22 @@ 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 + + @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 + + update_system_buffer!(system.buffer, semi) + + return system +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4f409716e6..28b7ea0c8e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -290,6 +290,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) 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 From 280f795602eda14e00dc079af0d47544b9f4861f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 28 Nov 2025 09:57:05 +0100 Subject: [PATCH 11/34] rm variable --- src/schemes/boundary/open_boundary/system.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index d0c793d60f..ca85fa7b5e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -563,7 +563,6 @@ end @inline calculate_neighbor_count!(system, ::Nothing, u, semi) = system @inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi) - (; particle_spacing) = system.initial_condition (; neighbor_counter) = system.cache set_zero!(neighbor_counter) From 0597b8f35d76df7c9000a2471cad4cb7131fa0eb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 28 Nov 2025 11:06:22 +0100 Subject: [PATCH 12/34] check if update is necessary --- src/general/neighborhood_search.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index ea3046de6d..f19cbd0c12 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -24,7 +24,9 @@ function deactivate_out_of_bounds_particles!(system, ::SystemBuffer, end end - update_system_buffer!(system.buffer, semi) + if count(system.buffer.active_particle) != system.buffer.active_particle_count[] + update_system_buffer!(system.buffer, semi) + end return system end From f8ac0cb34d03d9e641907425a48192e201152e5b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 28 Nov 2025 21:15:43 +0100 Subject: [PATCH 13/34] fix --- src/general/neighborhood_search.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index f19cbd0c12..6a58f7c3f9 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -15,11 +15,16 @@ 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) + if !all(min_corner_ .<= particle_position .<= max_corner_) deactivate_particle!(system, particle, u) end end From a51c03e12ff5ac56ed7fa8f81baaa57b5e59f646 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 2 Dec 2025 16:48:52 +0100 Subject: [PATCH 14/34] implement suggestions --- src/schemes/boundary/open_boundary/boundary_zones.jl | 4 ++-- src/schemes/boundary/open_boundary/system.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index fd396617ef..5c288c7dc1 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -446,13 +446,13 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi) # may leave a particle unassigned. Potential causes for failure: # - `face_normal` is not exactly normal to the `boundary_face` # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926). - # - Large downstream expansion can shift an inflow particle to the zone edge; + # - Large downstream domain expansion can shift an inflow particle to the zone edge; # even after upstream adjustment it may remain outside # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). # - Floating-point rounding when a particle lies almost exactly on the `boundary_face` # during transition, causing a reset just outside the zone # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). - @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle" + @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle" end return system diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 303d4c0049..0c8d200748 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -425,7 +425,7 @@ end # Particles lying almost exactly on the `boundary_face` might be reset # slightly outside the boundary zone. # Thus, we use a shortened vector to account for numerical precision issues. - reset_dist = boundary_zone.zone_width - sqrt(eps(boundary_zone.zone_width)) + reset_dist = boundary_zone.zone_width - eps(boundary_zone.zone_width) reset_vector = -boundary_zone.face_normal * reset_dist for dim in 1:ndims(system) u[dim, particle] += reset_vector[dim] From 526d5cc1f880590f93135d454a5e890594a3380f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 3 Dec 2025 09:37:16 +0100 Subject: [PATCH 15/34] add flag to deactivate --- src/schemes/boundary/open_boundary/system.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index ca85fa7b5e..8dc5ae32cb 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -64,6 +64,7 @@ end function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; fluid_system::AbstractFluidSystem, buffer_size::Integer, boundary_model, + deactivate_isolated_particles=false, pressure_acceleration=boundary_model isa BoundaryModelDynamicalPressureZhang ? fluid_system.pressure_acceleration_formulation : @@ -82,7 +83,9 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; mass = copy(initial_conditions.mass) volume = similar(initial_conditions.density) - cache = (; neighbor_counter=zeros(Int, nparticles(initial_conditions)), + 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_)...) @@ -565,6 +568,9 @@ end @inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi) (; neighbor_counter) = system.cache + # Skip when `deactivate_isolated_particles` ist not activated + isnothing(neighbor_counter) && return system + set_zero!(neighbor_counter) # Loop over all pairs of particles and neighbors within the kernel cutoff. @@ -587,6 +593,9 @@ end (; particle_spacing) = system.initial_condition (; neighbor_counter) = system.cache + # Skip when `deactivate_isolated_particles` ist not activated + isnothing(neighbor_counter) && return false + min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, compact_support(system, system)) / 10 From 43e5068891aa596083f4bcb6c656ef903af62013 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 3 Dec 2025 11:15:35 +0100 Subject: [PATCH 16/34] fix typo --- src/schemes/boundary/open_boundary/system.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 8dc5ae32cb..06bd0b39c3 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -568,7 +568,7 @@ end @inline function calculate_neighbor_count!(system, ::AbstractShiftingTechnique, u, semi) (; neighbor_counter) = system.cache - # Skip when `deactivate_isolated_particles` ist not activated + # Skip when `deactivate_isolated_particles` is not activated isnothing(neighbor_counter) && return system set_zero!(neighbor_counter) @@ -593,7 +593,7 @@ end (; particle_spacing) = system.initial_condition (; neighbor_counter) = system.cache - # Skip when `deactivate_isolated_particles` ist not activated + # Skip when `deactivate_isolated_particles` is not activated isnothing(neighbor_counter) && return false min_neighbors = ideal_neighbor_count(Val(ndims(system)), particle_spacing, From 63301183b959e914d44b8cdef95670991092cc79 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Dec 2025 08:12:22 +0100 Subject: [PATCH 17/34] ramp up shifting velocity --- src/schemes/boundary/open_boundary/system.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4f409716e6..332382adac 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -538,13 +538,16 @@ end dist_free_surface = boundary_zone.zone_width - dist_to_transition if dist_free_surface < compact_support(fluid_system, fluid_system) - # Disable shifting for this particle. + # Ramp up shifting velocity for this particle. # Note that Sun et al. 2017 propose a more sophisticated approach with a transition phase # where only the component orthogonal to the surface normal is kept and the tangential - # component is set to zero. However, we assume laminar flow in the boundary zone, - # so we simply disable shifting completely. + # component is set to zero. + kernel_max = smoothing_kernel(system, 0, particle) + support_gap = compact_support(fluid_system, fluid_system) - dist_free_surface + shifting_weight = smoothing_kernel(system, support_gap, particle) / kernel_max + delta_v_ramped = delta_v(system, particle) * shifting_weight for dim in 1:ndims(system) - cache.delta_v[dim, particle] = zero(eltype(system)) + cache.delta_v[dim, particle] = delta_v_ramped[dim] end end end From 83192220896fb7c56dfec5ffbd7e8e4cda48ec14 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Dec 2025 08:26:17 +0100 Subject: [PATCH 18/34] add comment --- src/schemes/boundary/open_boundary/system.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 332382adac..8434538c49 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -541,7 +541,8 @@ end # Ramp up shifting velocity for this particle. # Note that Sun et al. 2017 propose a more sophisticated approach with a transition phase # where only the component orthogonal to the surface normal is kept and the tangential - # component is set to zero. + # component is set to zero. However, our experiments showed that this approach + # leads to stronger particle disorder in the boundary zone. kernel_max = smoothing_kernel(system, 0, particle) support_gap = compact_support(fluid_system, fluid_system) - dist_free_surface shifting_weight = smoothing_kernel(system, support_gap, particle) / kernel_max From 2f38bc18f226dfe39194dc5995022166c08a94fb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Dec 2025 14:38:08 +0100 Subject: [PATCH 19/34] update comment --- src/schemes/boundary/open_boundary/system.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 8434538c49..96fa441b13 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -538,11 +538,10 @@ end dist_free_surface = boundary_zone.zone_width - dist_to_transition if dist_free_surface < compact_support(fluid_system, fluid_system) - # Ramp up shifting velocity for this particle. - # Note that Sun et al. 2017 propose a more sophisticated approach with a transition phase - # where only the component orthogonal to the surface normal is kept and the tangential - # component is set to zero. However, our experiments showed that this approach - # leads to stronger particle disorder in the boundary zone. + # Ramp shifting velocity near the free surface using a kernel-weighted transition. + # According to our experiment, alternative approaches lead to particle disorder anyway: + # - Sun et al. 2017: project onto surface-tangential component + # - Zhang et al. 2025: disable shifting entirely kernel_max = smoothing_kernel(system, 0, particle) support_gap = compact_support(fluid_system, fluid_system) - dist_free_surface shifting_weight = smoothing_kernel(system, support_gap, particle) / kernel_max From 398dd3dd7b671231dff07e178d472624770d4ebc Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Dec 2025 16:44:57 +0100 Subject: [PATCH 20/34] implement suggestions --- src/schemes/boundary/open_boundary/system.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 96fa441b13..4a09b28abe 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -539,12 +539,14 @@ end if dist_free_surface < compact_support(fluid_system, fluid_system) # Ramp shifting velocity near the free surface using a kernel-weighted transition. - # According to our experiment, alternative approaches lead to particle disorder anyway: - # - Sun et al. 2017: project onto surface-tangential component + # According to our experiment, alternative approaches lead to particle disorder: + # - Sun et al. 2017: only use surface-tangential component # - Zhang et al. 2025: disable shifting entirely kernel_max = smoothing_kernel(system, 0, particle) - support_gap = compact_support(fluid_system, fluid_system) - dist_free_surface - shifting_weight = smoothing_kernel(system, support_gap, particle) / kernel_max + dist_from_cutoff = compact_support(fluid_system, fluid_system) - + dist_free_surface + shifting_weight = smoothing_kernel(system, dist_from_cutoff, particle) / + kernel_max delta_v_ramped = delta_v(system, particle) * shifting_weight for dim in 1:ndims(system) cache.delta_v[dim, particle] = delta_v_ramped[dim] From 8f40241e8e3038ae9ae3ece2fa473be5570d8b8c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 6 Dec 2025 08:25:33 +0100 Subject: [PATCH 21/34] modify comment --- src/schemes/boundary/open_boundary/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4a09b28abe..4abdf15fb5 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -539,7 +539,7 @@ end if dist_free_surface < compact_support(fluid_system, fluid_system) # Ramp shifting velocity near the free surface using a kernel-weighted transition. - # According to our experiment, alternative approaches lead to particle disorder: + # According to our experiments, the proposed alternative approaches lead to particle disorder: # - Sun et al. 2017: only use surface-tangential component # - Zhang et al. 2025: disable shifting entirely kernel_max = smoothing_kernel(system, 0, particle) From 76efc2db3b4e6d66f5f908a02e2e4c481d3e7672 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 8 Dec 2025 09:01:19 +0100 Subject: [PATCH 22/34] poc --- .../boundary/open_boundary/pressure_model.jl | 16 ++++++++-------- .../boundary/open_boundary/pressure_model.jl | 6 ++---- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl index 92a6b575f2..561f3e91ef 100644 --- a/src/schemes/boundary/open_boundary/pressure_model.jl +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -103,16 +103,16 @@ function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, s 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 to compute the velocity average. - # 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) + # Compute volumetric flow rate: Q = ∫ v ⋅ n dA + current_flow_rate = sum(candidates) do particle + v_n = dot(current_velocity(v, system, particle), -face_normal) + V_particle = hydrodynamic_mass(system, particle) / + current_density(v, system, particle) + A_particle = V_particle / boundary_zone.zone_width + + return v_n * A_particle end - # Compute volumetric flow rate: Q = v * A - current_flow_rate = velocity_avg * cross_sectional_area - previous_pressure = pressure[] previous_flow_rate = flow_rate[] flow_rate[] = current_flow_rate diff --git a/test/schemes/boundary/open_boundary/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl index 64f03f2578..c33b567e09 100644 --- a/test/schemes/boundary/open_boundary/pressure_model.jl +++ b/test/schemes/boundary/open_boundary/pressure_model.jl @@ -71,13 +71,11 @@ characteristic_resistance=R1, cross_sectional_area=1.0) reference_velocity(pos, t) = SVector(func(t), 0.0) - # Use an arbitrarily wide `boundary_face` to verify that the computed flow rate - # is invariant with respect to the boundary-zone front. - boundary_zone = BoundaryZone(; boundary_face=([0.0, -2.5], [0.0, 3.0]), + boundary_zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), 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=7, rest_pressure=p_0) system = OpenBoundarySystem(boundary_zone; buffer_size=0, boundary_model=nothing, From 5990c42cb7a4f2507a003bd709771496ae310802 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 8 Dec 2025 14:50:58 +0100 Subject: [PATCH 23/34] remove cross sectional area kwarg --- .../boundary/open_boundary/pressure_model.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl index 561f3e91ef..8e81d6e8ec 100644 --- a/src/schemes/boundary/open_boundary/pressure_model.jl +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -19,7 +19,6 @@ It is derived from an electrical circuit analogy and consists of three elements: to store and release volume; in other words, it models the "stretchiness" of the vessel walls. Analogous to a capacitor in an electrical circuit, it absorbs blood when pressure rises and releases it during diastole. The presence of ``C`` smooths pulsatile flow and produces a more uniform outflow profile. -- `cross_sectional_area`: Representative cross-sectional area through which the volumetric flow is evaluated. Lumped-parameter models for the vascular system are well described in the literature (e.g. [Westerhof2008](@cite)). A practical step-by-step procedure for identifying the corresponding model parameters is provided by [Gasser2021](@cite). @@ -33,19 +32,17 @@ struct RCRWindkesselModel{ELTYPE <: Real, P, FR} <: AbstractPressureModel characteristic_resistance :: ELTYPE peripheral_resistance :: ELTYPE compliance :: ELTYPE - cross_sectional_area :: ELTYPE pressure :: P flow_rate :: FR end # The default constructor needs to be accessible for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. -function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance, - cross_sectional_area) +function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance) pressure = Ref(zero(compliance)) flow_rate = Ref(zero(compliance)) return RCRWindkesselModel(characteristic_resistance, peripheral_resistance, compliance, - cross_sectional_area, pressure, flow_rate) + pressure, flow_rate) end function Base.show(io::IO, ::MIME"text/plain", pressure_model::RCRWindkesselModel) @@ -94,7 +91,7 @@ end function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system, boundary_zone, v, u, dt) - (; characteristic_resistance, peripheral_resistance, compliance, cross_sectional_area, + (; characteristic_resistance, peripheral_resistance, compliance, flow_rate, pressure) = pressure_model (; face_normal) = boundary_zone @@ -105,12 +102,13 @@ function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, s # Compute volumetric flow rate: Q = ∫ v ⋅ n dA current_flow_rate = sum(candidates) do particle - v_n = dot(current_velocity(v, system, particle), -face_normal) + vn = dot(current_velocity(v, system, particle), -face_normal) V_particle = hydrodynamic_mass(system, particle) / current_density(v, system, particle) + # Cross-sectional area represented by the particle in the boundary zone A_particle = V_particle / boundary_zone.zone_width - return v_n * A_particle + return vn * A_particle end previous_pressure = pressure[] From b043bd11e4cc43cebb23bdba1388f5174f59a653 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 8 Dec 2025 17:52:26 +0100 Subject: [PATCH 24/34] add `coordinates_eltype` --- src/io/read_vtk.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl index 326773d4b4..884588033b 100644 --- a/src/io/read_vtk.jl +++ b/src/io/read_vtk.jl @@ -33,8 +33,9 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu")) └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -function vtk2trixi(file; element_type=Float64) +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, ...) @@ -43,7 +44,7 @@ function vtk2trixi(file; element_type=Float64) # Retrieve fields ndims = first(ReadVTK.get_data(field_data["ndims"])) - coordinates = convert.(ELTYPE, 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}}() From 0359425131784b5ebc3064a15fdc03b2613d615f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 09:30:19 +0100 Subject: [PATCH 25/34] first prototype --- .../boundary/open_boundary/boundary_zones.jl | 23 ++++- src/schemes/boundary/open_boundary/system.jl | 99 ++++++++++++++++--- 2 files changed, 102 insertions(+), 20 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 7088fea58e..be5979ba6e 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -147,7 +147,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 +156,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 +168,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=nothing, rest_pressure=zero(eltype(density)), reference_density=nothing, reference_pressure=nothing, reference_velocity=nothing) @@ -262,10 +263,12 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, ic.velocity .= stack(velocity_ref.(coordinates_svector, 0)) end + cache = (; create_cache_boundary_zone(ic, 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) @@ -301,6 +304,16 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone) end end +create_cache_boundary_zone(initial_condition, sample_points::Nothing) = (;) + +function create_cache_boundary_zone(initial_condition, sample_points::Matrix) + # TODO: Check matrix (ndims etc.) + shepard_coefficient = zeros(eltype(initial_condition), axes(sample_points, 2)) + dA = initial_condition.particle_spacing + return (; sample_points=sample_points, sample_velocity=copy(sample_points), + shepard_coefficient, dA) +end + function set_up_boundary_zone(boundary_face, face_normal, density, particle_spacing, initial_condition, extrude_geometry, open_boundary_layers, boundary_type) @@ -453,7 +466,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi) # - Floating-point rounding when a particle lies almost exactly on the `boundary_face` # during transition, causing a reset just outside the zone # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). - @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle" + @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle" end return system diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9396c3f1fb..2c4c5c6c0c 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -63,7 +63,7 @@ end function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; fluid_system::AbstractFluidSystem, buffer_size::Integer, - boundary_model, + boundary_model, calculate_flow_rate=false, pressure_acceleration=boundary_model isa BoundaryModelDynamicalPressureZhang ? fluid_system.pressure_acceleration_formulation : @@ -85,7 +85,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; cache = (; 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 +110,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 +133,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 +142,15 @@ 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) + 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 + boundary_zones_flow_rate = zeros(ELTYPE, 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 +158,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 +179,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 +191,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 @@ -305,6 +307,10 @@ 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) + # 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, u, 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 +329,28 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system +function calculate_flow_rate!(system, v, u, v_ode, u_ode, semi) + system.cache.calculate_flow_rate || return system + + for boundary_zone in system.boundary_zones + interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi) + end + + foreach_enumerate(system.boundary_zones) do (zone_id, boundary_zone) + (; face_normal) = boundary_zone + (; sample_velocity, dA) = boundary_zone.cache + # Compute volumetric flow rate: Q = ∫ v ⋅ n dA + current_flow_rate = sum(axes(sample_velocity, 2)) do point + vn = dot(current_velocity(sample_velocity, system, point), -face_normal) + return vn * dA + end + + system.cache.boundary_zones_flow_rate[zone_id] = current_flow_rate + end + + return system +end + function check_domain!(system, v, u, v_ode, u_ode, semi) (; boundary_zones, boundary_candidates, fluid_candidates, fluid_system) = system @@ -570,3 +598,44 @@ end return system end + +function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, v, u, + 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 From fd15671fcb56a9a37a276d75a396457b10b94cd5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 11:03:16 +0100 Subject: [PATCH 26/34] write Q --- src/io/write_vtk.jl | 8 ++++++++ src/schemes/boundary/open_boundary/boundary_zones.jl | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 07fa559d2e..d1146e67b1 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -409,6 +409,14 @@ 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 + for i in eachindex(system.cache.boundary_zones_flow_rate) + vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i] + end + + vtk["Q_total"] = sum(system.cache.boundary_zones_flow_rate) + 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 be5979ba6e..2bde137093 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -466,7 +466,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi) # - Floating-point rounding when a particle lies almost exactly on the `boundary_face` # during transition, causing a reset just outside the zone # (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997). - @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle" + @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle" end return system From 364e24ea6894836aee569d2e11d7a934c0a7642b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 11:42:06 +0100 Subject: [PATCH 27/34] gpu fix --- src/io/write_vtk.jl | 6 ++++-- src/schemes/boundary/open_boundary/system.jl | 20 +++++++++++++------- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index d1146e67b1..ff2f2f85eb 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -410,11 +410,13 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem) 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] + vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i][] + Q_total += system.cache.boundary_zones_flow_rate[i][] end - vtk["Q_total"] = sum(system.cache.boundary_zones_flow_rate) + vtk["Q_total"] = Q_total end return vtk diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 2c4c5c6c0c..12cde7de5e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -147,7 +147,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit velocity_reference_values=velocity_reference_values, calculate_flow_rate) if calculate_flow_rate - boundary_zones_flow_rate = zeros(ELTYPE, length(boundary_zones)) + boundary_zones_flow_rate = ntuple(i -> Ref(zero(ELTYPE)), + Val(length(boundary_zones))) cache = (; boundary_zones_flow_rate, cache...) end @@ -329,23 +330,28 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system -function calculate_flow_rate!(system, v, u, v_ode, u_ode, semi) - system.cache.calculate_flow_rate || return system +function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode, + u_ode, semi) where {ELTYPE, NDIMS} + (; calculate_flow_rate, boundary_zones_flow_rate) = system.cache + calculate_flow_rate || return system for boundary_zone in system.boundary_zones interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi) end - foreach_enumerate(system.boundary_zones) do (zone_id, boundary_zone) + foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate) + boundary_zone = system.boundary_zones[zone_id] (; face_normal) = boundary_zone (; sample_velocity, dA) = boundary_zone.cache + # Compute volumetric flow rate: Q = ∫ v ⋅ n dA - current_flow_rate = sum(axes(sample_velocity, 2)) do point - vn = dot(current_velocity(sample_velocity, system, point), -face_normal) + velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity) + current_flow_rate = sum(velocities) do velocity + vn = dot(velocity, -face_normal) return vn * dA end - system.cache.boundary_zones_flow_rate[zone_id] = current_flow_rate + boundary_zone_flow_rate[] = current_flow_rate end return system From c1b0abf6536801fa6669908a5f32030b31d693f5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 14:46:11 +0100 Subject: [PATCH 28/34] add checks --- .../boundary/open_boundary/boundary_zones.jl | 54 ++++++++++++++++--- src/schemes/boundary/open_boundary/system.jl | 14 +++-- 2 files changed, 56 insertions(+), 12 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 2bde137093..bfefefca4c 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -263,7 +263,8 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, ic.velocity .= stack(velocity_ref.(coordinates_svector, 0)) end - cache = (; create_cache_boundary_zone(ic, sample_points)...) + 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, @@ -300,18 +301,55 @@ 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 -create_cache_boundary_zone(initial_condition, sample_points::Nothing) = (;) +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 + @warn "The sampled area of the boundary face " * + "($(discrete_face_area)) is larger than the actual face area ($(face_area)). " + end -function create_cache_boundary_zone(initial_condition, sample_points::Matrix) - # TODO: Check matrix (ndims etc.) - shepard_coefficient = zeros(eltype(initial_condition), axes(sample_points, 2)) - dA = initial_condition.particle_spacing - return (; sample_points=sample_points, sample_velocity=copy(sample_points), - shepard_coefficient, dA) + 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, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 12cde7de5e..4bed038600 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -147,6 +147,12 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit 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...) @@ -332,8 +338,8 @@ update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = syst function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode, u_ode, semi) where {ELTYPE, NDIMS} - (; calculate_flow_rate, boundary_zones_flow_rate) = system.cache - calculate_flow_rate || return system + system.cache.calculate_flow_rate || return system + (; boundary_zones_flow_rate) = system.cache for boundary_zone in system.boundary_zones interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi) @@ -342,13 +348,13 @@ function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate) boundary_zone = system.boundary_zones[zone_id] (; face_normal) = boundary_zone - (; sample_velocity, dA) = boundary_zone.cache + (; sample_velocity, area_increment) = boundary_zone.cache # 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 * dA + return vn * area_increment end boundary_zone_flow_rate[] = current_flow_rate From eeebcc22a371b01d4d51d50ddf90ad31c3ae7a40 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 15:14:33 +0100 Subject: [PATCH 29/34] add docs --- src/schemes/boundary/open_boundary/boundary_zones.jl | 6 +++++- src/schemes/boundary/open_boundary/system.jl | 4 +++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index bfefefca4c..043beb8196 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`) @@ -168,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(), sample_points=nothing, + boundary_type=BidirectionalFlow(), sample_points=:default, rest_pressure=zero(eltype(density)), reference_density=nothing, reference_pressure=nothing, reference_velocity=nothing) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 4bed038600..36889bc0e2 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). From d21a5492311db0c13e75aa7fa77866cef49814d5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 10 Dec 2025 17:29:05 +0100 Subject: [PATCH 30/34] adapt pressure model --- .../boundary/open_boundary/pressure_model.jl | 39 +++++-------------- src/schemes/boundary/open_boundary/system.jl | 4 ++ src/util.jl | 14 +++++++ .../boundary/open_boundary/pressure_model.jl | 6 +-- 4 files changed, 30 insertions(+), 33 deletions(-) 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 36889bc0e2..7038f48534 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -144,6 +144,10 @@ 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) diff --git a/src/util.jl b/src/util.jl index 18389cd44b..be032a8b95 100644 --- a/src/util.jl +++ b/src/util.jl @@ -11,6 +11,20 @@ end @inline foreach_noalloc(func, collection::Tuple{}) = nothing +@inline function foreach_noalloc(func, collection1, collection2) + element1 = first(collection1) + remaining_collection1 = Base.tail(collection1) + element2 = first(collection2) + remaining_collection2 = Base.tail(collection2) + + func((element1, element2)) + + # Process remaining collection + foreach_noalloc(func, remaining_collection1, remaining_collection2) +end + +@inline foreach_noalloc(func, collection1::Tuple{}, collection2::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) diff --git a/test/schemes/boundary/open_boundary/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl index f35c2e867f..c0e1baffe7 100644 --- a/test/schemes/boundary/open_boundary/pressure_model.jl +++ b/test/schemes/boundary/open_boundary/pressure_model.jl @@ -81,15 +81,13 @@ 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 From b4020c8df060d51bb3f012ef4c977aabe14c7834 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 11 Dec 2025 08:07:32 +0100 Subject: [PATCH 31/34] use `foreach_no_alloc` --- examples/fluid/poiseuille_flow_2d.jl | 3 +- src/schemes/boundary/open_boundary/system.jl | 40 ++++++++++---------- src/util.jl | 14 +++++++ 3 files changed, 37 insertions(+), 20 deletions(-) diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl index dd324a5bca..1afacd81a8 100644 --- a/examples/fluid/poiseuille_flow_2d.jl +++ b/examples/fluid/poiseuille_flow_2d.jl @@ -29,7 +29,7 @@ open_boundary_layers = 10 # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 2.0) +tspan = (0.0, 1.0) wcsph = true domain_size = (flow_length, wall_distance) @@ -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/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 36889bc0e2..648a97cb03 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -318,7 +318,7 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode # 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, u, v_ode, u_ode, semi) + calculate_flow_rate!(system, v_ode, u_ode, semi) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) @@ -338,28 +338,30 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system -function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode, - u_ode, semi) where {ELTYPE, NDIMS} +function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, + v_ode, u_ode, semi) where {ELTYPE, NDIMS} system.cache.calculate_flow_rate || return system - (; boundary_zones_flow_rate) = system.cache - for boundary_zone in system.boundary_zones - interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi) - end + @trixi_timeit timer() "flow rate calculation" begin + (; boundary_zones) = system + (; boundary_zones_flow_rate) = system.cache - foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate) - boundary_zone = system.boundary_zones[zone_id] - (; face_normal) = boundary_zone - (; sample_velocity, area_increment) = boundary_zone.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 - # 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 + interpolate_velocity!(system, boundary_zone, v_ode, u_ode, semi) - boundary_zone_flow_rate[] = current_flow_rate + # 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 @@ -613,7 +615,7 @@ end return system end -function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, v, u, +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) diff --git a/src/util.jl b/src/util.jl index 18389cd44b..be032a8b95 100644 --- a/src/util.jl +++ b/src/util.jl @@ -11,6 +11,20 @@ end @inline foreach_noalloc(func, collection::Tuple{}) = nothing +@inline function foreach_noalloc(func, collection1, collection2) + element1 = first(collection1) + remaining_collection1 = Base.tail(collection1) + element2 = first(collection2) + remaining_collection2 = Base.tail(collection2) + + func((element1, element2)) + + # Process remaining collection + foreach_noalloc(func, remaining_collection1, remaining_collection2) +end + +@inline foreach_noalloc(func, collection1::Tuple{}, collection2::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) From 55be6a601973cc7a9eff94cafe9d6b0a5114feaf Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 11 Dec 2025 08:28:07 +0100 Subject: [PATCH 32/34] fix unit tests --- .../boundary/open_boundary/boundary_zones.jl | 2 +- .../boundary/open_boundary/boundary_zone.jl | 62 ++++++++++++++----- 2 files changed, 49 insertions(+), 15 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 043beb8196..55b86be737 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma face_area = norm(v2 - v1) end - if discrete_face_area > face_area + 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 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 From a431873732519aaae7b4d748d612e3cea2b5e7be Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 11 Dec 2025 08:28:07 +0100 Subject: [PATCH 33/34] fix unit tests --- examples/fluid/poiseuille_flow_2d.jl | 2 +- .../boundary/open_boundary/boundary_zones.jl | 2 +- .../boundary/open_boundary/boundary_zone.jl | 62 ++++++++++++++----- 3 files changed, 50 insertions(+), 16 deletions(-) diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl index 1afacd81a8..5e57c7c0bd 100644 --- a/examples/fluid/poiseuille_flow_2d.jl +++ b/examples/fluid/poiseuille_flow_2d.jl @@ -29,7 +29,7 @@ open_boundary_layers = 10 # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 1.0) +tspan = (0.0, 2.0) wcsph = true domain_size = (flow_length, wall_distance) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 043beb8196..55b86be737 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma face_area = norm(v2 - v1) end - if discrete_face_area > face_area + 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 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 From 9339bcbe7a217f84e59a819ca91b17e8fbefb647 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 11 Dec 2025 08:28:07 +0100 Subject: [PATCH 34/34] fix unit tests --- examples/fluid/poiseuille_flow_2d.jl | 2 +- .../boundary/open_boundary/boundary_zones.jl | 2 +- .../boundary/open_boundary/boundary_zone.jl | 46 ++++++++++++++++++- 3 files changed, 46 insertions(+), 4 deletions(-) diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl index 1afacd81a8..5e57c7c0bd 100644 --- a/examples/fluid/poiseuille_flow_2d.jl +++ b/examples/fluid/poiseuille_flow_2d.jl @@ -29,7 +29,7 @@ open_boundary_layers = 10 # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 1.0) +tspan = (0.0, 2.0) wcsph = true domain_size = (flow_length, wall_distance) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 043beb8196..55b86be737 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma face_area = norm(v2 - v1) end - if discrete_face_area > face_area + 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 diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index 619c048bbe..8542fd554d 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