Dynamical pressure (open) BC

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

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)

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.
 

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`)

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

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