Add RCR-Windkessel model

docs/src/refs.bib

@@ -330,6 +330,15 @@ @Article{Ganzenmueller2015
   publisher = {Elsevier {BV}},
 }
 
+@Book{Gasser2021,
+  author    = {Gasser, T. Christian},
+  title     = {Vascular Biomechanics: Concepts, Models, and Applications},
+  year      = {2021},
+  isbn      = {9783030709662},
+  doi       = {10.1007/978-3-030-70966-2},
+  publisher = {Springer International Publishing},
+}
+
 @Article{Giles1990,
   author    = {Giles, Michael B.},
   title     = {Nonreflecting boundary conditions for {Euler} equation calculations},
@@ -828,6 +837,20 @@ @Article{Wendland1995
   publisher = {Springer Science and Business Media LLC},
 }
 
+@Article{Westerhof2008,
+  author    = {Westerhof, Nico and Lankhaar, Jan-Willem and Westerhof, Berend E.},
+  journal   = {Medical \& Biological Engineering \& Computing},
+  title     = {The arterial Windkessel},
+  year      = {2008},
+  issn      = {1741-0444},
+  month     = jun,
+  number    = {2},
+  pages     = {131--141},
+  volume    = {47},
+  doi       = {10.1007/s11517-008-0359-2},
+  publisher = {Springer Science and Business Media LLC},
+}
+
 @Article{Zhang2025,
   author    = {Zhang, Shuoguo and Fan, Yu and Wu, Dong and Zhang, Chi and Hu, Xiangyu},
   journal   = {Physics of Fluids},

docs/src/systems/boundary.md

@@ -342,3 +342,9 @@ without the need to specifically identify those near the free surface.
 To further handle incomplete kernel support, for example in the viscous term of the momentum equation,
 the updated velocity of particles within the [`BoundaryZone`](@ref) is projected onto the face normal,
 so that only the component in flow direction is kept.
+
+# Pressure Models
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "boundary", "open_boundary", "pressure_model.jl")]
+```

src/TrixiParticles.jl

@@ -90,6 +90,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx
 export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring
 export HertzContactModel, LinearContactModel
 export PrescribedMotion, OscillatingMotion2D
+export RCRWindkesselModel
 export examples_dir, validation_dir
 export trixi2vtk, vtk2trixi
 export RectangularTank, RectangularShape, SphereShape, ComplexShape

src/general/gpu.jl

@@ -20,6 +20,7 @@ Adapt.@adapt_structure TotalLagrangianSPHSystem
 Adapt.@adapt_structure BoundaryZone
 Adapt.@adapt_structure SystemBuffer
 Adapt.@adapt_structure OpenBoundarySystem
+Adapt.@adapt_structure RCRWindkesselModel
 
 KernelAbstractions.get_backend(::PtrArray) = KernelAbstractions.CPU()
 function KernelAbstractions.get_backend(system::AbstractSystem)

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -203,7 +203,7 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
     end
 
     if !(reference_pressure isa Function || reference_pressure isa Real ||
-         isnothing(reference_pressure))
+         reference_pressure isa AbstractPressureModel || isnothing(reference_pressure))
         throw(ArgumentError("`reference_pressure` must be either a function mapping " *
                             "each particle's coordinates and time to its pressure, " *
                             "or a scalar"))
@@ -213,10 +213,15 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
             if length(test_result) != 1
                 throw(ArgumentError("`reference_pressure` function must be a scalar function"))
             end
+            pressure_ref = reference_pressure
+        elseif reference_pressure isa AbstractPressureModel
+            pressure_ref = reference_pressure
+            pressure_ref.pressure[] = rest_pressure
+        else
+            # We need this dummy for type stability reasons
+            pressure_dummy = convert(ELTYPE, Inf)
+            pressure_ref = wrap_reference_function(reference_pressure, pressure_dummy)
         end
-        # We need this dummy for type stability reasons
-        pressure_dummy = convert(ELTYPE, Inf)
-        pressure_ref = wrap_reference_function(reference_pressure, pressure_dummy)
     end
 
     if !(reference_density isa Function || reference_density isa Real ||

src/schemes/boundary/open_boundary/open_boundary.jl

@@ -1,3 +1,4 @@
+include("pressure_model.jl")
 include("boundary_zones.jl")
 include("mirroring.jl")
 include("method_of_characteristics.jl")

src/schemes/boundary/open_boundary/pressure_model.jl

@@ -0,0 +1,139 @@
+abstract type AbstractPressureModel end
+
+"""
+    RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance)
+
+The `RCRWindkessel` model is a biomechanical lumped-parameter representation
+that captures the relationship between pressure and flow in pulsatile systems, e.g., vascular systems
+and is used to compute the pressure in a [`BoundaryZone`](@ref).
+It is derived from an electrical circuit analogy and consists of three elements:
+
+- characteristic resistance (``R_1``): Represents the proximal resistance at the vessel entrance.
+  It models the immediate pressure drop that arises at the entrance of a vessel segment,
+  due either to a geometric narrowing or to a mismatch in characteristic impedance between adjacent segments.
+  A larger ``R_1`` produces a sharper initial pressure rise at the onset of flow.
+- peripheral resistance (``R_2``): Represents the distal resistance,
+  which controls the sustained outflow into the peripheral circulation and thereby determines the level of the mean pressure.
+  A high ``R_2`` maintains a higher pressure (reduced outflow), whereas a low ``R_2`` allows a faster pressure decay.
+- compliance (``C``): Connected in parallel with ``R_2`` and represents the capacity of elastic walls
+  to store and release volume; in other words, it models the "stretchiness" of 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.
+
+ 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).
+
+# Keywords
+- `characteristic_resistance`: characteristic resistance (``R_1``)
+- `peripheral_resistance`: peripheral resistance (``R_2``)
+- `compliance`: compliance (``C``)
+"""
+struct RCRWindkesselModel{ELTYPE <: Real, P, FR} <: AbstractPressureModel
+    characteristic_resistance :: ELTYPE
+    peripheral_resistance     :: ELTYPE
+    compliance                :: 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)
+    pressure = Ref(zero(compliance))
+    flow_rate = Ref(zero(compliance))
+    return RCRWindkesselModel(characteristic_resistance, peripheral_resistance, compliance,
+                              pressure, flow_rate)
+end
+
+function Base.show(io::IO, ::MIME"text/plain", pressure_model::RCRWindkesselModel)
+    @nospecialize pressure_model # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, pressure_model)
+    else
+        summary_header(io, "RCRWindkesselModel")
+        summary_line(io, "characteristic_resistance",
+                     pressure_model.characteristic_resistance)
+        summary_line(io, "peripheral_resistance",
+                     pressure_model.peripheral_resistance)
+        summary_line(io, "compliance", pressure_model.compliance)
+        summary_footer(io)
+    end
+end
+
+function update_pressure_model!(system, v, u, semi, dt)
+    # Skip update for the first time step
+    dt < sqrt(eps()) && return system
+
+    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)
+        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)
+    end
+
+    return system
+end
+
+function calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, v, u, 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
+    (; 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[]
+    flow_rate[] = current_flow_rate
+
+    # Calculate new pressure according to eq. 22 in Zhang et al. (2025)
+    R1 = characteristic_resistance
+    R2 = peripheral_resistance
+    C = compliance
+
+    term_1 = (1 + R1 / R2) * flow_rate[]
+    term_2 = C * R1 * (flow_rate[] - previous_flow_rate) / dt
+    term_3 = C * previous_pressure / dt
+    divisor = C / dt + 1 / R2
+
+    pressure_new = (term_1 + term_2 + term_3) / divisor
+
+    pressure[] = pressure_new
+
+    return system
+end
+
+function (pressure_model::RCRWindkesselModel)(x, t)
+    return pressure_model.pressure[]
+end

src/schemes/boundary/open_boundary/system.jl

@@ -72,9 +72,8 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                                                BoundaryModelDynamicalPressureZhang ?
                                                shifting_technique(fluid_system) : nothing)
     boundary_zones_ = filter(bz -> !isnothing(bz), boundary_zones)
-    reference_values_ = map(bz -> bz.reference_values, boundary_zones_)
 
-    initial_conditions = union((bz.initial_condition for bz in boundary_zones)...)
+    initial_conditions = union((bz.initial_condition for bz in boundary_zones_)...)
 
     buffer = SystemBuffer(nparticles(initial_conditions), buffer_size)
 
@@ -86,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,
-                                        reference_values_)...)
+                                        boundary_zones_)...)
 
     fluid_system_index = Ref(0)
 
@@ -101,6 +100,8 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
     # Create new `BoundaryZone`s with `reference_values` set to `nothing` for type stability.
     # `reference_values` are only used as API feature to temporarily store the reference values
     # in the `BoundaryZone`, but they are not used in the actual simulation.
+    # The reference values are extracted above in the "create cache" function
+    # and then stored in `system.cache` as a `Tuple`.
     boundary_zones_new = map(zone -> BoundaryZone(zone.initial_condition,
                                                   zone.spanning_set,
                                                   zone.zone_origin,
@@ -113,7 +114,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                                                   zone.prescribed_density,
                                                   zone.prescribed_pressure,
                                                   zone.prescribed_velocity),
-                             boundary_zones)
+                             boundary_zones_)
 
     return OpenBoundarySystem(boundary_model, initial_conditions, fluid_system,
                               fluid_system_index, smoothing_kernel, smoothing_length, mass,
@@ -130,8 +131,9 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
-function create_cache_open_boundary(boundary_model, fluid_system,
-                                    initial_condition, reference_values)
+function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
+                                    boundary_zones)
+    reference_values = map(bz -> bz.reference_values, boundary_zones)
     ELTYPE = eltype(initial_condition)
 
     # Separate `reference_values` into pressure, density and velocity reference values
@@ -217,6 +219,10 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySystem)
             summary_line(io, "density diffusion", density_diffusion(system))
             summary_line(io, "shifting technique", shifting_technique(system))
         end
+        for (i, pm) in enumerate(system.cache.pressure_reference_values)
+            !isa(pm, AbstractPressureModel) && continue
+            summary_line(io, "pressure model", type2string(pm) * " (in boundary zone $i)")
+        end
         summary_footer(io)
     end
 end
@@ -300,6 +306,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode
 
     @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
 
+    update_pressure_model!(system, v, u, semi, integrator.dt)
+
     # Update density, pressure and velocity based on the specific boundary model
     @trixi_timeit timer() "update boundary quantities" begin
         update_boundary_quantities!(system, boundary_model, v, u, v_ode, u_ode, semi, t)

src/util.jl

@@ -11,6 +11,21 @@ end
 
 @inline foreach_noalloc(func, collection::Tuple{}) = nothing
 
+# Same as `foreach(enumerate(something))`, but without allocations.
+# Note that compile times may increase if this is used with big tuples.
+@inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1)
+@inline foreach_enumerate(func, collection::Tuple{}, index) = nothing
+
+@inline function foreach_enumerate(func, collection, index)
+    element = first(collection)
+    remaining_collection = Base.tail(collection)
+
+    func((index, element))
+
+    # Process remaining collection
+    foreach_enumerate(func, remaining_collection, index + 1)
+end
+
 # Returns `functions[index](args...)`, but in a type-stable way for a heterogeneous tuple `functions`
 @inline function apply_ith_function(functions, index, args...)
     if index == 1

test/schemes/boundary/open_boundary/open_boundary.jl

@@ -1,3 +1,4 @@
 include("characteristic_variables.jl")
 include("mirroring.jl")
 include("boundary_zone.jl")
+include("pressure_model.jl")