From 6882e58c85234eff9315bb7cee18273a5cc9c6e4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 17 Jul 2025 22:21:06 +0200 Subject: [PATCH] first working prototype --- src/TrixiParticles.jl | 2 +- src/general/semidiscretization.jl | 24 ++- src/schemes/boundary/boundary.jl | 5 +- .../dummy_particles/dummy_particles.jl | 17 +- .../boundary/open_boundary/boundary_zones.jl | 4 +- .../open_boundary/dynamic_pressure.jl | 199 ++++++++++++++++++ src/schemes/boundary/open_boundary/system.jl | 74 ++++++- 7 files changed, 303 insertions(+), 22 deletions(-) create mode 100644 src/schemes/boundary/open_boundary/dynamic_pressure.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 5e79798f35..f7917d70a5 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -80,7 +80,7 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr export tensile_instability_control export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni, - BernoulliPressureExtrapolation + BoundaryModelZhang, BernoulliPressureExtrapolation export HertzContactModel, LinearContactModel export BoundaryMovement export examples_dir, validation_dir diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 3e033e447b..3115d12a36 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -154,6 +154,12 @@ end return 0.0 end +@inline function compact_support(system::OpenBoundarySPHSystem{<:BoundaryModelZhang}, + neighbor::OpenBoundarySPHSystem{<:BoundaryModelZhang}) + # Use the compact support of the fluid + return compact_support(system.fluid_system, neighbor.fluid_system) +end + @inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem) # This NHS is never used return 0.0 @@ -683,7 +689,9 @@ function update_nhs!(neighborhood_search, end function update_nhs!(neighborhood_search, - system::FluidSystem, neighbor::BoundarySPHSystem, + system::Union{FluidSystem, + OpenBoundarySPHSystem{<:BoundaryModelZhang}}, + neighbor::BoundarySPHSystem, u_system, u_neighbor, semi) # Boundary coordinates only change over time when `neighbor.ismoving[]` update!(neighborhood_search, @@ -718,6 +726,20 @@ function update_nhs!(neighborhood_search, semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) end +function update_nhs!(neighborhood_search, + system::OpenBoundarySPHSystem{<:BoundaryModelZhang}, + neighbor::OpenBoundarySPHSystem{<:BoundaryModelZhang}, + u_system, u_neighbor, semi) + # The current coordinates of both open boundaries and fluids change over time. + + # TODO: Update only `active_coordinates` of open boundaries. + # Problem: Removing inactive particles from neighboring lists is necessary. + update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) +end + function update_nhs!(neighborhood_search, system::OpenBoundarySPHSystem, neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi) diff --git a/src/schemes/boundary/boundary.jl b/src/schemes/boundary/boundary.jl index a09545a9b7..6ddd7250c4 100644 --- a/src/schemes/boundary/boundary.jl +++ b/src/schemes/boundary/boundary.jl @@ -1,9 +1,10 @@ -include("dummy_particles/dummy_particles.jl") -include("system.jl") include("open_boundary/boundary_zones.jl") include("open_boundary/mirroring.jl") include("open_boundary/method_of_characteristics.jl") include("open_boundary/system.jl") +include("open_boundary/dynamic_pressure.jl") +include("dummy_particles/dummy_particles.jl") +include("system.jl") # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem` # and the `TotalLagrangianSPHSystem` and are therefore included later. diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 9912b91726..5d78fd29ad 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -470,7 +470,9 @@ end end @inline function boundary_pressure_extrapolation!(parallel::Val{true}, boundary_model, - system, neighbor_system::FluidSystem, + system, + neighbor_system::Union{FluidSystem, + OpenBoundarySPHSystem{<:BoundaryModelZhang}}, system_coords, neighbor_coords, v, v_neighbor_system, semi) (; pressure, cache, viscosity, density_calculator) = boundary_model @@ -492,7 +494,9 @@ end # Note that this needs to be serial, as we are writing into the same # pressure entry from different loop iterations. @inline function boundary_pressure_extrapolation!(parallel::Val{false}, boundary_model, - system, neighbor_system::FluidSystem, + system, + neighbor_system::Union{FluidSystem, + OpenBoundarySPHSystem{<:BoundaryModelZhang}}, system_coords, neighbor_coords, v, v_neighbor_system, semi) (; pressure, cache, viscosity, density_calculator) = boundary_model @@ -513,7 +517,10 @@ end end @inline function boundary_pressure_inner!(boundary_model, boundary_density_calculator, - system, neighbor_system::FluidSystem, v, + system, + neighbor_system::Union{FluidSystem, + OpenBoundarySPHSystem{<:BoundaryModelZhang}}, + v, v_neighbor_system, particle, neighbor, pos_diff, distance, viscosity, cache, pressure, pressure_offset) @@ -525,8 +532,8 @@ end neighbor) # Hydrostatic pressure term from fluid and boundary acceleration - resulting_acceleration = neighbor_system.acceleration - - @inbounds current_acceleration(system, particle) + resulting_acceleration = @inbounds current_acceleration(system, particle) - + @inbounds current_acceleration(system, particle) hydrostatic_pressure = dot(resulting_acceleration, density_neighbor * pos_diff) # Additional dynamic pressure term (only with `BernoulliPressureExtrapolation`) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 30e589e1c9..97d4a254f8 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -262,7 +262,7 @@ end end function remove_outside_particles(initial_condition, spanning_set, zone_origin) - (; coordinates, density, particle_spacing) = initial_condition + (; coordinates, velocity, density, particle_spacing) = initial_condition in_zone = fill(true, nparticles(initial_condition)) @@ -274,5 +274,5 @@ function remove_outside_particles(initial_condition, spanning_set, zone_origin) end return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), - particle_spacing) + velocity=velocity[:, in_zone], particle_spacing) end diff --git a/src/schemes/boundary/open_boundary/dynamic_pressure.jl b/src/schemes/boundary/open_boundary/dynamic_pressure.jl new file mode 100644 index 0000000000..d06b0a81e5 --- /dev/null +++ b/src/schemes/boundary/open_boundary/dynamic_pressure.jl @@ -0,0 +1,199 @@ +struct BoundaryModelZhang{EOS, BP} + state_equation :: EOS + boundary_pressure :: BP +end + +function BoundaryModelZhang(; state_equation, boundary_pressure) + + # Theoretically, spatially non-uniform boundary pressure could also be + # imposed since Eq. (12) does not require pb to be constant along the + # channel cross section, although it is generally not necessary. + return BoundaryModelZhang(state_equation, boundary_pressure) +end + +@inline function v_nvariables(system, boundary_model::BoundaryModelZhang) + return ndims(system) + 1 +end + +@inline function current_density(v, ::BoundaryModelZhang, system::OpenBoundarySPHSystem) + # When using `ContinuityDensity`, the density is stored in the last row of `v` + return view(v, size(v, 1), :) +end + +function write_v0!(v0, system::OpenBoundarySPHSystem, ::BoundaryModelZhang) + (; reference_density, initial_condition) = system + + coords_svector = reinterpret(reshape, SVector{ndims(system), eltype(system)}, + initial_condition.coordinates) + + # Note that `.=` is very slightly faster, but not GPU-compatible + v0[end, + :] = copy(reference_value.(reference_density, initial_condition.density, + coords_svector, 0)) + + return v0 +end + +@inline impose_new_density!(v, u, system, boundary_model, particle, t) = v + +function impose_new_density!(v, u, system, boundary_model::BoundaryModelZhang, particle, t) + (; boundary_pressure, boundary_pressure_function, boundary_pressure) = system.cache + (; state_equation) = boundary_model + + density = current_density(v, boundary_model, system) + + particle_coords = current_coords(u, system, particle) + p_current = boundary_pressure[particle] + + p_boundary = reference_value(boundary_pressure_function, p_current, particle_coords, t) + + # the density of the newly populated (actually recycled) particles in the bidirectional + # in-/outflow buffer is obtained following the boundary pressure and EoS + @inbounds density[particle] = inverse_state_equation(state_equation, p_boundary) + + return v +end + +function update_final!(system, boundary_model::BoundaryModelZhang, v, u, v_ode, u_ode, + semi, t) + (; boundary_pressure, boundary_pressure_function, prescribed_pressure) = system.cache + + @threaded semi for particle in eachparticle(system) + particle_coords = current_coords(u, system, particle) + p_current = boundary_pressure[particle] + + boundary_pressure[particle] = reference_value(boundary_pressure_function, p_current, + particle_coords, t) + + if !(prescribed_pressure) + rho = current_density(v, system, particle) + system.pressure[particle] = boundary_model.state_equation(rho) + end + end + + return system +end + +function update_boundary_quantities!(system, boundary_model::BoundaryModelZhang, v, u, + v_ode, u_ode, semi, t) + (; cache, pressure, boundary_zone, reference_density, + reference_velocity, reference_pressure) = system + (; prescribed_density, prescribed_pressure, prescribed_velocity) = cache + + density = current_density(v, boundary_model, system) + + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + if prescribed_pressure + p_current = current_pressure(v, system, particle) + pressure[particle] = reference_value(reference_pressure, p_current, + particle_coords, t) + end + + if prescribed_velocity + v_current = current_velocity(v, system, particle) + v_ref = reference_value(reference_velocity, v_current, particle_coords, t) + @inbounds for dim in eachindex(v_ref) + v[dim, particle] = v_ref[dim] + end + end + + if prescribed_density + rho_current = current_density(v, system, particle) + density[particle] = reference_value(reference_density, rho_current, + particle_coords, t) + end + + project_velocity_on_plane_normal!(v, system, particle, boundary_zone, + boundary_model) + end + + return system +end + +function project_velocity_on_plane_normal!(v, system, particle, boundary_zone, + boundary_model) + return v +end + +function project_velocity_on_plane_normal!(v, system, particle, boundary_zone, + boundary_model::BoundaryModelZhang) + # Project `vel` on the normal direction of the boundary zone + # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: + # "Because flow from the inlet interface occurs perpendicular to the boundary, + # only this component of interpolated velocity is kept [...]" + vel = current_velocity(v, system, particle) + vel_ = dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal + + @inbounds for dim in eachindex(vel) + v[dim, particle] = vel_[dim] + end + + return v +end + +# Interaction of boundary with other systems +function interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + particle_system::OpenBoundarySPHSystem{<:BoundaryModelZhang}, + neighbor_system, semi) + (; fluid_system, cache) = particle_system + sound_speed = system_sound_speed(fluid_system) + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_neighbor(particle_system, neighbor_system, + system_coords, neighbor_system_coords, semi; + points=each_moving_particle(particle_system)) do particle, + neighbor, + pos_diff, + distance + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds current_density(v_particle_system, particle_system, particle) + rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor) + + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + p_a = current_pressure(v_particle_system, particle_system, particle) + p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor) + + # "To avoid the lack of support near the buffer surface entirely, one may use the + # angular momentum conservative form." + dv_pressure = inter_particle_averaged_pressure(m_a, m_b, rho_a, rho_b, + p_a, p_b, grad_kernel) + + # This vanishes for particles with full kernel support + p_boundary = cache.boundary_pressure[particle] + dv_pressure_missing = 2 * p_boundary * (m_b / (rho_a * rho_b)) * grad_kernel + + # Propagate `@inbounds` to the viscosity function, which accesses particle data + dv_viscosity_ = @inbounds dv_viscosity(viscosity_model(fluid_system, + neighbor_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) + + for i in 1:ndims(particle_system) + @inbounds dv[i, + particle] += dv_pressure[i] + dv_viscosity_[i] + + dv_pressure_missing[i] + end + + # Continuity equation + vdiff = current_velocity(v_particle_system, particle_system, particle) - + current_velocity(v_neighbor_system, neighbor_system, neighbor) + + @inbounds dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel) + end + + return dv +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 14034fedc2..c27369d340 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -173,6 +173,7 @@ function create_cache_open_boundary(boundary_model, initial_condition, reference_density, reference_velocity, reference_pressure) ELTYPE = eltype(initial_condition) + NDIMS = ndims(initial_condition) prescribed_pressure = isnothing(reference_pressure) ? false : true prescribed_velocity = isnothing(reference_velocity) ? false : true @@ -184,6 +185,32 @@ function create_cache_open_boundary(boundary_model, initial_condition, prescribed_velocity=prescribed_velocity) end + if boundary_model isa BoundaryModelZhang + (; boundary_pressure) = boundary_model + + if boundary_pressure isa Function + test_result = boundary_pressure(zeros(NDIMS), 0) + if length(test_result) != 1 + throw(ArgumentError("`boundary_pressure` function in `BoundaryModelZhang` must be a scalar function")) + end + end + + boundary_pressure_function = wrap_reference_function(boundary_pressure, Val(NDIMS)) + + coords_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, + initial_condition.coordinates) + + boundary_pressure = copy(reference_value.(boundary_pressure_function, + initial_condition.pressure, + coords_svector, 0)) + + return (; prescribed_pressure=prescribed_pressure, + prescribed_density=prescribed_density, + prescribed_velocity=prescribed_velocity, + boundary_pressure_function=boundary_pressure_function, + boundary_pressure=boundary_pressure) + end + characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) @@ -229,6 +256,16 @@ end return ELTYPE end +function system_smoothing_kernel(system::OpenBoundarySPHSystem) + system.fluid_system.smoothing_kernel +end + +@inline function v_nvariables(system::OpenBoundarySPHSystem) + return v_nvariables(system, system.boundary_model) +end + +@inline v_nvariables(system, boundary_model) = ndims(system) + function reset_callback_flag!(system::OpenBoundarySPHSystem) system.update_callback_used[] = false @@ -248,6 +285,10 @@ end @inline hydrodynamic_mass(system::OpenBoundarySPHSystem, particle) = system.mass[particle] @inline function current_density(v, system::OpenBoundarySPHSystem) + return current_density(v, system.boundary_model, system) +end + +@inline function current_density(v, boundary_model, system::OpenBoundarySPHSystem) return system.density end @@ -255,6 +296,10 @@ end return system.pressure end +@inline function current_acceleration(system::OpenBoundarySPHSystem, particle) + return system.fluid_system.acceleration +end + function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t; update_from_callback=false) if !update_from_callback && !(system.update_callback_used[]) @@ -270,6 +315,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) + @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi, t) + # Update density, pressure and velocity based on the characteristic variables. # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, @@ -278,15 +325,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ v_ode, u_ode, semi, t) - - @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) - return system end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) = system -function check_domain!(system, v, u, v_ode, u_ode, semi) +function check_domain!(system, v, u, v_ode, u_ode, semi, t) (; boundary_zone, boundary_candidates, fluid_candidates) = system fluid_system = corresponding_fluid_system(system, semi) @@ -316,7 +360,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle_new = available_fluid_particles[i] convert_particle!(system, fluid_system, boundary_zone, particle, particle_new, - v, u, v_fluid, u_fluid) + v, u, v_fluid, u_fluid, t) end update_system_buffer!(system.buffer, semi) @@ -345,19 +389,22 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle_new = available_boundary_particles[i] convert_particle!(fluid_system, system, boundary_zone, particle, particle_new, - v, u, v_fluid, u_fluid) + v, u, v_fluid, u_fluid, t) end update_system_buffer!(system.buffer, semi) update_system_buffer!(fluid_system.buffer, semi) + # Since particles have been transferred, the neighborhood searches must be updated + update_nhs!(semi, u_ode) + return system end # Outflow particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, boundary_zone::BoundaryZone{OutFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) + particle_new, v, u, v_fluid, u_fluid, t) deactivate_particle!(system, particle, u) return system @@ -366,7 +413,7 @@ end # Inflow particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, boundary_zone::BoundaryZone{InFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) + particle_new, v, u, v_fluid, u_fluid, t) (; spanning_set) = boundary_zone # Activate a new particle in simulation domain @@ -383,7 +430,7 @@ end # Buffer particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, boundary_zone::BoundaryZone{BidirectionalFlow}, - particle, particle_new, v, u, v_fluid, u_fluid) + particle, particle_new, v, u, v_fluid, u_fluid, t) relative_position = current_coords(u, system, particle) - boundary_zone.zone_origin # Check if particle is in- or outside the fluid domain. @@ -402,13 +449,15 @@ end u[dim, particle] += boundary_zone.spanning_set[1][dim] end + impose_new_density!(v, u, system, system.boundary_model, particle, t) + return system end # Fluid particle is in boundary zone @inline function convert_particle!(fluid_system::FluidSystem, system, boundary_zone, particle, particle_new, - v, u, v_fluid, u_fluid) + v, u, v_fluid, u_fluid, t) # Activate particle in boundary zone transfer_particle!(system, fluid_system, particle, particle_new, v, u, v_fluid, u_fluid) @@ -445,6 +494,8 @@ function write_v0!(v0, system::OpenBoundarySPHSystem) indices = CartesianIndices(system.initial_condition.velocity) copyto!(v0, indices, system.initial_condition.velocity, indices) + write_v0!(v0, system, system.boundary_model) + return v0 end @@ -513,7 +564,8 @@ end end # When the neighbor is an open boundary system, just use the viscosity of the fluid `system` instead -@inline viscosity_model(system, neighbor_system::OpenBoundarySPHSystem) = system.viscosity +@inline viscosity_model(system, + neighbor_system::OpenBoundarySPHSystem) = neighbor_system.fluid_system.viscosity function system_data(system::OpenBoundarySPHSystem, v_ode, u_ode, semi) v = wrap_v(v_ode, system, semi)