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/semidiscretization.jl

@@ -1149,5 +1149,6 @@ function set_system_links(system::OpenBoundarySystem, semi)
                               system.buffer,
                               system.pressure_acceleration_formulation,
                               system.shifting_technique,
+                              system.pressure_model_values,
                               system.cache)
 end

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, PM}
     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
+    pressure_model    :: PM
     # 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
@@ -168,6 +169,7 @@ 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(),
+                      pressure_model=nothing,
                       rest_pressure=zero(eltype(density)),
                       reference_density=nothing, reference_pressure=nothing,
                       reference_velocity=nothing)
@@ -246,7 +248,11 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
     coordinates_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, ic.coordinates)
 
     if prescribed_pressure
-        ic.pressure .= pressure_ref.(coordinates_svector, 0)
+        if isnothing(pressure_model)
+            ic.pressure .= pressure_ref.(coordinates_svector, 0)
+        else
+            throw(ArgumentError("Setting prescribed pressure together with a pressure model is not supported."))
+        end
     end
     if prescribed_density
         ic.density .= density_ref.(coordinates_svector, 0)
@@ -258,8 +264,8 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
 
     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)
+                        pressure_model, average_inflow_velocity, prescribed_density,
+                        prescribed_pressure, prescribed_velocity)
 end
 
 function boundary_type_name(boundary_zone::BoundaryZone)
@@ -291,6 +297,10 @@ 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 !isnothing(boundary_zone.pressure_model) &&
+           boundary_zone.pressure_model.is_prescribed
+            summary_line(io, "pressure model", type2string(boundary_zone.pressure_model))
+        end
         summary_footer(io)
     end
 end
@@ -494,7 +504,11 @@ function reference_pressure(boundary_zone, v, system, particle, pos, t)
         # `pressure_reference_values[zone_id](pos, t)`, but in a type-stable way
         return apply_ith_function(pressure_reference_values, zone_id, pos, t)
     else
-        return current_pressure(v, system, particle)
+        pressure = current_pressure(v, system, particle)
+
+        # Return either the current pressure or a pressure computed by a pressure model
+        return imposed_pressure(system, system.pressure_model_values, boundary_zone,
+                                pressure, particle)
     end
 end
 

src/schemes/boundary/open_boundary/dynamical_pressure.jl

@@ -120,7 +120,9 @@ function reference_pressure(boundary_zone, v,
         # `pressure_reference_values[zone_id](pos, t)`, but in a type-stable way
         return apply_ith_function(pressure_reference_values, zone_id, pos, t)
     else
-        return rest_pressure
+        # Return either the rest pressure or a pressure computed by a pressure model
+        return imposed_pressure(system, system.pressure_model_values, boundary_zone,
+                                rest_pressure, particle)
     end
 end
 

src/schemes/boundary/open_boundary/mirroring.jl

@@ -95,14 +95,14 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelMirror
 
     @threaded semi for particle in each_integrated_particle(system)
         boundary_zone = current_boundary_zone(system, particle)
-        (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone
+        (; prescribed_density, prescribed_velocity) = boundary_zone
 
         particle_coords = current_coords(u, system, particle)
 
-        if prescribed_pressure
-            pressure[particle] = reference_pressure(boundary_zone, v, system, particle,
-                                                    particle_coords, t)
-        end
+        # The pressure can be prescribed, left unprescribed (in which case the current
+        # pressure is returned), or determined by a boundary model.
+        pressure[particle] = reference_pressure(boundary_zone, v, system, particle,
+                                                particle_coords, t)
 
         if prescribed_density
             density[particle] = reference_density(boundary_zone, v, system, particle,

src/schemes/boundary/open_boundary/open_boundary.jl

@@ -4,3 +4,4 @@ include("method_of_characteristics.jl")
 include("system.jl")
 include("dynamical_pressure.jl")
 include("rhs.jl")
+include("pressure_model.jl")

src/schemes/boundary/open_boundary/pressure_model.jl

@@ -0,0 +1,138 @@
+"""
+    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.
+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}
+    characteristic_resistance :: ELTYPE
+    peripheral_resistance     :: ELTYPE
+    compliance                :: ELTYPE
+    is_prescribed             :: Bool
+end
+
+function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance)
+    return RCRWindkesselModel(characteristic_resistance, peripheral_resistance, compliance,
+                              true)
+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::OpenBoundarySystem, v, u, semi, dt)
+    isnothing(system.pressure_model_values) && return system
+
+    calculate_flow_rate_and_pressure!(system, v, u, dt)
+
+    return system
+end
+
+function calculate_flow_rate_and_pressure!(system, v, u, dt)
+    for (zone_id, boundary_zone) in enumerate(system.boundary_zones)
+        if boundary_zone.pressure_model.is_prescribed
+            calculate_flow_rate_and_pressure!(boundary_zone.pressure_model, system,
+                                              boundary_zone, zone_id, v, u, dt)
+        end
+    end
+
+    return system
+end
+
+function calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone,
+                                           zone_id, v, u, dt)
+    dt < sqrt(eps()) && return pressure_model
+    (; particle_spacing) = system.initial_condition
+    (; characteristic_resistance, peripheral_resistance, compliance) = pressure_model
+    (; flow_rate, pressure) = system.pressure_model_values[zone_id]
+    (; face_normal, zone_origin) = boundary_zone
+
+    # Use a thin slice for the flow rate calculation
+    slice = particle_spacing * 125 / 100
+
+    # Find particles within a thin slice near the boundary face for flow rate computation
+    candidates = findall(x -> dot(x - zone_origin, -face_normal) <= slice,
+                         reinterpret(reshape, SVector{ndims(system), eltype(u)}, u))
+
+    particles_in_zone = findall(particle -> boundary_zone ==
+                                            current_boundary_zone(system, particle),
+                                each_integrated_particle(system))
+
+    intersect!(candidates, particles_in_zone)
+
+    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 imposed_pressure(system, pressure_model_values, boundary_zone,
+                          pressure, particle)
+    boundary_zone.pressure_model.is_prescribed || return pressure
+
+    zone_id = system.boundary_zone_indices[particle]
+    return pressure_model_values[zone_id].pressure[]
+end
+
+function imposed_pressure(system, pressure_model_values::Nothing, boundary_zone,
+                          pressure, particle)
+    return pressure
+end