# A high-order space-time Fourier continuation approach for one-dimensional hemodynamics and wave propagation in the entire human circulatory system

Faisal Amlani, Niema M Pahlevan

► To cite this version:

Faisal Amlani, Niema M Pahlevan. A high-order space-time Fourier continuation approach for one-dimensional hemodynamics and wave propagation in the entire human circulatory system. Physics of Fluids, 2025, 37 (7), pp.071921. 10.1063/5.0273041. hal-05184930

HAL Id: hal-05184930

https://hal.science/hal-05184930v1

Submitted on 24 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A high-order space-time Fourier continuation approach for one-dimensional hemodynamics and wave propagation in the entire human circulatory system

Faisal Amlani ; Niema M. Pahlevan

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for updates

Physics of Fluids 37, 071921 (2025)

https://doi.org/10.1063/5.0273041

<figure><img src="image_5.png" /></figure>

CHORUS

<figure><img src="image_6.png" /></figure>

View
Online

<figure><img src="image_7.png" /></figure>

Export
Citation

## Articles You May Be Interested In

* Effects of boundary conditions on the transmitral pressure gradient for numerical simulation of flow in a left heart model

* Physics of Fluids (November 2024)

* Stokesian peristaltic pumping in a three-dimensional tube with a phase-shifted asymmetry

* Physics of Fluids (August 2011)

* Kinematic responses of an autorotating samara to concentrated crosswind

* Physics of Fluids (November 2022)

<figure><img src="image_8.png" /></figure>

AIP
Publishing

Physics of Fluids

Special Topics Open
for Submissions

Learn More

<figure><img src="image_9.png" /></figure>

AIP
Publishing
# A high-order space-time Fourier continuation approach for one-dimensional hemodynamics and wave propagation in the entire human circulatory system

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

Cite as: Phys. Fluids **37**, 071921 (2025); doi: 10.1063/5.0273041

Submitted: 28 March 2025 · Accepted: 4 June 2025 ·

Published Online: 24 July 2025

<figure><img src="image_12.png" /></figure>

View Online

<figure><img src="image_13.png" /></figure>

Export Citation

<figure><img src="image_14.png" /></figure>

CiteMark

Faisal Amlani<fnref n="1" /> ID and Niema M. Pahlevan<fnref n="2" />,<fnref n="3" />,a ID

## AFFILIATIONS

## ABSTRACT

Accurate and robust numerical simulation of hemodynamics is of interest for the identification and investigation of important physical or physiological quantities and their relationships to various cardiac and vascular functions. Comprehensive analysis requires efficiently simulating pressure and flow waves that travel and reflect through a complex network of vessels with varying geometric/material properties. This work introduces a new numerical methodology for modeling such wave propagation throughout the (closed-loop) circulation, employing a fast high-order (pseudo)spectral approach for resolving the well-established reduced-order one-dimensional Navier-Stokes hemodynamics formulations (coupled to hyperelastic tube laws) that govern the corresponding fluid-structure dynamics in each vascular segment. The model includes both systemic and pulmonary circulations and four heart chambers and four heart valves. Together with a correspondingly high-order treatment of multiscale zero-dimensional boundary conditions based on time-dependent ordinary differential equations, the overall solver has a number of attractive qualities: high-order accuracy in time and space; fast Fourier transform (FFT)-level computational efficiency; little to no numerical pollution errors (faithfully preserving the diffusion and dispersion characteristics of the underlying continuous operators); relatively mild Courant-Friedrichs-Lewy constraints for explicit temporal integration methods; robustness to extreme physiological parameters; and stable incorporation of the nonlinear and nonstationary coupling to other cardiovascular system components (e.g., heart chambers, valves, and microvasculature). The convergence properties, computational performance, and physiological accuracy of the proposed framework are demonstrated through a variety of numerical experiments that include applications to community benchmark problems previously proposed for mutual validation with other solvers (and three-dimensional or in vitro reference solutions).

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0273041

## I. INTRODUCTION

Comprehensive modeling of blood flow is essential for capturing arterial wave dynamics and understanding hemodynamic behavior across the cardiovascular system. It plays a critical role in uncovering key physiological mechanisms and detecting pathophysiological changes that drive cardiovascular function and dysfunction. Maximizing the physiological relevancy of computational models entails simulating pressure and flow waves within a complete circulatory system, where such waves reflect and attenuate through a complex

vascular network as they encounter various bifurcated or trifurcated junctions and changes in geometry or material properties. In order to make this endeavor tractable, several well-established reduced-order mathematical models<fnref n="1" />⁻⁸ have been developed. These models use cross section-averaged, axisymmetric Navier-Stokes partial differential equation (PDE) formulations coupled to ordinary differential equation (ODE)-based boundary conditions to represent the effects of end organs, valves, and microvasculature. Such approaches are noted for their computational efficiency, which is achieved by decomposing the

<footnote n="1">¹Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS – Laboratoire de Mécanique Paris-Saclay, 91190 Gif-sur-Yvette, France</footnote>

<footnote n="2">²Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, California 90089, USA</footnote>

<footnote n="3">³Division of Cardiovascular Medicine, Department of Medicine of the Keck School of Medicine, University of Southern California, Los Angeles, California 90017-2395, USA</footnote>

<footnote>a) Author to whom correspondence should be addressed: pahlevan@usc.edu</footnote>
domain into segments of tubes of different lengths whose properties are described by single axial coordinates.¹⁻⁵

A number of well-established and successful numerical methods for resolving the corresponding PDE systems have been proposed, including those based on the method of characteristics,⁹ finite differences¹,²,¹⁰⁻¹² (FD), finite volumes (FV),¹³,¹⁴ and finite elements (FE).¹⁵⁻²² Many commonly used solvers—such as those based on FE/FV/FD—are known to suffer from significant numerical pollution²³ (e.g., numerical dispersion²⁴ or diffusion²⁵): phase errors accumulate with each cycle of a propagating wavetrain, necessitating an ever-increasing dense spatial resolution (i.e., more points per wavelength) to maintain accuracy over longer distances and time periods.²⁴,²⁶ As a result, the computational cost can rise substantially, undermining the efficiency typically expected from these formulations. For further details on both the theoretical (analytical) and numerical aspects of reduced-order hemodynamic wave propagation, the reader is referred to some excellent overviews found elsewhere.⁸,¹⁴,²⁷⁻²⁹

Indeed, even for one-dimensional (1D) formulations, common solution strategies¹¹,¹⁴,²⁹ are often low-order in either space or time (or both) in the interest of computational speed. For similar reasons, those that are ostensibly high-order (or can easily achieve higher orders) sacrifice such properties (i.e., are tuned down) as configurations become more complex—particularly as the number of modeled vascular segments becomes large.¹⁴ Hence, this work introduces a new pseudospectral numerical approach that is high-order in both space and time (including in its treatment of ODE-based boundary conditions) for complex, multi-segment open-loop and closed-loop 1D hemodynamics models. Such a methodology is based on the concept of Fourier continuation (FC)²⁶,³⁰,³¹ for spatial resolution, which can produce highly-accurate trigonometric interpolations (Fourier series representations) of non-periodic functions while avoiding the notorious Gibb's "ringing effect." FC-based techniques expand the applicability of Fourier-based PDE solvers toward general (physical) boundary conditions²⁴⁻²⁶ and computational domains.²⁶,³¹ Corresponding solutions provide fast Fourier transform (FFT)-speed high-order accuracy and faithfully capture the dispersion or diffusion characteristics of the underlying continuous problems (i.e., effectively no pollution errors). A number of both implicit- and explicit-in-time FC-based solvers have been constructed for a variety of physical equations including those for compressible and incompressible fluids,³¹,³² Euler equations³³,³⁴ Cauchy–Navier elastic wave equations,²⁶ and advection–diffusion equations²⁵ (among others). For hemodynamics problems, the FC method has been applied only to a single segment with a single chamber heart model and a linearized set of fluid–structure equations,²⁴ although to good effect. The complete circulation solver proposed herein represents a significant and expansive extension of an FC methodology to highly nonlinear modeling of the entire (closed-loop) system (which can also treat multi-segment open-loop models).

Pulse waves in the arterial network (which can be measured *in vivo* using both invasive and noninvasive devices²⁷) are heavily affected by the heart and vascular system. Reduced-order 1D modeling approaches are well-validated vs 3D or *in vivo* data⁶,¹⁰,¹⁷,³⁵ and can be formulated to be patient-specific.³⁶ They have been successfully employed (mostly as open-loop models) in a variety of applications,²⁷ including for studying the impact of aging on hemodynamics (where arterial stiffness can change³⁷), for providing mechanistic insight into physiology or pathophysiology (such as heart–aorta–brain coupling³⁸),

for detecting conditions such as aneurysms or stenosis,³⁹ and for investigating the efficacy of interventions or medical devices (such as grafts⁴⁰). For an excellent overview of their utility, the reader is referred to a comprehensive VascAgeNet review paper.²⁷ The FC-based methodology introduced in this work further extends the applicability of such 1D models by providing resolution for a wide range of physical parameters without the need for any numerical tuning and by providing, for the first time, numerically dispersionless resolution of the complete circulatory system. The latter quality enables the accurate simulation of long-time pulse wave propagation that is suitable for incorporating feedback regulatory models at different timescales (i.e., control loops) and for iteratively solving inverse problems related to arterial wave reflections.⁴¹,⁴² Additionally, the fast, robust nature of the approach may also enable coupling to multiphysical configurations (e.g., mass transport models) as well as also enable treatment of broader, multiscale contexts (such as facilitating the use of an entire closed-loop circulatory system model as a physiological boundary condition for a 3D model⁴³).

The proposed solver and its validation are inspired by two (recent) seminal works in the 1D hemodynamics modeling community: the closed-loop (entire circulation) formulations proposed by Mynard and collaborators⁸ (also invoked by others⁷,³⁵,⁴⁴,⁴⁵) and the open-loop (arterial) benchmark formulations of Alastruey and collaborators¹⁴ [including their powerful discontinuous Galerkin (DCG) code¹⁶ based on the open-source Nektar++⁴⁷,⁴⁸]. Indeed, in the spirit of further contributing to these endeavors, this manuscript details an alternative numerical (and algorithmic) approach that ultimately provides the highest space-time accuracy proposed thus far for both open-loop¹⁴ and closed-loop formulations,⁸ with consistent and stable coupling to microvascular, coronary, hepatic, and four-chamber heart models. The solver enables rapid, faithful resolution of complex 1D configurations, with the ultimate goal of enabling, through efficient and high-order computational tools, fast parametric exploration of physiological parameters as well as robust coupling to higher-dimensional (3D) or multiphysical configurations.⁴³ Section II presents the general forms of the governing fluid–structure PDE systems as well as boundary conditions that include junctions between vessel segments, various vascular bed ODE models, and coupled four-chamber heart-valve ODE models with septal interactions⁷ (where the lattermost is presented in a newly proposed ODE formulation that avoids the need for low-order numerical integration to determine the volume as in other contributions⁷,⁸,⁴⁴,⁴⁵,⁴⁹). Section III describes the overall solver based on FC and on high-order temporal integration for all components of the modeling framework. A number of performance studies attesting to convergence, computational cost, and physiological accuracy of the new solver are presented in Sec. IV, which includes dispersion experiments and a suite of community benchmark problems¹⁴ with comparisons to other methods as well as reference data. Section V presents examples of simulated closed-loop configurations illustrating robustness of the solver with respect to varying physiological parameters of the heart, before concluding with Sec. VI.

## II. GOVERNING EQUATIONS

The complex, multiscale modeling framework of the entire closed-loop circulation, replete with a four-chamber heart model, is directly inspired by well-validated, physiologically-accurate formulations.⁷,⁸,⁴⁴ Other researchers have adopted similar systems (minus a four-chamber heart model) for open-loop studies (modeling only the
forward circulation and not the venous return)—benchmark cases derived from these models are considered in Sec. IV C, all of which are part of a comprehensive community study paper.¹⁴ This work also uniquely formulates heart components in differential form to enable high-order integration (avoiding the use of a two-point trapezoidal integration for computing chamber volume⁷,⁸,⁴⁴,⁴⁵,⁴⁹). The subsequent high-order numerical treatment of all solver elements is detailed in Sec. III.

## A. 1D (2D-axisymmetric) fluid-structure dynamics in a vessel segment

Propagation of flow and pressure waves in the circulatory system is commonly modeled by decomposing the domain into segments of tubes of different lengths whose properties are described by single axial coordinates.¹⁻⁵,⁸,¹⁴ Along a segment of length $L_{seg}$, for a mean pressure $P = P(x, t)$ and a mean velocity $U = U(x, t)$ over a cross-sectional area $A = A(x, t)$ (yielding a flow rate of $Q(x, t) = AU$), such propagation can be expressed as a reduced-order nonlinear Navier-Stokes system in 1D given by⁸,¹⁴

$$ \left\{ \begin{aligned} -\frac{\partial A}{\partial t} &= \frac{\partial A}{\partial x} U + \frac{\partial U}{\partial x} A, \\ -\frac{\partial U}{\partial t} &= U \frac{\partial U}{\partial x} + \frac{1}{\rho} \frac{\partial P}{\partial x} + \frac{2(\zeta + 2)\mu\pi U}{\rho A}, \end{aligned} \right. \quad (1) $$

$\forall (x, t) \in [0, L_{seg}] \times \mathbb{R}^+$, where $\rho \in \mathbb{R}$ is a (constant) blood density, $\mu \in \mathbb{R}$ is a (constant) blood viscosity, and $\zeta \in \mathbb{R}$ is a given (constant) order of an assumed axisymmetric velocity profile.⁶,¹⁴,⁵⁰ The lattermost is assumed constant (as is standard in 1D hemodynamics models⁶,⁸,¹⁴,⁵⁰) since the vessel diameters considered in this work are no smaller than 1.352 mm in the open-loop configurations¹⁴ (Sec. IV C) and 0.500 mm in the closed-loop configurations⁸ (Sec. V). The source materials⁸,¹⁴ that propose these configurations also assume constant values, as the cell-free layer is primarily relevant to microvasculature⁵¹ (modeled here using the 0D lumped-parameter Windkessels described in Sec. II B 3).

The internal pressure $P = P(x, t)$, which accounts for the fluid-structure interaction of the problem, can be given by a number of constitutive tube (power) laws describing linear or nonlinear elastic wall deformation.⁸,¹⁴,²⁰,²¹ Following others,⁷,⁴⁹,⁵⁰,⁵²,⁵³ the solver proposed in this work closes the system by a generalized nonlinear (hyper)elastic power law relation given by

$$ P(x, t) = P_{\text{ext}}(x) + P_0(x) + \frac{2\rho c_0^2(x)}{b(x)} \left[ \left( \frac{A(x, t)}{A_0(x)} \right)^{b(x)/2} - 1 \right], \quad (2) $$

where $P_{\text{ext}}(x, t)$ is an external pressure constraint (if it exists), $P_0(x)$ is a reference (often diastolic) pressure, $A_0(x)$ is the reference area, and $c_0(x)$ is a reference wave speed in the vessel (which accounts for wall material properties of stiffness and thickness and can be written in such terms,¹⁴ see Sec. IV C). The corresponding instantaneous wave speed or local pulse wave velocity (PWV) is formulated based on these reference values and is given by

$$ c(x, t) = c_0(x) \left( \frac{A(x, t)}{A(x)} \right)^{b(x)/4}, \quad (3) $$

which is also a well-known biomarker (non-invasively measured⁵⁴) quantifying arterial stiffness and correlating with cardiovascular disease risk/progression.²⁷,³⁸,⁵⁴ The exponent $b(x)$ in Eqs. (2) and (3) is often taken to be equal to unity¹⁴⁻¹⁶,²⁰,⁵⁵ ($b(x) = 1 \forall x \in [0, L_{seg})$, corresponding to a thin, homogeneous wall with constant Young's modulus. Such a value is only employed in this work for the benchmark cases in Sec. IV C, which have been proposed¹⁴ with its usage since the corresponding (open-loop) configurations only concern the systemic arteries, and $b = 1$ adequately captures the associated pressure range⁴⁹ (which may not be true for the systemic venous return nor the pulmonary systems⁷,⁸). However, for the complete (closed-loop) circulation, following others,⁷,⁸,⁴⁴,⁴⁵ the solver proposed here adopts a more generalized exponent to account for collapse in vessels as well as wall distension over a wider range of pressures (i.e., more physiologically suitable for the additional modeled circulatory components⁷,⁸):

$$ b(x) = \frac{2\rho(c_0(x))^2}{P_0(x) - P_{\text{collapse}}}. \quad (4) $$

Here, $P_{\text{collapse}}$ is an estimated pressure for which $A = 0$ and is assumed to be $P_{\text{collapse}} = -10$ mm Hg throughout this work (following others⁷,⁸,⁴⁵), unless indicated otherwise.

## B. Boundary conditions

Boundary condition types include vessel junctions (branching), ODE-based heart and valve models, as well as ODE-based vascular bed models (including the hepatic/liver and coronary systems) that capture end-organ and microcirculation effects. The lattermost is used in full circulation models for connecting arterial to venous systems and in open circulation models to model terminal vasculature at the end of the systemic arterial system. Governing formulations for all cases are summarized in the following.

**Remark II.1.** Certain boundary models described in this section invoke the use of the term “0D,” which is common parlance in the hemodynamics community to refer to models that do not take into account any spatial variation, i.e., they are merely lumped-parameter time-dependent ODEs.

### 1. Riemann invariants and hyperbolic characteristics

In order to facilitate the application of the boundary conditions for vessel branching as well as the ODEs of the heart, valve, and vascular beds (Sec. II B), Riemann invariants are employed to enforce compatibility of propagating characteristics, derived from the hyperbolic nature of the governing PDE system.²,¹⁶,⁴⁴,⁴⁹,⁵⁵⁻⁵⁷ Such characteristic variables (backward-going, $-$, and forward-going, $+$) for Eq. (1) are given by⁴⁴,⁴⁹

$$ w_{\pm}(x, t) = u \pm \frac{4}{b(x)} (c(x, t) - c_0(x)), \quad (5) $$

where $c_0(x)$ is given for a segment and $c(x, t)$ and $b(x)$ are given by Eqs. (3) and (4), respectively. The corresponding solutions for velocity $U(x, t)$ and area $A(x, t)$ can be expressed in terms of such characteristics by⁴⁴,⁴⁹
$$ 
\begin{aligned}
U(x, t) &= \frac{w_{+}(x, t) + w_{-}(x, t)}{2}, \\
A(x, t) &= A_{0}(x) \left[ \frac{b(x)}{8c_{0}(x, t)} (w_{+}(x, t) - w_{-}(x, t)) + 1 \right]^{4/b(x)},
\end{aligned}
\quad (6) $$

where $A_0(x)$ is the initial cross-sectional area of the segment of interest.

## 2. Vessel branching

A multi-segment vascular network is of a fractal (tree) structure that includes branching points representing bifurcations or trifurcations in both directions (splitting and merging). Such junctions effectively act as mathematical discontinuities in material properties or cross-sectional area, where one must physically enforce conservation of mass as well as continuity of total pressure. In its most general form, given $N_p$ parent vessels $p_i$, $i = 1, ..., N_p$, which split into $N_d$ daughter vessels $d_j$, $j = 1, ..., N_d$, this leads to the expressions given by$^{14,44,45,49}$

$$ \sum_{i=1}^{N_p} \left( P_{p_i} + \frac{\rho}{2} U_{p_i}^2 \right) = \sum_{j=1}^{N_d} \left( P_{d_j} + \frac{\rho}{2} U_{d_j}^2 \right) \quad (7) $$

and

$$ \sum_{i=1}^{N_p} A_{p_i} U_{p_i} + \sum_{j=1}^{N_d} A_{d_j} U_{d_j} = 0, \quad (8) $$

which are evaluated at the respective boundaries of each vessel segment (end of a parent, beginning of a daughter). For closing (i.e., fully determining) the system, these equations are supplemented via the solution of the corresponding Riemann invariant problem that enforces compatibility of propagating characteristics$^{2,16,44,49,55-57}$ (Sec. II B 1) to provide additional expressions given by

$$ 
\begin{aligned}
w_{+,p_i} &= U_{p_i} + \frac{4}{b_{p_i}} \left( c \left( A_{p_i}^{1/4} \right) - c_{0,p_i} \right) \quad \text{and} \\
w_{-,d_j} &= U_{d_j} - \frac{4}{b_{d_j}} \left( c \left( A_{d_j}^{1/4} \right) - c_{0,d_j} \right),
\end{aligned}
\quad (9) $$

where $w_{+,p_i}$ represents the outgoing characteristic from the right boundary of a parent vessel $p_i$ and $w_{-,d_j}$ represents the outgoing characteristic from the left boundary of a daughter vessel $d_j$ [A and U can both be represented at the boundary in terms of both incoming and outgoing characteristics, see Eq. (6)].

## 3. OD Windkessel-based vascular bed modeling

Vascular beds representing organs and microcirculation employ well-established OD lumped-parameter models based on Windkessel formulations. For the complete circulation, such beds connect the systemic arterial system to the systemic venous system and the pulmonary arterial system to the pulmonary venous system. The solver proposed in this work directly employs a well-established formulation$^8$ that has been used extensively.$^{7,35,44,45}$ For a vascular bed connecting $N_{art}$ supplying arteries to $N_{ven}$ draining veins, the corresponding flows are expressed as ODEs governing their respective vascular bed pressures $P_{art}(t)$ and $P_{ven}(t)$ that are given by$^{7,8}$

$$ \sum_{j=1}^{N_{art}} Q_{j,\text{out}}(t) = Q_{\text{cap}}(t) + C_{\text{art}} \frac{dP_{\text{art}}}{dt}(t) \quad (10) $$

and

$$ \sum_{k=1}^{N_{ven}} Q_{k,\text{in}}(t) = Q_{\text{cap}}(t) - C_{\text{ven}} \frac{dP_{\text{ven}}}{dt}(t), \quad (11) $$

where $Q_{cap}(t)$ is the capillary flow (calculated via Ohm's law)$^7$ across the vascular bed resistance $R_{cap}(t)$ and the constants $C_{art}$, $C_{ven}$ are lumped arterial and venous compliances, respectively. As in other work,$^{7,35,44,45}$ $R_{cap}(t)$ depends on the transmural pressure $P_{tm}(t)$ and a zero-flow constant pressure $P_{zf}$:

$$ R_{\text{cap}}(t) = \begin{cases} R_0 \left( \frac{P_{\text{tm}}(0) - P_{\text{zf}}}{P_{\text{tm}}(t) - P_{\text{zf}}} \right), & p_{\text{tm}}(t) > P_{\text{zf}}, \\ \infty, & p_{\text{tm}}(t) \le P_{\text{zf}}. \end{cases} $$

From the solutions $P_{art}(t)$ and $P_{ven}(t)$ derived from their respective ODEs, the pressures and flows at each connecting outlet of artery j and inlet of vein k are governed by Ohm's law and are given by

$$ Q_{j,\text{out}}(t) = A_{j,\text{out}}(t) U_{j,\text{out}}(t) = \frac{P_{j,\text{out}}(t) - P_{\text{art}}(t)}{Z_{j,\text{out}}(t)} \quad (12) $$

and

$$ Q_{k,\text{in}}(t) = A_{k,\text{in}}(t) U_{k,\text{in}}(t) = \frac{P_{\text{ven}}(t) - P_{k,\text{in}}(t)}{Z_{k,\text{in}}(t)}, \quad (13) $$

where $Z_{j,\text{out}}(t)$, $Z_{k,\text{in}}(t)$ are the characteristic impedances of the outgo- and incoming vessels ($Z = \rho c/A$ at the outlet or inlet connecting point).

**Remark II.2.** For the open-loop benchmark studies of Sec. IV C, classical one-sided three-element Windkessels$^{14}$ terminating in a constant pressure, with constant resistances and compliances, are considered for boundary conditions at terminal segments of systemic arterial configurations (i.e., no connections to the venous systems). Each Windkessel model only has one supplying segment, whose outflow $Q(t)$ is related to vascular bed pressure by a similar ODE given by

$$ Q_{\text{out}}(t) = Q_{\text{cap}}(t) + C_{\text{art}} \frac{dP_{\text{art}}}{dt}(t), $$

where

$$ Q_{\text{cap}} = \frac{P_{\text{art}}(t) - P_0}{R_2} $$

for a constant bed resistance $R_2 \in \mathbb{R}^+$ and constant (given) terminal pressure $P_0$ at which flow to the microcirculation ceases.$^{14}$ The corresponding pressure and flow at the outlet of the terminal vessel is given by

$$ Q_{\text{out}}(t) = A_{\text{out}}(t) U_{\text{out}}(t) = \frac{P_{\text{out}}(t) - P_{\text{art}}(t)}{R_1}, $$

where $R_1 \in \mathbb{R}^+$ is a constant (given) characteristic impedance at the end point of the segment, i.e., in the notation above, $R_1 = Z_{\text{out}}$ for a constant $Z_{\text{out}}$. More details can be found elsewhere$^{14,49}$ and the
references therein. The numerical resolution of such a formulation is similarly treated to high-order through the methodology outlined in Sec. III B 2.

## 4. OD Windkessel-based hepatic bed modeling

Although similar in all other respects to the generic bed model, the hepatic (liver) bed model in the full circulation solver contains both arterial and venous inlets, whose pressure difference can be large.⁴⁴,⁴⁹ In order to achieve a common portal pressure in the lobules, such a configuration is modeled identically to the Windkessel formulation above, but with an additional resistance and compliance for the hepatic artery.⁸ The sum of the flows from the portal vein and the hepatic artery supplies the bed. For further technical details on the specific formulation of the hepatic bed model employed here, the reader is referred to the original reference works.⁸,⁴⁹

## 5. OD Windkessel-based coronary bed modeling

Similarly, formulations of intramyocardial coronary beds⁸,³⁵ are based on similar Windkessel elements directly adopted from other solvers,⁸,³⁵,⁴⁴,⁴⁵ consisting of three transmural layers (subendocardium, midwall, subepicardium). Each layer $ℓ$ contains two compartments, each with two compliances $C_{ℓ,1}$, $C_{ℓ,2}$ and two time-varying resistances $R_{ℓ,1}(t)$, $R_{ℓ,2}(t)$ that share a middle resistance $R_{ℓ,m}(t)$. The resistances are non-linearly dependent on volume and are given by

$$R_{\ell,i}(t) = R_{\ell,0,i} \frac{V_{\ell,0,i}^2}{V_{\ell,i}(t)^2}, \quad i = 1, 2 \quad (14)$$

and

$$R_{\ell,m}(t) = R_{\ell,0,m} \left( \gamma \frac{V_{\ell,0,1}^2}{V_{\ell,1}(t)^2} + (\gamma - 1) \frac{V_{\ell,0,2}^2}{V_{\ell,2}(t)^2} \right), \quad \gamma \in (0, 1) \quad (15)$$

where

$$V_{\ell,i}(t) = V_{\ell,0,i} + \int_{0}^{t} C_{\ell,i} \frac{dP_{\ell,tm,i}}{dt}(\tau) d\tau, \quad i = 1, 2,$$

and $R_{\ell,0,i}$, $R_{\ell,0,m}$, and $V_{\ell,0,i}$ are reference values. Here, $P_{\ell,tm,i}(t)$ is given by $P_{\ell,tm,i}(t) = P(t) - P_{\ell,im}(t)$, where the applied time-varying intramyocardial pressure $P_{\ell,im}(t)$ accounts for the change in volume by external compression of the heart and is the sum of a cavity-induced extracellular pressure and an intracellular pressure that is proportional to effective chamber elastance (in terms of ventricular pressure and volume). Further details can be found elsewhere.⁸,³⁵

## 6. OD valve modeling

The models for the aortic, pulmonary, mitral, and tricuspid valves employed by the solver are based on Bernoulli's equation and are those which have been proposed in the context of a simple closed-loop formulation by Mynard et al.⁷ (also used in other contributions⁸,⁴⁴,⁴⁵). The relationship between flow $Q_{valve}$ and pressure difference $\Delta P_{valve}(t) = P_{in}(t) - P_{out}(t)$ across a valve v is given by⁷

$$\Delta P_v(t) = B_v(t) Q_v(t) |Q_v(t)| + L_v \frac{dQ_v}{dt}(t) \quad (16)$$

for a Bernoulli resistance $B_v = \rho/2A_v^2$ and blood inertance $L_v = \rho\ell_v/A_v$, where $\ell_v$ is the effective valve length and $A_v(t)$ is a time-varying effective orifice area. The lattermost models the valve dynamics via a state variable $\xi_v(t)$ that is expressed in terms of a given minimum and maximum valve area by

$$A_v(t) = (A_{v,max} - A_{v,min}) \xi_v(t) + A_{v,min}.$$

Here, $\xi_v(t) \in [0, 1]$ represents the state of the valve, which is equal to unity when it is fully open and zero when it is fully closed. The corresponding valve dynamics during opening and closing are modeled by⁷

$$\begin{cases} \frac{d\xi_v}{dt}(t) = K_{v,open}(1 - \xi_v(t))\Delta P_v(t), & \Delta P_v(t) \ge 0, \\ \frac{d\xi_v}{dt}(t) = K_{v,close}\xi_v(t)\Delta P_v(t), & \Delta P_v(t) < 0, \end{cases} \quad (17)$$

for given rate coefficients $K_{v,open}$, $K_{v,closed} \in \mathbb{R}^+$.

## 7. OD heart modeling

The four-chamber heart model employed by the solver is inspired by Mynard and Smolich⁸ and includes source resistances, external constraints via the pericardium, and atrial-ventricular interactions via the septa. However, it is uniquely reformulated in differential equation form, as described in what follows.

For each chamber c, the time-varying free wall elastance $E_{c,FW}(t)$ is given by⁸

$$E_{c,FW}(t) = k_c \left( \frac{g_{c,1}(t)}{1 + g_{c,1}(t)} \right) \left( \frac{1}{1 + g_{c,1}(t)} \right) + E_{c,MIN}, \quad (18)$$

where $k_c \in \mathbb{R}$ is a scaling constant. The functions $g_{c,i}$, $i = 1, 2$ are given in terms of contraction onset time $t_{c,onset}$ by

$$g_{c,i}(t) = \left( \frac{t - t_{c,onset}}{\tau_{c,i}} \right)^{m_i},$$

where the constants $\tau_{c,1}$, $\tau_{c,2} \in \mathbb{R}$ are contraction and relaxation time offsets, respectively, and the constants $m_{c,1}$, $m_{c,2} \in \mathbb{R}$ are contraction and relaxation rates, respectively. The corresponding septal elastance is given in terms of free wall elastance of the chamber c (e.g., the left ventricle) and its contralateral chamber $c^*$ (e.g., the right ventricle) by⁸

$$E_{c,sep}(t) = \kappa_c E_{c,FW}(t) + \kappa_{c^*} E_{c^*,FW}(t) \quad (19)$$

for constants $\kappa_c$, $\kappa_{c^*} \in \mathbb{R}$. This leads to the native elastance given by

$$E_{c,nat}(t) = \frac{E_{c,FW}(t) E_{c,sep}(t)}{E_{c,FW}(t) + E_{c,sep}(t)} - \mu_{c,AV} Q_c(t), \quad (20)$$

where $Q_c(t)$ is the ventricular flow and $\mu_{c,AV} \in \mathbb{R}$ is an atrioventricular constant modeling piston-like effects.⁸ Defining a time-varying model of pressure from the pericardial sac as

$$P_{pc}(t) = K_{pc} \exp \left[ \frac{\sum_{c=1}^{4} V_c(t) + V_{myo} + V_{pcf} - V_{0,pc}}{\Phi_{pc}} \right] \quad (21)$$
for chamber volumes $V_c(t)$ and constants $V_{myo}$, $V_{pcf}$, $V_{0,pc}$, and $\Phi_{pc}$, the final pressure in each chamber is given by⁸

$$P_c(t) = P_{pc}(t) + E_{c,nat}(t)(V_c(t) - V_{c,0})(1 - K_{c,s}Q_c(t)) + \frac{E_{c,nat}(t)}{E_{c,sep}(t)}P_c^*(t). \quad (22)$$

Here, $K_{c,s}$ is a source flow resistance constant, $V_{c,0}$ is a given constant residual (zero pressure) chamber volume, and $P_c^*(t)$ is the pressure in the contralateral chamber [also determined by Eq. (22)]. Such a formulation for pressure has been employed by a number of closed-loop solvers.⁸,⁴⁴,⁴⁵ This work instead considers the time-derivative of Eq. (22), i.e.,

$$\frac{\partial P_c}{\partial t} = \frac{\partial}{\partial t} \left[ P_{pc}(t) + E_{c,nat}(t)(V_c(t) - V_{c,0})(1 - K_{c,s}Q_c(t)) + \frac{E_{c,nat}(t)}{E_{c,sep}(t)}P_c^*(t) \right], \quad (23)$$

which is solved as an initial value problem.

**Remark II.3.** In all other solvers, Eq. (22) is employed where, numerically, the volume $V_c(t + \Delta t)$ at the next timestep is calculated from the volume flow rate $Q = -dV/dt$ by a simple two-point (i.e., low-order) trapezoidal integration over the interval $[t, t + \Delta t]$ as

$$V_c(t + \Delta t) = V_c(t) - \frac{\Delta t}{2}(Q(t + \Delta t) - Q(t)).$$

In order to facilitate the high-order time integration described in Sec. III B, the solver introduced here employs Eq. (23), which can be expanded and algebraically rewritten to avoid explicit tracking and calculation of the volumes $V_c(t)$ (values of $dV/dt$ at $t$ in the differential form are replaced by the modeled chamber flow $dV/dt = -Q(t)$, including in the time derivative of the pericardial pressure $P_{pc}(t)$). Volume in each chamber can be recovered after evolving Eq. (23) to a time $t + \Delta t$ and substituting the new pressures back into Eq. (22) in order to construct a system to solve all chamber volumes simultaneously (since they are coupled via $P_{pc}(t)$ and the contralateral pressure $P_c^*(t)$) via a classical Newton-Raphson algorithm.

### C. Summary of geometric and physiological configurations

Figure 1 provides illustrative examples of the open- and closed-loop configurations (Secs. IV C and V, respectively) treated by the solver introduced in this work.

### III. SOLVER METHODOLOGY: FC AND HIGH-ORDER TIME INTEGRATION

This section describes the high-order Fourier continuation (FC) treatment²⁴⁻²⁶,³⁰ of spatial derivatives in the governing PDE system described by Eq. (1). Corresponding time integrations, both for the PDEs and for all 0D models described in Sec. II B, are similarly treated to high-order and described in Sec. III B.

### A. Spatial discretization: Accelerated FC(Gram)

Numerical treatment of the spatial variables and derivatives of Eq. (1) is achieved by considering point values $A(x_i)$ of a given

smooth function $A(x) : [0, 1] \to \mathbb{R}$ defined on a structured uniform discretization $x_i = i\Delta x$, $i = 0, ..., N - 1$, $\Delta x = 1/(N - 1)$. The FC method constructs a fast-converging interpolating trigonometric polynomial (Fourier series representation) $A_{cont} : [0, b] \to \mathbb{R}$ on a region $[0, b]$ that is slightly larger than the original physical domain of definition $[0, 1]$:

$$A_{cont}(x) = \sum_{k=-M}^{M} a_k \exp\left(\frac{2\pi ikx}{b}\right) \quad (24) \\ \text{s.t. } A_{cont}(x_i) = A(x_i), \quad i = 0, ..., N-1,$$

where $M = (N + N_{cont})/2$ is a full bandwidth parameter for a number of points $N_{cont}$ added to the original domain [such that $b = (N + N_{cont})\Delta x$]. The resulting FC function $A_{cont}$ extends the original function $A(x)$ into one that is discretely periodic, i.e., $A_{cont}(x)$ approximates the original $A(x)$ to very high-order in the original interval $[0, 1]$ but is approximately (and discretely) periodic on the slightly larger interval $[0, b]$, $b > 1$. Spatial derivatives of a PDE can then be produced by exact termwise differentiation of this series as

$$\frac{\partial A_{cont}}{\partial x}(x) = \sum_{k=-M}^{M} \left(\frac{2\pi ik}{b}\right) a_k \exp\left(\frac{2\pi ikx}{b}\right). \quad (25)$$

This ultimately provides the numerical derivatives of $A$ to high-order by restricting the domain of $\partial A_{cont}/\partial x$ to the original unit interval. Hence, the approximation rests in the construction of Eq. (24) from which the computation of the derivative Eq. (25) can be rapidly facilitated by the fast Fourier transform (FFT). Note here that $A$ has been defined on $[0, 1]$ without loss of generality (and can be extended to any interval via affine transformations³⁰). In the most straightforward intuitive manner,⁵⁸,⁵⁹ the coefficients $a_k$ of (24) can be found via the solution to the least squares problem given by

$$\min_{a_k} \sum_{i=0}^{N-1} |A_{cont}(x_i) - A(x_i)|^2 \quad (26)$$

through the Singular Value Decomposition (SVD). This can become rather costly for time-dependent solutions of complex boundary-valued PDEs such as those considered here (where each vascular segment requires applications of SVDs at each timestep).

An accelerated method,²⁵,²⁶,³⁰,³¹ known as FC(Gram), is applied by the solver introduced in this work and can circumvent such expense by employing small vectors of only a handful of function values near the left and right end points at $x = 0$ and $x = 1$ which can then be projected onto a Gram polynomial basis [whose continuations are pre-computed through solving the corresponding least squares problem Eq. (26) by variable high-precision SVD²⁶]. Only a subset of the given function values on small numbers $d_\ell$ and $d_r$ of matching points $\{x_0, ..., x_{d_\ell-1}\}$ and $\{x_{N-d_r}, ..., x_{N-1}\}$ contained in small subintervals on the left and right ends of the interval $[0, 1]$ are used to produce a discrete periodic extension of size $N_{cont}$. This is accomplished by projecting these end values onto a Gram polynomial basis of degree $d_\ell - 1$ (or $d_r - 1$) (producing a polynomial interpolant), where each basis element has its FC extensions precomputed with orthogonality enforced by the natural discrete scalar product defined by the discretization points. This effectively forms a “basis” of continuation functions with which to quickly and accurately extend any given function $A$ to
<figure><img src="image_15.png" /><figcaption>FIG. 1. Representative illustrations of the geometric and physiological configurations treated in this work. Left: the systemic arterial system. Yellow dots represent boundaries between vessel segments where branching conditions are employed. Right top: a single supplying artery (of Young's modulus E and thickness h corresponding to a pulse wave velocity c₀) with a one-sided Windkessel, used at the end of each terminal segment (e.g., in the left figure) for all open-loop benchmarks of Sec. IV C. Right middle: the corresponding two-sided Windkessel that is used to connect (possibly multiple) supplying arteries to draining veins for all closed-loop problems in Sec. V. Right bottom: a general overall schematic of such closed-loop setups, where the respective arterial and venous systems (including head/neck, coronary, and hepatic*) are represented by branching or merging networks that are similar to the left figure (detailed illustrations can be found elsewhere*). Arrows indicate heart chamber interactions modeled via the septum. MV: mitral valve, AV: aortic valve, TV: tricuspid valve, PV: pulmonary valve.</figcaption></figure>

provide a smooth transition from $A(x = 0)$ back to $A(x = 1)$ over the interval $[0, b]$.

The methodology for constructing such operators is briefly summarized as follows (for full technical details, the reader is referred to prior work$^{26}$). Defining the vectors of matching points for the left and right as

$$ 
\begin{aligned}
\mathbf{A}_{\ell} &= (A(x_0), A(x_1), \dots, A(x_{d_{\ell}-1}))^{\mathrm{T}}, \\
\mathbf{A}_{r} &= (A(x_{N-d_{\ell}}), A(x_{N-d_{\ell}+1}), \dots, A(x_{N-1}))^{\mathrm{T}},
\end{aligned}
\qquad (27) 
$$

the continuation operation can be expressed in a block matrix form as

$$ \mathbf{A}_{\text{cont}} = \left[ \mathbf{B}_{\ell} \mathbf{Q}_{\ell}^{\mathrm{T}} \mathbf{A}_{\ell} + \mathbf{B}_{r} \mathbf{Q}_{r}^{\mathrm{T}} \mathbf{A}_{r} \right] = \left[ \mathbf{C}_{\ell} \mathbf{A}_{\ell} + \mathbf{C}_{r} \mathbf{A}_{r} \right], \quad (28) $$

where $\mathbf{A} = (A(x_0), \dots, A(x_{N-1}))^T$ is a column vector containing the discrete point values of $A(x)$, $x \in [0, 1]$; $\mathbf{A}_{\text{cont}}$ is a vector of the $N + N_{\text{cont}}$

continued function values; $I$ is the $N \times N$ identity matrix; and $\mathbf{B}_{\ell}$, $\mathbf{B}_{r}$ contain the corresponding $N_{\text{cont}}$ values that blend (continue) to zero$^{26,31}$ from the left and the right (such that the sum of leftward and rightward continuations provides the necessary smooth transition). The columns of $\mathbf{Q} = \mathbf{Q}_{\ell}$, $\mathbf{Q}_{r}$ contain the $d = d_{\ell}$, $d_{r}$ point values of each element of the corresponding Gram polynomial basis, produced from a symbolic QR decomposition of a Vandermonde matrix $\mathbf{V}$ of monomials, i.e.,

$$ \mathbf{V} = \begin{pmatrix} 1 & x_0 & (x_0)^2 & \dots & (x_0)^{d-1} \\ 1 & x_1 & (x_1)^2 & \dots & (x_1)^{d-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{d-1} & (x_{d-1})^2 & \dots & (x_{d-1})^{d-1} \end{pmatrix} = \mathbf{Q}\mathbf{R}. \quad (29) $$

Each corresponding column of $\mathbf{B}_{\ell}$ (respectively $\mathbf{B}_{r}$) can then be found$^{26,31}$ by constructing an appropriate blend-to-zero function via an auxiliary band-limited trigonometric polynomial that is given by
$$f_{\text{aux}}(x) = \sum_{k=-M}^{M} a_k \exp\left(\frac{2\pi ikx}{(d + 2C + Z - 1)\Delta x}\right) \quad (30)$$

and that closely matches (in a least squares sense) the Gram polynomials (given by the columns $q_j$ of $Q$) and a number $Z$ of vanishing function values (where $M = (d + 2C + Z)/2$). The corresponding coefficients $a_k$ for each blending function are found by solving a minimization problem given by

$$\min_{\mathbf{a}=(a_{-M},...,a_M)^T} \left\| \mathbf{D}\mathbf{a} - \left( \begin{pmatrix} q_j \\ 0 \end{pmatrix} \right) \right\|_2, \quad (31)$$

where $\mathbf{D}$ is a matrix of function values of Eq. (30) evaluated at the matching and zero-valued intervals. Such a system is underdetermined and possibly ill-conditioned; hence, it is solved by a high-precision SVD (for its stability and accuracy in treating such problems$^{60}$). The computational cost is not carried over to the solver: such a (pre-computation) procedure for constructing $\mathbf{B}_\ell, \mathbf{B}_r$ is performed only once, with the resulting operators stored in a file to be loaded whenever the solver is invoked. For complete technical details on the overall FC(Gram) algorithm, the reader is referred elsewhere.$^{26}$

Figure 2 illustrates an example Fourier continuation of a non-periodic function. The resulting continued vector $\mathbf{f}_{\text{cont}}$ can be

interpreted as a set of discrete values of a smooth and periodic function that can be approximated to high-order via an FFT on an interval of size $(N + N_{\text{cont}})\Delta x$. Following others,$^{24-26,34}$ the hemodynamics solver introduced in this work employs values of $N_{\text{cont}} = 25, Z = 12$, and $d_\ell, d_r = 5$.

## B. Time discretization

### 1. Governing PDEs

Similarly to other FC-based solvers,$^{24,26,31,34}$ an explicit multi-step fourth-order Adams-Bashforth (AB4) method is employed for time marching of the PDEs in each vascular segment. Considering the spatially discretized governing system [Eq. (1)] in the form given by

$$\left[ \frac{\partial \mathbf{A}}{\partial t}(t) \quad \frac{\partial \mathbf{U}}{\partial t}(t) \right] = \left[ \mathbf{f}_1(t) \quad \mathbf{f}_2(t) \right], \quad (32)$$

the corresponding formulation for evolving from a time $t$ to a time $t + \Delta t$ at a uniform timestep $\Delta t$ is given by

$$\left[ \begin{array}{c} \mathbf{A}(t + \Delta t) \\ \mathbf{U}(t + \Delta t) \end{array} \right] = \left[ \begin{array}{c} \mathbf{A}(t) + \frac{\Delta t}{24}(55\mathbf{f}_1(t) - 59\mathbf{f}_1(t - \Delta t) + 37\mathbf{f}_1(t - 2\Delta t) - 9\mathbf{f}_1(t - 3\Delta t)) \\ \mathbf{U}(t) + \frac{\Delta t}{24}(55\mathbf{f}_2(t) - 59\mathbf{f}_2(t - \Delta t) + 37\mathbf{f}_2(t - 2\Delta t) - 9\mathbf{f}_2(t - 3\Delta t)) \end{array} \right]. \quad (33)$$

Such an explicit approach entails only one evaluation of the right-hand-sides $f_1, f_2$ (which, although of FFT-speed, are the most computationally expensive components of the overall solver). The boundary junctions and ODEs described in Sec. II B are applied as Dirichlet conditions, imposed after each timestep $t \to t + \Delta t$ by direct injection via the hyperbolic characteristics.

**Remark III.1.** Other high-order explicit schemes can also be used, such as the fourth-order Runge-Kutta (RK4) method.$^{61}$ Indeed, for FC-based methods, both RK4 and AB4 provide adequate regions of stability.$^{26,31}$ However, the timestep in both the open-loop and closed-loop configurations considered in this work is determined by the numerically stiff boundary ODEs of Sec. II B that govern the microvasculature, endorgan effects, as well as the coupled heart and valve models. As such, the subsequently required temporal discretization in the overall solver is observed to fall well within both regions of stability for the FC-based resolution of the PDE systems. Although RK4 and AB4 in this case (at the same limited timestep) are equivalently accurate,$^{25}$ RK4 is significantly more expensive, entailing four evaluations of the right-hand-sides (whose complexity here is dominated by the FFT). Additionally, spatial boundary value enforcement at intermediate RK steps may be problematic,$^{62,63}$ especially for such highly time-dependent (ODE-based) boundary conditions.

A classical spatial spectral filter is additionally employed at each timestep to control error growth in unresolved modes.$^{24-26,31}$ For Fourier coefficients $\hat{F}_k$ of a function $F(x)$, this is given by

$$\sum_{k=-\frac{N}{2}}^{\frac{N}{2}} \hat{F}_k \exp(ikx) \to \sum_{k=-\frac{N}{2}}^{\frac{N}{2}} \exp(-\alpha(2k/N)^p) \hat{F}_k \exp(ikx), \quad (34)$$

where $p \in \mathbb{Z}^+$ determines the rate of coefficient decay and where the $\alpha \in \mathbb{R}^+$ determines the level of suppression (e.g., the highest-frequency modes are multiplied by $e^{-\alpha x}$). For all numerical simulations presented in this work, a “light” choice of values that does not completely eliminate the highest-frequency modes, $(p, \alpha) = (8, -0.5 \cdot \ln(10^{-2}))$, is utilized. For further discussion on filtering practices for FC, the reader is referred elsewhere.$^{24,26,64}$

### 2. ODE-type boundary conditions

In other 1D hemodynamics solvers for the full circulation, low-order time integration is employed for the ODEs,$^{7,8,44,45}$ even when the spatial PDE scheme is of higher order.$^{14,44}$ In stark contrast, the fully high-order solver proposed in this work treats all ODEs involved in modeling the heart, valves, and vascular beds (Sec. II B) by a fourth-order explicit RK4$^{61}$ method. Since such ODEs are strictly
<figure><img src="image_16.png" /><figcaption>FIG. 2. An illustrative example of a Fourier continuation applied to the (discretized) non-periodic function $A(x) = (1 - x^2) \sin(2\pi x - \pi/2)$, where $x \in [0, 1]$. The matrix operators $C_t$, $C_r$ blend the left and right intervals $A_t$, $A_r$ to zero via a projection onto a basis of precomputed extensions based on an interpolating (Gram) polynomial basis. The new discrete function values $A_{\text{cont}}$ represent a periodic function on $[0, b]$.</figcaption></figure>

scalar-valued initial value problems, computation of the associated right-hand-sides is not expensive (cf. Remark III.1). Additionally, the valve model⁷ described in Sec. II B 6 utilizes a particular ODE for valve opening and another for valve closing [Eq. (17)], and hence the form of the right-hand-side can switch suddenly from one timestep to another (although the overall solution may be smooth; such a switch is facilitated by a threshold detector based on the continuous pressure difference across the valve⁷). This can be problematic for methods such as AB4, which require and utilize knowledge of the right-hand-sides evaluated at previous timesteps. Thus, for consistency, all 0D ODE-based models employ the RK4 scheme (taking sequential forward intermediate steps), which is given for a general form $\partial p/\partial t = f(p, x, t)$ by

$$p(t + \Delta t) = p(t) + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4), \quad (35)$$

where

$$\begin{align*} k_1 &= f(t, p(t)), \\ k_2 &= f\left(t + \frac{\Delta t}{2}, p(t) + \Delta t \frac{k_1}{2}\right), \\ k_3 &= f\left(t + \frac{\Delta t}{2}, p(t) + \Delta t \frac{k_2}{2}\right), \\ k_4 &= f(t + \Delta t, p(t) + \Delta t k_3). \end{align*} \quad (36)$$

In addition to the fourth-order nature of such temporal integration, a particular novelty lies in its use for also evolving the four heart chamber models presented in Sec. II B 7 (as opposed to the low-order volume extrapolation normally used in non-differential form, cf. Remark II.3).

### C. Application of boundary conditions

ODE-based boundary conditions at vessel inlets and outlets are applied in terms of incoming and outgoing characteristics [Eq. (6)]

by finding the corresponding cross-sectional area. For example, at the outlet of a vessel j with a corresponding characteristic impedance $Z_{j,out}$ connected to a vascular bed model, this is accomplished by finding the corresponding cross-sectional area $A_{j,out}$ and applying incoming characteristics at the outlet (where the outgoing characteristics, which come from within the spatial domain, are naturally evolved via the governing hyperbolic PDEs). That is, from a new value of pressure $P = P_{\text{art}}$ obtained at time $t + \Delta t$ by evolving Eq. (10), the corresponding flow at the outlet of vessel j is given by Eq. (12) as

$$\begin{align} Q_{j,out}(t + \Delta t) &= A_{j,out}(t + \Delta t)U_{j,out}(t + \Delta t) \nonumber \\ &= \frac{P_{j,out}(t + \Delta t) - P(t + \Delta t)}{Z_{j,out}}, \quad (37) \end{align}$$

which, replacing $U_{j,out}(t + \Delta t)$ with its form in terms of the outgoing characteristics via Eq. (6) and noting that $P_{j,out}(t + \Delta t) = P_{j,out}(A_{j,out}(t + \Delta t))$, can be expressed as

$$\begin{align} A_{j,out}(t + \Delta t) &= \left( w_{+,out} - \frac{4}{b} \left( c(A_{j,out}(t + \Delta t)) - c_0 \right) \right) \nonumber \\ &= \frac{P_{j,out}(A_{j,out}(t + \Delta t)) - P(t + \Delta t)}{Z_{j,out}}, \quad (38) \end{align}$$

where b is the vessel wall exponent in Eq. (2) evaluated at the outlet. Such a system can be solved iteratively by a classical Newton-Raphson method,⁴⁴,⁴⁹,⁵⁵ after which the incoming characteristics can be constructed and the boundary conditions can be applied through direct injection of values from Eq. (6). ODE-derived quantities at the venous side, at the valves, and at the heart are straightforwardly analogous. Similarly, junction values can be applied by solving the system described in Sec. II B 2 for area $A(t + \Delta t)$ and velocity $U(t + \Delta t)$ via Newton-Raphson and subsequently enforcing the corresponding incoming characteristics. For further details, the reader is referred elsewhere.⁸,⁵⁵

### IV. PERFORMANCE AND BENCHMARK STUDIES

This section includes a variety of validating studies attesting to high-order convergence, limited numerical dispersion, and physiological accuracy/relevancy of the FC-based hemodynamics solver and framework proposed in this contribution.

#### A. Error and convergence analysis

In order to verify numerical accuracy, convergence, and, ultimately, code implementation of the proposed FC-based solver for nonlinear hemodynamic wave propagation, the classical method of manufactured solutions (MMS) is employed.²¹,²²,²⁴⁻²⁶,⁴³,⁶⁵⁻⁶⁷ In such an exercise, smooth solutions are postulated where the corresponding right-hand-side forcing terms and boundary conditions are algebraically derived such that the proposed solutions are exact for the PDE system. In the following experiment, highly-dynamic traveling solutions to the area $A(x, t)$ and velocity $U(x, t)$ are proposed as

$$\begin{pmatrix} A(x, t) \\ U(x, t) \end{pmatrix} = \begin{pmatrix} \alpha_A \sin(k_{1,A}x - \omega_{1,A}t) \cos(k_{2,A}x - \omega_{2,A}t) + \beta_A \\ \alpha_U \sin(k_{1,U}x - \omega_{1,U}t) \cos(k_{2,U}x - \omega_{2,U}t) + \beta_U \end{pmatrix} \quad (39)$$
for {$\alpha_A, k_{1,A}, \omega_{1,A}, k_{2,A}, \omega_{2,A}, \beta_A$} = {2.5 \cdot 10^{-4}, 5\pi/4, 6\pi, 7\pi/4, 2\pi, 2.5 \cdot 10^{-3}}$ and {$\alpha_U, k_{1,U}, \omega, 1, U, k_{2,U}, \omega_{2,U}, \beta_U$} = {2.5 \cdot 10^{-1}, -5\pi/4, 4\pi, -7\pi/4, 6\pi, 2.5 \cdot 10^{-1}}. The considered domain is a vessel of length $l=100$ cm without external pressure constraints ($P_{\text{ext}}=0$ mm Hg), a reference pressure $P_0(x)=71$ mm Hg, and non-uniform material parameters corresponding to

$$c_0^2(x) = \frac{2}{3\rho} \left( k_1 \exp \left( k_2 \sqrt{\frac{A(x, 0)}{\pi}} \right) + k_3 \right) \quad (40)$$

for $k_1 = 3 \cdot 10^6$ g/(s²cm), $k_2 = -9$ cm⁻¹, and $k_3 = 33.7 \times 10^4$ g/(s²cm) (parameters that are inspired by those employed for systemic arteries in the full circulation⁸) and fluid density $\rho = 1.06$ g/cm³, viscosity $\mu = 4 \cdot 10^{-3}$ Pa s, and velocity profile order $\xi = 9$. This yields an average pulse wave velocity of $c_0 = 4.6$ m/s over the domain and an average hyperelastic exponent of $b = 4.16$, where the subsequent pressure, $P(x, t)$, is given by Eq. (2). Snapshots at various times of the solutions for both area $A(x, t)$ and velocity $U(x, t)$ are illustrated in the top row of Fig. 3. Using progressively larger discretization sizes of $N = 10, 20, 40, 80, 160$, and 320 points, the FC-based simulation is advanced for 40 000 timesteps in all cases at a time discretization corresponding to $\Delta t = 1.25 \times 10^{-4}$ s (which is chosen to be sufficiently small so that errors are dominated by the spatial discretization sizes). The maximum normalized absolute $L^\infty$ errors (i.e., over all space and all time) are presented in the bottom row of Fig. 3 for both area and velocity, where overlaid slopes demonstrate that the solver achieves an expected fifth-order accuracy based on the chosen FC parameters ($d_l, d_r = 5$ matching points).

## B. Numerical dispersion and stability analysis

Numerical dispersion and stability qualities for the proposed FC-based solver can be analyzed through test problems for which waves propagate over longer distances or times (which, for the former, is the equivalent of increasing the number of wavelengths in a domain over correspondingly increasing time intervals for a crest to travel from one end to the other²⁴,²⁶). Here, such computational experiments are applied to fluid parameters $\rho = 1.04$ g/cm³, $\mu = 1 \times 10^{-3}$ Pa s, and $\xi = 2$, with wall material parameters given by the same as above [Eq. (40)]. The considered traveling manufactured solution is given by

$$\begin{pmatrix} A(x, t) \\ U(x, t) \end{pmatrix} = \begin{pmatrix} \alpha_A \sin(2\pi(x - vt)) + \beta_A \\ \alpha_U \sin(2\pi(x - vt)) + \beta_U \end{pmatrix} \quad (41)$$

for a wave speed $v = 100$ cm/s, where $\alpha_A = 1$ cm², $\beta_A = 4$ cm², $\alpha_U = 100$ cm/s, and $\beta_U = 25$ cm/s. In order to test numerical dispersion characteristics, a sequence of increasingly larger domains are simulated, i.e., for $x \in [0, N_\lambda \cdot 100]$ cm, where the associated time intervals correspond to the time for a crest of a wave to pass from one end to the other, i.e., $t \in [0, N_\lambda]$ s. One can interpret $N_\lambda$ as a parameter representing the number of wavelengths in the domain. Figure 4 provides temporal snapshots of the solutions $A(x, t)$ and $U(x, t)$ for examples of $N_\lambda = 1$ and $N_\lambda = 10$. Although such solutions are periodic, the FC-based solver still continues the domain and creates a periodic extension over a larger interval.

Propagation over long distances (and hence times) is simulated by increasing the size of the domain through multiples defined by $N_\lambda \in \mathbb{Z}^+$. Choosing a small enough timestep so that errors are dominated by those arising from the spatial discretization, each wavelength of a domain is discretized by a fixed density of $d_\lambda = 15$,

<figure><img src="image_17.png" /><figcaption>FIG. 3. Top row: time snapshots along the spatial domain of area $A(x, t)$ (left) and velocity $U(x, t)$ (right) corresponding to the manufactured analytical solution given by Eq. (39). Bottom row: the normalized $L^\infty$ errors (i.e., maximum errors over all space and time) after 40 000 timesteps for area $A(x, t)$ (left) and velocity $U(x, t)$ (right) produced by the FC-based solver applied to successively finer discretizations. The overlaid dashed lines demonstrate the fifth-order convergence expected from the FC parameters employed in the proposed solver ($d_l, d_r = 5$ matching points).</figcaption></figure>
$d_\lambda = 30$, and $d_\lambda = 45$ points per wavelength (PPW). The maximum (normalized) $L^\infty$ errors over all space and time, as functions of increasing domain length $N_\lambda$, are presented in the top row of Fig. 5. For reference, low-order finite difference/element/volume methods are well-known to be (often linearly) dispersive or “pollutive”: errors increase as a function of domain size for fixed points per wavelength $d_\lambda^{23-26,68-70}$— necessitating the use of increasing numbers of PPW for increasing domain lengths or simulation times. In contrast, the FC-based solver incurs errors that remain effectively constant as the size of the domain grows by a factor of $N_\lambda$, i.e., the errors are independent of frequency for a fixed density of $d_\lambda$ PPW in the overall spatial domain (essentially no numerical dispersion errors). Similar results have been reported for other FC-based solvers, including in comparison to high-order Padé-like methods and various orders of finite difference or finite volume schemes.²⁴,²⁶,³¹,⁶⁴,⁷¹,⁷² Such results imply that the proposed FC-based solver is well-adapted for modeling hemodynamic wave propagation throughout the vasculature for large numbers of cycles.

Stability can be similarly studied by considering Eq. (41) for a fixed spatial domain size $N_\lambda$ over increasing numbers of temporal periods or cycles (time intervals), where, again, the length of a cycle is the time required for any one crest to travel through the entire physical domain. For example, the bottom row of Fig. 5 presents the maximum normalized $L^\infty$ errors as a function of the number of cycles for an $N_\lambda = 10$ wavelength solution [corresponding to a spatial domain length of size $N_\lambda \cdot 100$ cm for Eq. (41)], where each cycle of length 10 s is discretized by 750, 1 500, and 2 250 timesteps for $d_\lambda = 15, 30$, and 45 PPW, respectively (with up to 150 000, 300 000, and 450 000 total timesteps in the simulations, respectively). The maximum error in each case remains constant for arbitrarily long times and hence implies absolute long-term stability of the FC-based nonlinear hemodynamics solver.

### C. Benchmark cases

In what follows, physiological accuracy of the proposed FC-based methodology is established, and its numerical implementation further validated, by treating a variety of published benchmark problems presented by Boileau et al.,¹⁴ based on the configurations discussed in Xiao et al.⁶ Such problems are “open-loop” by nature (i.e., with one-sided vascular beds as discussed in Remark II.2) but enable direct comparisons between FC and solutions provided by various other methods supplied as supplementary material in Boileau et al.:¹⁴ the discontinuous Galerkin method (hereafter denoted as DCG), the locally conservative Galerkin (LCG), the finite element method (FEM), the finite volume method (FVM), the McCormack method (McC), and the simplified trapezium rule method (STM). In all examples,¹⁴ the exponent in the corresponding tube law is taken to be unity, i.e., $b(x) = 1 \forall x$ in Eq. (2). That is, throughout each vessel segment $i = 1, \dots, N_{\text{seg}}$, for $N_{\text{seg}}$ quantifying the total vascular segments assigned to a problem, material properties are given in terms of a length $L_i$, elastic modulus $E_i(x)$ and wall thickness $h_i(x)$ such that the corresponding pulse wave velocity is given by¹⁴

$$c_i(x, t) = \sqrt{\frac{2\sqrt{\pi E_i(x) h_i(x)}}{3\rho A_{d,i}(x)}} A_i(x, t)^{1/4}, \quad (x, t) \in [0, L_i] \times \mathbb{R}^+. \quad (42)$$

Following the reference work,¹⁴ Secs. IV C 2–IV C 5 (Benchmarks 2–5) compare the time signals at a location $x = x_j$ (usually the midpoint) of a 1D segment for pressure $P(x, t)$ and flow $Q(x, t) = A(x, t)U(x, t)$ with respect to reference solutions $P^*(x, t)$ and $Q^*(x, t)$ (defined as either a 3D solution or in vivo experimental data¹⁴). For the pressure, such comparisons are made via various proposed error metrics defined as follows:¹⁴

<figure><img src="image_18.png" /><figcaption>FIG. 4. Time snapshots of the area solution $A(x, t)$ (top row) and the velocity solution $U(x, t)$ (bottom row) given by Eq. (41) for domain lengths of size $N_\lambda = 1$ (left) and $N_\lambda = 10$ (right).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 5. Top row: the normalized $L^\infty$ errors [i.e., maximum errors over all space and time relative to the solution given by Eq. (41)] as a function of domain size $N_\lambda$ (corresponding to a physical length of $N_\lambda \cdot 100$ cm) for area $A(x, t)$ (left) and velocity $U(x, t)$ (right) over a complete cycle (of length $N_\lambda$ s). Errors remain effectively constant for arbitrarily long distances. Bottom row: The normalized $L^\infty$ errors as a function of temporal cycles for area $A(x, t)$ (left) and velocity $U(x, t)$ (right) in a vessel that is $N_j = 10$ wavelengths in size. Each cycle spans 750, 1 500, and 2 250 timesteps for $d_\lambda = 15, 30, \text{ and } 45$ PPW, respectively (i.e., up to 150 000, 300 000, and 450 000 total timesteps in the simulations, respectively). Errors remain effectively constant, demonstrating long-term stability.</figcaption></figure>

$$ \varepsilon_P^{\text{RMS}}(x_j) = \sqrt{\frac{1}{N_t} \sum_{k=1}^{N_t} \left( \frac{P(x_j, t_k) - P^*(x_j, t_k)}{P^*(x_j, t_k)} \right)^2}, $$

$$ \varepsilon_P^{\text{MAX}}(x_j) = \max_k \left| \frac{P(x_j, t_k) - P^*(x_j, t_k)}{P^*(x_j, t_k)} \right|, $$

$$ \varepsilon_P^{\text{SYS}}(x_j) = \frac{\max_k P(x_j, t_k) - \max_k P^*(x_j, t_k)}{\max_k P^*(x_j, t_k)}, $$

$$ \varepsilon_P^{\text{DIA}}(x_j) = \frac{\min_k P(x_j, t_k) - \min_k P^*(x_j, t_k)}{\min_k P^*(x_j, t_k)}. $$

Here, $\varepsilon_P^{\text{RMS}}(x_j)$, $\varepsilon_P^{\text{MAX}}(x_j)$, $\varepsilon_P^{\text{SYS}}(x_j)$, and $\varepsilon_P^{\text{DIA}}(x_j)$ are the root mean square, the maximum value (over a cardiac cycle), the maximum systolic error, and the maximum diastolic error, respectively. The corresponding flow, vessel radius difference from inlet to outlet, and pressure difference from inlet to outlet, are considered alternatively (in order to prevent a divide by zero) by the expressions given by¹⁴

$$ \varepsilon_F^{\text{RMS}}(x_j) = \sqrt{\frac{1}{N_t} \sum_{k=1}^{N_t} \left( \frac{F(x_j, t_k) - F^*(x_j, t_k)}{\max_k F^*(x_j, t_k)} \right)^2}, $$

$$ \varepsilon_F^{\text{MAX}}(x_j) = \max_k \left| \frac{F(x_j, t_k) - F^*(x_j, t_k)}{F^*(x_j, t_k)} \right|, $$

$$ \varepsilon_F^{\text{SYS}}(x_j) = \frac{\max_k F(x_j, t_k) - \max_k F^*(x_j, t_k)}{\max_k F^*(x_j, t_k)}, $$

$$ \varepsilon_F^{\text{DIA}}(x_j) = \frac{\min_k F(x_j, t_k) - \min_k F^*(x_j, t_k)}{\min_k F^*(x_j, t_k)}. $$

**Remark IV.1.** Employing the open benchmark data of Boileau et al.,¹⁴ all error tables in Secs. IVC2–IVC4 have been re-calculated with an additional digit of precision in order to verify the self-implementation of the error analysis metrics (the originals are reported up to two decimal places). In doing so, it appears that the DCG results are shifted by a single timestep of $\Delta t = 1$ ms in order to arrive at the same published errors (such a value corresponds to the provided sampling rate of all simulation data for the different numerical methods¹⁴). This shift is additionally evident in the published waveform curves.¹⁴ Indeed, in a simple test (not discussed here), re-shifting the DCG data gives values that are closer to all the other methods of the reference work. However, in order to be consistent with the published data, the original shift is maintained for the re-computed error metrics for DCG. In order to account for this and provide a more complete comparison with the solver proposed herein, the FC-based results in all the corresponding benchmark error tables are provided for both cases: with and without the same timestep shift as DCG. It should be noted that this does not alter any conclusions of the original reference work:¹⁴ results from all computations and re-computations are all within the same order of accuracy with respect to the reference data (which is either fully 3D or based on in vivo experimental data). This implies that all methods—regardless of the shift—are more-or-less equally competitive (save for a few exceptions that are evident in the figures of this section) in terms of adequately capturing the 1D reduced-order approximations of the full-order or real-world data. The principal differences are related to the theoretical order of accuracy and subsequent computational cost.

**Remark IV.2.** Computational cost is compared with a DCG solver, generously provided by the authors,¹⁶ based on the open-source Nektar++ libraries. All such comparisons between the FC solver and the DCG solver are based on simulations conducted on the same
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>FIG. 6. (Left) Solution snapshots for the pressure wave solution of the single-pulse case considering both inviscid ($\mu = 0$ mPa s, solid blue) and viscous blood ($\mu = 4$ mPa s, dashed red). The solid black line denotes the theoretical rate of decay for the viscous case [Eq. (44)]. (Right) Two snapshots of the inviscid case. $^1$Boileau et al.$^{14}$</figcaption></figure>

MacBook Air M1 2020 with 16GB RAM (Rosetta 2 emulation of x86_64 MATLAB 2023a$^{73}$ for FC; self-compiled Nektar++ for DCG also using Rosetta 2 emulation). For the multi-segment cases, the MATLAB implementation of the FC-solver is supplemented, partially, with an automatic conversion of a couple of functions (notably, those requiring large loops: the vessel segment PDE evolution and the vascular beds) into un-optimized C++-based pre-compiled MATLAB MEX files that are constructed via the versatile MATLAB Coder (this is indicated where applicable: “MATLAB Classic” or “MATLAB+Coder”).

## 1. Benchmark 1: Single pulse

As a first benchmark,$^{14}$ a simple Gaussian pulse is propagated through a single vessel ($N_{seg} = 1$) of length $L = 10$ m (several times larger than the “effective pulse width”). Uniform material parameters are considered:$^{14}$ wall thickness $h = 1.5$ mm, blood density $\rho = 1050$ kg/m$^3$, velocity profile order $\zeta = 9$, and Young’s modulus $E = 400$ kPa. Initial area and velocity are taken to be $A(x, 0) = \pi$ cm$^2$ and $U(x, 0) = 0$ m/s, respectively. (Such parameters yield a corresponding initial pulse wave velocity of $c_0 \approx 6.17$ m/s.) Diastolic and external pressure are taken to be zero, i.e., $P_0 = P_{ext} = 0$. The incident pulse is generated by prescribing the volume inflow rate at the inlet via the expression given by$^{14}$

$$Q_{\text{in}}(t) = 10^{-6} e^{-10^{4}(t-0.05)^{2}} \text{ m}^{3}/\text{s}, \qquad (43)$$

and the simulation is advanced to a final time of $t_f = 2$ s.

Figure 6 (left) presents solution snapshots, at various times $t$, of the pressure (normalized by peak amplitude) generated by the FC-solver applied to both inviscid blood ($\mu = 0$ mPa s) and viscous blood ($\mu = 4$ mPa s). Results demonstrate that the corresponding simulations capture the indicated amplitude decay of the peak pressure (or lack thereof in the inviscid case), which is given$^{14}$ theoretically by

$$P_{\text{peak,viscous}}(x) = P_{\text{peak,inviscid}}e^{-(\zeta+2)\pi\mu x/(\rho c_0 A_0)}. \qquad (44)$$

The right panel of Fig. 6 presents a zoom on two such snapshots in time for both solutions generated by FC as well as for those generated by other methods surveyed in published benchmarks$^{14}$—some of which incur numerical dispersion/diffusion (“pollution”) errors (e.g., LCG, FVM) or require fine discretizations (e.g., FEM, STM). Indeed, a comparison of the orders-of-accuracy and subsequent choices of temporal and spatial discretizations in Table I demonstrates the competitive qualities of an FC-based approach: stable resolution at the highest order in space and time; a reasonably large timestep; and the use of among the fewest spatial solution points (with only FVM being lower—albeit with one of the largest evident discrepancies/errors in Fig. 6, likely due to numerical dispersion errors).

<table><caption>TABLE I. Benchmark 1: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 6 (right). Boldface denotes values associated with the FC solver introduced in this work.</caption><thead><tr><th>Method</th><th>Δt</th><th>N</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>2.500 × 10<sup>-4</sup></td><td>1551</td><td>5th</td><td>4th</td><td>1.63<sup>b</sup> (MATLAB Classic)</td></tr><tr><td>DCG<sup>a</sup></td><td>1.000 × 10<sup>-4</sup></td><td>2500<sup>a</sup></td><td>5th<sup>a</sup></td><td>2nd<sup>a</sup></td><td>16.86<sup>b</sup> (C++)</td></tr><tr><td>LCG<sup>a</sup></td><td>1.000 × 10<sup>-5</sup></td><td>2000<sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>FEM<sup>a</sup></td><td>1.000 × 10<sup>-7</sup></td><td>4001<sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>1.457 × 10<sup>-3</sup></td><td>1000<sup>a</sup></td><td>4th<sup>a</sup></td><td>4th<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>2.500 × 10<sup>-4</sup></td><td>3201<sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>1.000 × 10<sup>-4</sup></td><td>10 001<sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

<footnote>$^a$Boileau et al.$^{14}$</footnote>

<footnote>$^b$Remark IV.2.</footnote>
**Remark IV.3.** The spatial discretization size *N* in Table I is defined here as the total number of resolved unknown points in space (i.e., degrees-of-freedom); its value is not provided directly in the text of Boileau et al.<sup>14</sup> (which instead reports a “characteristic element size” whose definition differs by method of discretization), but can be identified from element size and order (additionally verified by the simulation data files provided in the supplementary material of that paper<sup>14</sup>). Hence, it is used in Table I as a more indicative and consistent computational cost/resolution measure for the different element-based (e.g., DCG, LCG, FEM), cell-based (e.g., FVM), and stencil-based (e.g., FC, McC/FD, STM) methods considered.

**2. Benchmark 2: Common carotid artery**

A second benchmark<sup>14</sup> concerns the modeling of a single segment ($N_{seg} = 1$) representing the common carotid artery. Such a domain is of length $L = L_{seg} = 126$ mm and consists of uniform material parameters with a wall thickness $h = h_i(x) = 0.3$ mm and a Young's modulus $E = E_i(x) = 700$ kPa. Fluid parameters correspond to blood density $\rho = 1060$ kg/m<sup>3</sup>, blood viscosity $\mu = 4$ mPa·s, and velocity profile order $\zeta = 2$. Area and pressure at diastole are taken to be $A_0(x) = 0.28274$ cm<sup>2</sup> and $P_0(x) = 10.933$ kPa, respectively. Initial area and velocity are taken to be $A(x, 0) = A_d(x)$ and $U(x, 0) = 0$ m/s, respectively. External pressure is taken to be $P_{ext} = 0$ kPa, where the outlet is terminated by a three-element Windkessel with parameters $R_1 = 2.4875 \times 10^8$ Pa·s/m<sup>3</sup>, $C = 1.7529 \times 10^{-10}$ m<sup>3</sup>/Pa, $R_1 = 1.8697 \times 10^9$ Pa·s/m<sup>3</sup>, and $P_{out} = 0$ kPa. A flow profile $Q_{in}(t)$ (given by a closed-form expression provided in the supplementary material of Boileau et al.<sup>14</sup>) is prescribed at the inlet, and the simulation is advanced to a final time of $t_f = 10$ s. A total of $N = 5$ discretization points is employed in the FC solver.

Figure 7 presents the time evolution, for the last complete cardiac cycle (length $T_{per} = 1.1$ s), of FC-generated pressure and flow at the midpoint $x = x^*$ of the modeled carotid artery (overlaid with corresponding benchmark solutions<sup>14</sup>). Results demonstrate excellent agreement between all simulations. Indeed, when compared to reference 3D solutions of pressure $P^*(x^*, t)$, flow $Q^*(x^*, t)$, pressure difference

$\Delta P^*(t)$, and vessel radius difference $\Delta r^*(t)$, the various error metrics in Table IV of Appendix A provide a mutual validation of the FC solver with other approaches. The qualities of the FC-based approach in terms of order-of-accuracy (the highest), timestep (among the highest), and computational cost (competitive with DCG) is additionally demonstrated in Table X of Appendix B.

**3. Benchmark 3: Aortic bifurcation**

A third benchmark<sup>14</sup> concerns the modeling of three segments ($N_{seg} = 3$): an aorta and its bifurcation into iliac arteries. All model parameters can be found in the reference work.<sup>14</sup> A flow profile $Q_{in}(t)$, given by a closed-form expression provided in the supplementary material of Boileau et al.,<sup>14</sup> is prescribed at the inlet, and the simulation is advanced to a final time of $t_f = 30$ s. A total of $N = 15$ discretization points (five in each segment) is employed in the FC solver.

The top and bottom rows of Fig. 8 present the time evolution, for the last complete cardiac cycle (length $T_{per} = 1.1$ s), of FC-generated pressure and flow at the midpoint $x = x^*$ of the aorta and the bifurcated iliac arteries, respectively (overlaid with corresponding benchmark solutions<sup>14</sup>). Results demonstrate excellent agreement between all simulations. Only minimal differences can be observed between any two methods in the line plots. Indeed, when compared to reference 3D solutions of pressure $P^*(x^*, t)$, flow $Q^*(x^*, t)$, pressure difference $\Delta P^*(t)$, and vessel radius difference $\Delta r^*(t)$, the various error metrics in Tables V and VI of Appendix A provide a mutual validation of the FC solver with other approaches. The qualities of the FC-based approach in terms of order-of-accuracy (the highest) and computational cost (competitive with DCG) is additionally demonstrated in Table XI of Appendix B.

**4. Benchmark 4: 37-artery network**

A fourth benchmark<sup>14</sup> concerns the modeling of thirty-seven segments ($N_{seg} = 37$). All model parameters can be found in the reference work.<sup>14</sup> A flow profile $Q_{in}(t)$, given by a closed-form

<figure><img src="image_22.png" /><figcaption>FIG. 7. Time evolution at the midpoint of the common carotid artery for pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right). FC-generated solutions are overlaid with those generated by other methods. Solid lines are plotted at a resolution of $\Delta t = 1$ ms, and markers are plotted every twenty timesteps (for easier visualization). <sup>14</sup>Boileau et al.<sup>14</sup> "Without DCG shift (see Remark IV.1).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**FIG. 8.** Time evolution of pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right) at the midpoints of the aorta (top row) and the iliac artery (bottom row). FC-generated solutions are overlaid with those generated by other methods. Solid lines are plotted at a resolution of $\Delta t = 1$ ms, and markers are plotted every twenty timesteps (for easier visualization).

<sup>†</sup>Boileau *et al.*<sup>14</sup> *Without DCG shift* (see Remark IV.1).</figcaption></figure>

expression provided in the supplementary material of Boileau *et al.*,¹⁴ is prescribed at the inlet, and the simulation is advanced to a final time of $t_f = 15$ s. A spatial discretization of $\Delta x = 1$ cm is employed in the FC solver.

Figure 9 presents the time evolution, for the last complete cardiac cycle (length $T_{\text{per}} = 0.827$ s), of FC-generated pressure and flow at the midpoint $x = x^*$ of three randomly-selected segments: the left subclavian I, the right iliac-femoral II, and the splenic segment, respectively (overlaid with corresponding benchmark solutions¹⁴). Correspondence between all methods is mixed, although the majority agree well with FC. Indeed, when compared to experimental reference solutions of pressure $P^*(x^*, t)$, flow $Q^*(x^*, t)$, pressure difference $\Delta P^*(t)$, and vessel radius difference $\Delta r^*(t)$, the various error metrics in Tables VII–IX of Appendix A provide a mutual validation of the FC solver with other approaches (although FC, in general, has lower errors relative to the reference data for a number of metrics). The qualities of the FC-based approach in terms of order-of-accuracy (the highest) and computational cost (faster than DCG) is additionally demonstrated in Table XII of Appendix B.

## 5. Benchmark 5: 77-artery network (ADAN56)

A fifth benchmark¹⁴ concerns the modeling of seventy-seven segments ($N_{\text{seg}} = 77$). All model parameters can be found in the reference work.¹⁴ A flow profile $Q_{\text{in}}(t)$, given discretely (and interpolated in time) from the supplementary material of Boileau *et al.*,¹⁴ is prescribed at the inlet, and the simulation is advanced to a final time of $t_f = 10$ s. A spatial discretization of $\Delta x = 0.66$ cm is employed in the FC solver.

Figure 10 presents the time evolution, for the last complete cardiac cycle (length $T_{\text{per}} = 1$ s), of FC-generated pressure and flow at the midpoint $x = x^*$ of three randomly selected segments: the thoracic aorta III, the right internal carotid artery, and the right internal iliac artery, respectively (overlaid with corresponding benchmark solutions¹⁴). Most of the solvers agree well with each other and FC (except for STM). No reference solutions (and hence no error tables) are considered in the original work¹⁴ and hence are not reported here. The qualities of the FC-based approach in terms of order-of-accuracy (the highest) and computational cost (faster than DCG) are additionally demonstrated in Table XIII of Appendix B.

## 6. Summary of benchmarks

Figure 11 (left) presents the average error rank across all columns of Tables IV–IX in Appendix A for the seven considered methods, where a value of one is assigned to the lowest (best) absolute-value error and a value of seven is assigned to the highest (worst) absolute-value error (FC is considered here without the time shift discussed in Remark IV.1). The solver of this work performs, in general, better than the other methods, but it should be noted that errors are computed with respect to 3D or in vivo data (not the exact solution of the 1D equations), and so they cannot reasonably be expected to converge to the reference data. Indeed, the authors of the original benchmark paper¹⁴ suggest that such problems have been proposed mostly to validate method implementations. Figure 11 (right) additionally presents a summary of the numerical parameter tables (Table I and Tables X–XIII in Appendix B) in terms of the average order-of-accuracy in space and in time
<figure><img src="image_24.png" /><figcaption>**FIG. 9.** Time evolution of pressure $P(x*, t)$ (left) and flow $Q(x*, t)$ (right) at the midpoints of (from top to bottom) the left subclavian I, the right iliac-femoral II, and the splenic segment. FC-generated solutions are overlaid with those generated by other methods. Solid lines are plotted at a resolution of $\Delta t = 1$ ms, and markers are plotted every twenty timesteps (for easier visualization). †Boileau et al. ¹⁴ ²Without DCG shift (see Remark IV.1).</figcaption></figure>

applied by each of the methods to the benchmarks in Secs. IVC1-IVC5. Here, the qualities of the FC-based solver are evident: the highest order in both. Such accuracy, as demonstrated in Secs. IVB and V, enables robust and dispersionless resolution of closed-loop configurations that are valuable for simulating long-time behavior.

# V. FC-BASED SIMULATIONS OF CLOSED-LOOP SYSTEMS

## A. A simple circulation configuration with valves and coupled heart chambers

In order to verify the efficacy of the solver in treating closed-loop formulations and, in particular, the proposed ODE-based formulation of the fully-coupled heart (Sec. II B 7), a simplified configuration is considered.⁷ Such a model consists of one segment each that models each of the four major parts of the circulation: the systemic arterial system, the systemic venous system, the pulmonary arterial system, and the pulmonary venous system (with two vascular beds: one between the systemic segments and one between the

pulmonary). The model additionally considers four valves: aortic, pulmonary, tricuspid, and mitral. All geometric and material parameters are adopted directly from Mynard et al.⁷ However, in that contribution, a decoupled heart model is considered [i.e., without septal interactions between ventricles and atria of Eq. (22)] that also neglects the pericardium (i.e., no external pressure constraints on the heart). Hence, all heart, valve, and pericardium parameters are instead adopted from follow-up work by the same authors,⁸ which incorporates such coupling between the four heart chambers as well as the pericardial pressure source. The reader is reminded that the heart model herein is solved in a unique way through the FC solver by employing an ODE-based approach to facilitate high-order accuracy (instead of the usual low-order trapezoidal integration of volumes,⁷,⁸,⁴⁴ see Remark II.3).

With fluid parameters corresponding to $\rho = 1060$ kg/m³, $\mu = 3.5 \times 10^{-3}$ Pa s, and $\xi = 9$, the simulation is advanced to 40 cardiac cycles (at a heart rate of 75 beats per minute [BPM]) at a timestep of $\Delta t = 5 \times 10^{-5}$ s, where the multi-segment domain is discretized by $N = 35$ points on each of the two systemic segments and $N = 18$ on each of the pulmonary segments. Figure 12 presents
<figure><img src="image_25.png" /><figcaption>Benchmark 5: midpoint flow (ThoAortaIII)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Benchmark 5: midpoint pressure (RightICA)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Benchmark 5: midpoint flow (RightICA)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Benchmark 5: midpoint pressure (RightInternalIliac)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Benchmark 5: midpoint flow (RightInternalIliac)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**FIG. 10.** Time evolution of pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right) at the midpoint of (from top to bottom) the thoracic aorta III, the right internal carotid artery, and the right internal iliac artery. FC-generated solutions are overlaid with those generated by other methods. Solid lines are plotted at a resolution of $\Delta t = 1$ ms, and markers are plotted every twenty timesteps (for easier visualization). \*Boileau et al. <sup>14</sup> \*Without DCG shift (see Remark IV.1).</figcaption></figure>

the time evolution, for the first 15 cycles simulated by the FC solver, of pressure $P(x^*, t)$ and flow $Q(x^*, t)$ at the inlets of the systemic artery (connected to the aortic valve) and pulmonary artery (connected to the pulmonary valve). In both cases, convergence of the physiological waveforms as a function of cardiac cycle is achieved rather quickly. Table II additionally presents the integrated flow rate over a cycle (i.e., cardiac output) for the same locations at the 20th, 30th, and 40th cycles (beats), demonstrating conservation of flow as well as lack of compounding pollution errors over longer time intervals (up to 640 000 total timesteps). Indeed, due to the high-order and numerically dispersionless nature of the FC-based solver, confident and stable predictions of cardiac output (which is a well-known and clinically common measure of cardiovascular health) can be obtained for increasingly larger numbers of cycles, indicating a robustness toward simulating long-time hemodynamics that is useful for coupling with physiological feedback mechanisms.<sup>41,42</sup> To the best of the authors' knowledge, this has not been explicitly demonstrated for other open-loop or closed-loop solvers (which are also often run only until a periodic steady-state is achieved, defined by

when the beat-to-beat pressure variations are within 1%, often around ten cycles<sup>7,8</sup>). All such results additionally validate the proposed ODE-based heart formulation as well as its corresponding code implementation and numerical performance.

## B. A comprehensive circulation configuration

A full model is directly inspired by and adopted from the seminal work of Mynard and Smolich<sup>8</sup> (the predominant source for the model components described in Sec. II and elsewhere<sup>44,45</sup>). The configuration employed here entails 394 possibly non-uniform segments (i.e., with vessel tapering) representing the systemic and pulmonary systems, with a plethora of diverging and merging junctions/bifurcations/trifurcations and multi-input vascular beds (including the hepatic and coronary systems described in Secs. II B 4 and II B 5, respectively). All parameters are adopted from the same source,<sup>8</sup> including the tree structure, vessel geometric/material properties, vascular bed network, and heart/valve parameters. Fluid parameters
<figure><img src="image_31.png" /><figcaption>FIG. 11. Left: the average error rank across all columns from Tables IV–IX of Appendix A, where one is the lowest (best) absolute value error and seven is the highest (worst). Right: the average order-of-accuracy in space and in time invoked by each method across all the benchmarks of Secs. IV C 1–IV C 5. ¹Boileau et al.¹⁴ ²Without time shift (see Remark IV.1).</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>FIG. 12. Time evolution at the inlets of the systemic artery (top) and the pulmonary artery (bottom) for pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right).</figcaption></figure>

corresponding to $\rho = 1060 \text{ kg/m}^3$, $\mu = 3.5 \times 10^{-3} \text{ Pa s}$, and $\xi = 9$ are employed, with segment material properties determined according to Eq. (40). The simulation is advanced to 10 cardiac cycles at a heart rate of 70 BPM and at a timestep of $\Delta t =$

<table><caption>TABLE II. Simple circulation: cardiac output at different cardiac cycles (heartbeats) for the left ventricle (LV) and the right ventricle (RV). Values remain essentially constant as the number of cycles increases, implying robust numerical performance.</caption><thead><tr><th>Segment</th><th>Cycle 20 (l/min)</th><th>Cycle 30 (l/min)</th><th>Cycle 40 (l/min)</th></tr></thead><tbody><tr><td>LV</td><td>5.678</td><td>5.676</td><td>5.673</td></tr><tr><td>RV</td><td>5.679</td><td>5.677</td><td>5.674</td></tr></tbody></table>

$5 \times 10^{-5} \text{ s. Figure 13 presents the time evolution over the last complete cycle of pressure } P(x^*, t) \text{ and flow } Q(x^*, t) \text{ at the inlet of the left anterior descending artery, the hepatic portal vein, the coronary sinus, and the main pulmonary artery, respectively. Such figures illustrate that the FC-based solver is able to capture physiological morphologies.}$^{41,54,74,75}$

In order to demonstrate robustness to different ranges of important clinical parameters, Fig. 14 presents the aortic pressure and flow waveforms for a low heart rate of 40 BPM (bradycardia), a normal heart rate of 70 BPM (baseline), and a high heart rate of 120 BPM (tachycardia), which are, again, physiological.⁴¹,⁵⁴,⁷⁴,⁷⁵ The dashed line in each pressure waveform demarcates the location of the dirotic notch,²⁴ which represents the time of aortic valve closure and which
<figure><img src="image_33.png" /><figcaption>FIG. 13. Time evolution (baseline) of pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right) at the inlets of (from top to bottom) the left anterior descending artery, the hepatic portal vein, the coronary sinus, and the main pulmonary artery.</figcaption></figure>

separates the so-called systolic phase (when the valve is open) from the so-called diastolic phase (when the valve is closed). Each simulation employs identical physical parameters, and no numerical tuning is necessary to achieve such solutions.

Another important clinical parameter relates to left ventricular (LV) contractility (how well the heart muscles contract), which is used to assess LV function and overall cardiac performance.⁴¹,⁷⁵,⁷⁶ The direct, load-independent gold standard measure of LV contractility is the end-systolic elastance⁷⁵ ($\max_t E_{FW}$ in the heart model employed in this work), where abnormally low or high values indicate dysfunction. However, obtaining such measurements requires highly-invasive and impractical hemodynamics manipulation through a conductance catheter.⁷⁵ A more indirect measure that is often used is the ejection fraction (EF), which can be acquired non-invasively through ultrasound

imaging⁴¹ (echocardiography). EF is a surrogate for LV contractility that is used in clinical practice⁴¹ and is defined as the ratio of stroke volume (the amount of blood pumped out of the left ventricle during a cardiac cycle) and end-diastolic volume (the amount of blood left in the left ventricle immediately before the next contraction), where values of EF less than 0.4 are strongly indicative of heart failure.⁴¹,⁷⁶ Table III presents ejection fractions calculated from simulations of varying levels of LV contractility at a normal heart rate of 75 BPM (represented by nominal maximum free wall elastance $\max E_{FW}$, with a baseline corresponding to $\max E_{FW} = 2.5$ mm Hg/ml). FC-based solutions correctly capture the relationship between contractility and EF,⁷⁶ including at the baseline (yielding an expected healthy EF of $\approx 0.5$⁷⁶). It is notable here that the same numerical parameters are used for all simulations (i.e., no tuning) and that EF is not an input
<figure><img src="image_34.png" /><figcaption>FIG. 14. Time evolution of pressure $P(x^*, t)$ (top) and flow $Q(x^*, t)$ (bottom) at the inlet of the aortic root for heart rates (left to right) of 40, 70 (baseline), and 120 BPM. Dashed lines in the pressure waveforms indicate the respective locations of their dirotic notches (corresponding to the times of aortic valve closure, marking the ends of the systolic phases).</figcaption></figure>

parameter for such simulations, but an outcome of the global solution for the entire circulatory system (i.e., a final end-product⁴¹ which strongly depends on the input contractility parameter, $\max_t E_{FW}$, and weakly on preload and afterload⁷).

## VI. CONCLUSIONS

This work introduces a Fourier continuation-based (pseudo)spectral methodology and framework for the accurate and robust numerical modeling of hemodynamic wave propagation in the entire human circulatory system. The overall closed-loop solver, and its open-loop counterpart in Sec. IV C, is based on well-established theoretical formulations⁸,¹⁴ and is uniquely of high-

<table><caption>TABLE III. Full circulation: ejection fraction (EF) as a function of LV contractility [$\max_t E_{FW}(t)$], demonstrating the physiological accuracy of the solver over a wide range of (possibly extreme) contractilities. FC-based simulations correctly capture the relationship between contractility and EF, including at a healthy baseline corresponding to $\max_t E_{FW}(t)$ (yielding a healthy EF of ≈0.5).</caption><thead><tr><td>$\max_t E_{FW}(t)$ (mm Hg/ml)</td><td>EF</td></tr></thead><tbody><tr><td>0.7</td><td>0.167</td></tr><tr><td>1.5</td><td>0.365</td></tr><tr><td>2.0</td><td>0.451</td></tr><tr><td>2.5</td><td>0.515</td></tr><tr><td>3.5</td><td>0.600</td></tr><tr><td>5.0</td><td>0.674</td></tr></tbody></table>

order (in both space and time) for all constituent components (including the heart). Results demonstrate efficient resolution of the Navier–Stokes PDE systems that govern fluid–structure blood flow in vascular segments as well as the (non)linear ODE systems that govern the valves, end-organ effects, microcirculation, and the four chambers of the heart (including their coupling through the septum). A variety of presented numerical experiments, including a suite of single- and multi-segment open-loop benchmark cases¹⁴ (with both 3D and experimental reference data), attest to the high-order convergence, stability, physiological accuracy, and low computational cost of the proposed approach. The overall solver incurs no numerical dispersion or “pollution” errors, similarly to other FC-based solvers,²⁴⁻²⁶,³¹ offering a very competitive approach for long distance propagation through complex vasculature over many cardiac cycles. Such a fast, robust methodology may be well-suited for exploring wide ranges of physiological parameters and for coupling with 3D models (where, for the latter, such a solver can potentially provide more comprehensive 1D physiological boundary conditions that incorporate the entire closed-loop circulation as opposed to the more common lumped parameter 0D models⁴³).

## ACKNOWLEDGMENTS

Support is acknowledged from the Agence Nationale de la Recherche (ANR) Grant No. ANR-23-CE46-0008 (F.A.). This study was partially funded by the National Science Foundation CAREER Grant No. 2145890 (N.M.P.).
<table><caption>TABLE IV. Benchmark 2: a comparison of relative errors (%) with 3D reference data for FC and other benchmarked solvers of the results plotted in Fig. 7.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>RMS</sup><sub>Δr</sub></th><th>ε<sup>RMS</sup><sub>ΔP</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>Δr</sub></th><th>ε<sup>MAX</sup><sub>ΔP</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>Δr</sub></th><th>ε<sup>SYS</sup><sub>ΔP</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>Δr</sub></th><th>ε<sup>DIA</sup><sub>ΔP</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>0.237</td><td>0.283</td><td>1.007</td><td>4.270</td><td>0.509</td><td>1.119</td><td>1.922</td><td>15.898</td><td>-0.260</td><td>-0.525</td><td>-1.620</td><td>-15.138</td><td>0.269</td><td>0.260</td><td>0.062</td><td>4.982</td></tr><tr><td>FC<sup>c</sup></td><td>0.186</td><td>0.385</td><td>0.941</td><td>4.431</td><td>0.298</td><td>1.205</td><td>1.680</td><td>16.644</td><td>-0.258</td><td>-0.525</td><td>-1.616</td><td>-15.138</td><td>0.272</td><td>0.261</td><td>0.066</td><td>4.982</td></tr><tr><td>DCG<sup>a</sup></td><td>0.198</td><td>0.393</td><td>0.956</td><td>4.418</td><td>0.343</td><td>1.230</td><td>1.749</td><td>16.693</td><td>-0.271</td><td>-0.547</td><td>-1.646</td><td>-15.046</td><td>0.278</td><td>0.247</td><td>0.080</td><td>4.992</td></tr><tr><td>LCG<sup>a</sup></td><td>0.242</td><td>0.256</td><td>0.999</td><td>4.346</td><td>0.517</td><td>1.015</td><td>1.926</td><td>16.172</td><td>-0.269</td><td>-0.456</td><td>-1.626</td><td>-15.373</td><td>0.282</td><td>0.221</td><td>0.105</td><td>4.761</td></tr><tr><td>FEM<sup>a</sup></td><td>0.235</td><td>0.284</td><td>0.996</td><td>4.269</td><td>0.497</td><td>1.120</td><td>1.880</td><td>15.877</td><td>-0.262</td><td>-0.541</td><td>-1.614</td><td>-15.160</td><td>0.271</td><td>0.264</td><td>0.075</td><td>5.031</td></tr><tr><td>FVM<sup>a</sup></td><td>0.280</td><td>0.287</td><td>1.068</td><td>4.210</td><td>0.660</td><td>1.070</td><td>2.277</td><td>15.611</td><td>-0.285</td><td>-0.575</td><td>-1.682</td><td>-14.825</td><td>0.290</td><td>0.249</td><td>0.111</td><td>4.908</td></tr><tr><td>McC<sup>a</sup></td><td>0.242</td><td>0.302</td><td>1.007</td><td>4.465</td><td>0.516</td><td>1.200</td><td>1.928</td><td>16.576</td><td>-0.266</td><td>-0.621</td><td>-1.627</td><td>-15.733</td><td>0.277</td><td>0.234</td><td>0.086</td><td>4.447</td></tr><tr><td>STM<sup>a</sup></td><td>0.242</td><td>0.288</td><td>0.978</td><td>4.473</td><td>0.525</td><td>1.130</td><td>1.961</td><td>16.900</td><td>-0.265</td><td>-0.547</td><td>-1.630</td><td>-15.764</td><td>0.269</td><td>0.268</td><td>0.077</td><td>4.855</td></tr></tbody></table>

<table><caption>TABLE V. Benchmark 3 (aorta): a comparison of relative errors (%) with 3D reference data for FC and other benchmarked solvers of the results plotted in the top row of Fig. 8.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>RMS</sup><sub>Δr</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>Δr</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>Δr</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>Δr</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>0.364</td><td>0.833</td><td>2.472</td><td>0.624</td><td>2.512</td><td>3.993</td><td>-0.517</td><td>-2.510</td><td>-3.928</td><td>0.485</td><td>1.160</td><td>-1.405</td></tr><tr><td>FC<sup>c</sup></td><td>0.401</td><td>0.852</td><td>2.489</td><td>0.697</td><td>2.593</td><td>3.942</td><td>-0.517</td><td>-2.510</td><td>-3.928</td><td>0.485</td><td>1.160</td><td>-1.405</td></tr><tr><td>DCG<sup>a</sup></td><td>0.368</td><td>0.894</td><td>2.476</td><td>0.606</td><td>2.512</td><td>4.045</td><td>-0.532</td><td>-2.467</td><td>-3.959</td><td>0.499</td><td>1.157</td><td>-1.388</td></tr><tr><td>LCG<sup>a</sup></td><td>0.388</td><td>0.712</td><td>2.426</td><td>0.673</td><td>2.236</td><td>3.898</td><td>-0.539</td><td>-2.227</td><td>-3.861</td><td>0.511</td><td>1.029</td><td>-1.370</td></tr><tr><td>FEM<sup>a</sup></td><td>0.368</td><td>0.820</td><td>2.546</td><td>0.649</td><td>2.530</td><td>4.035</td><td>-0.514</td><td>-2.516</td><td>-3.991</td><td>0.473</td><td>1.158</td><td>-1.488</td></tr><tr><td>FVM<sup>a</sup></td><td>0.404</td><td>0.812</td><td>2.486</td><td>0.676</td><td>2.442</td><td>4.012</td><td>-0.551</td><td>-2.442</td><td>-3.996</td><td>0.531</td><td>1.155</td><td>-1.345</td></tr><tr><td>McC<sup>a</sup></td><td>0.378</td><td>0.832</td><td>2.484</td><td>0.663</td><td>2.561</td><td>3.980</td><td>-0.526</td><td>-2.554</td><td>-3.942</td><td>0.487</td><td>1.171</td><td>-1.399</td></tr><tr><td>STM<sup>a</sup></td><td>0.376</td><td>0.823</td><td>2.446</td><td>0.655</td><td>2.523</td><td>3.937</td><td>-0.514</td><td>-2.522</td><td>-3.927</td><td>0.481</td><td>1.168</td><td>-1.406</td></tr></tbody></table>

<table><caption>TABLE VI. Benchmark 3 (iliac): a comparison of relative errors (%) with 3D reference data for FC and other benchmarked solvers of the results plotted in the bottom row of Fig. 8.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>RMS</sup><sub>Δr</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>Δr</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>Δr</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>Δr</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>0.442</td><td>0.739</td><td>4.282</td><td>0.837</td><td>1.951</td><td>7.273</td><td>-0.820</td><td>-1.788</td><td>-7.217</td><td>0.536</td><td>1.151</td><td>-2.375</td></tr><tr><td>FC<sup>c</sup></td><td>0.463</td><td>0.586</td><td>4.287</td><td>0.893</td><td>2.090</td><td>7.381</td><td>-0.819</td><td>-1.788</td><td>-7.217</td><td>0.536</td><td>1.151</td><td>-2.375</td></tr><tr><td>DCG<sup>a</sup></td><td>0.452</td><td>0.920</td><td>4.292</td><td>0.844</td><td>2.328</td><td>7.295</td><td>-0.831</td><td>-1.684</td><td>-7.253</td><td>0.547</td><td>1.168</td><td>-2.361</td></tr><tr><td>LCG<sup>a</sup></td><td>0.454</td><td>0.659</td><td>4.303</td><td>0.864</td><td>2.003</td><td>7.364</td><td>-0.830</td><td>-1.836</td><td>-7.257</td><td>0.555</td><td>1.191</td><td>-2.368</td></tr><tr><td>FEM<sup>a</sup></td><td>0.441</td><td>0.647</td><td>4.349</td><td>0.840</td><td>2.013</td><td>7.370</td><td>-0.811</td><td>-1.811</td><td>-7.275</td><td>0.535</td><td>1.166</td><td>-2.447</td></tr><tr><td>FVM<sup>a</sup></td><td>0.471</td><td>0.675</td><td>4.311</td><td>0.920</td><td>1.803</td><td>7.475</td><td>-0.854</td><td>-1.575</td><td>-7.308</td><td>0.573</td><td>1.184</td><td>-2.349</td></tr><tr><td>McC<sup>a</sup></td><td>0.447</td><td>0.665</td><td>4.305</td><td>0.855</td><td>2.013</td><td>7.346</td><td>-0.825</td><td>-1.886</td><td>-7.246</td><td>0.544</td><td>1.192</td><td>-2.382</td></tr><tr><td>STM<sup>a</sup></td><td>0.443</td><td>0.662</td><td>4.257</td><td>0.848</td><td>2.016</td><td>7.337</td><td>-0.814</td><td>-1.794</td><td>-7.222</td><td>0.537</td><td>1.176</td><td>-2.410</td></tr></tbody></table>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>
<table><caption>TABLE VII. Benchmark 4 (left subclavian I): a comparison of relative errors (%) with experimental reference data for FC and other benchmarked solvers of the results plotted in the top row of Fig. 9.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>2.715</td><td>14.201</td><td>5.694</td><td>38.076</td><td>-3.534</td><td>-1.943</td><td>-3.343</td><td>-11.358</td></tr><tr><td>FC<sup>c</sup></td><td>2.673</td><td>14.202</td><td>5.491</td><td>37.826</td><td>-3.534</td><td>-1.943</td><td>-3.255</td><td>-11.358</td></tr><tr><td>DCG<sup>a</sup></td><td>3.093</td><td>14.324</td><td>6.131</td><td>38.392</td><td>-3.920</td><td>-2.115</td><td>-3.974</td><td>-11.568</td></tr><tr><td>LCG<sup>a</sup></td><td>3.125</td><td>13.884</td><td>6.126</td><td>34.758</td><td>-3.376</td><td>-2.922</td><td>-4.711</td><td>-11.326</td></tr><tr><td>FEM<sup>a</sup></td><td>3.054</td><td>14.188</td><td>6.096</td><td>37.891</td><td>-3.913</td><td>-2.534</td><td>-3.719</td><td>-10.858</td></tr><tr><td>FVM<sup>a</sup></td><td>3.124</td><td>14.260</td><td>6.293</td><td>38.124</td><td>-3.961</td><td>-2.248</td><td>-3.927</td><td>-11.052</td></tr><tr><td>McC<sup>a</sup></td><td>3.108</td><td>14.326</td><td>6.351</td><td>38.248</td><td>-4.094</td><td>-2.084</td><td>-3.692</td><td>-10.920</td></tr><tr><td>STM<sup>a</sup></td><td>3.108</td><td>14.464</td><td>6.131</td><td>38.965</td><td>-3.800</td><td>-1.184</td><td>-3.974</td><td>-10.833</td></tr></tbody></table>

<table><caption>TABLE VIII. Benchmark 4 (right iliac-femoral II): a comparison of relative errors (%) with experimental reference data for FC and other benchmarked solvers of the results plotted in the middle row of Fig. 9.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>3.411</td><td>24.460</td><td>8.449</td><td>62.056</td><td>-1.866</td><td>51.922</td><td>-3.637</td><td>-39.799</td></tr><tr><td>FC<sup>c</sup></td><td>3.456</td><td>24.715</td><td>8.210</td><td>62.871</td><td>-1.866</td><td>51.922</td><td>-3.637</td><td>-39.799</td></tr><tr><td>DCG<sup>a</sup></td><td>3.819</td><td>24.480</td><td>8.999</td><td>59.633</td><td>-2.205</td><td>50.203</td><td>-4.133</td><td>-39.719</td></tr><tr><td>LCG<sup>a</sup></td><td>3.975</td><td>24.160</td><td>9.685</td><td>61.143</td><td>-2.191</td><td>49.226</td><td>-5.575</td><td>-36.922</td></tr><tr><td>FEM<sup>a</sup></td><td>3.694</td><td>23.888</td><td>8.656</td><td>59.277</td><td>-2.316</td><td>49.995</td><td>-3.668</td><td>-37.176</td></tr><tr><td>FVM<sup>a</sup></td><td>3.750</td><td>24.178</td><td>9.085</td><td>60.613</td><td>-2.263</td><td>51.424</td><td>-4.037</td><td>-39.062</td></tr><tr><td>McC<sup>a</sup></td><td>3.651</td><td>24.788</td><td>8.946</td><td>61.610</td><td>-2.334</td><td>52.408</td><td>-3.733</td><td>-40.206</td></tr><tr><td>STM<sup>a</sup></td><td>3.749</td><td>24.106</td><td>9.317</td><td>60.408</td><td>-2.373</td><td>51.295</td><td>-4.171</td><td>-38.256</td></tr></tbody></table>

<table><caption>TABLE IX. Benchmark 4 (splenic): a comparison of relative errors (%) with experimental reference data for FC and other benchmarked solvers of the results plotted in the bottom row of Fig. 9.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>2.265</td><td>9.317</td><td>5.655</td><td>23.401</td><td>-0.429</td><td>0.234</td><td>0.292</td><td>-5.828</td></tr><tr><td>FC<sup>c</sup></td><td>2.271</td><td>9.435</td><td>5.689</td><td>23.556</td><td>-0.429</td><td>0.234</td><td>0.153</td><td>-5.828</td></tr><tr><td>DCG<sup>a</sup></td><td>2.345</td><td>9.433</td><td>6.254</td><td>23.740</td><td>-0.792</td><td>0.089</td><td>-0.275</td><td>-5.975</td></tr><tr><td>LCG<sup>a</sup></td><td>2.336</td><td>9.787</td><td>6.228</td><td>22.817</td><td>0.092</td><td>3.586</td><td>-0.766</td><td>-6.753</td></tr><tr><td>FEM<sup>a</sup></td><td>2.280</td><td>9.016</td><td>5.945</td><td>22.832</td><td>-0.900</td><td>-1.318</td><td>-0.097</td><td>-5.458</td></tr><tr><td>FVM<sup>a</sup></td><td>2.331</td><td>9.212</td><td>6.116</td><td>23.199</td><td>-0.872</td><td>-0.583</td><td>-0.171</td><td>-5.948</td></tr><tr><td>McC<sup>a</sup></td><td>2.217</td><td>9.030</td><td>5.566</td><td>23.305</td><td>-0.968</td><td>-1.474</td><td>-0.316</td><td>-5.215</td></tr><tr><td>STM<sup>a</sup></td><td>2.359</td><td>9.558</td><td>5.961</td><td>24.468</td><td>-0.378</td><td>1.047</td><td>-0.168</td><td>-6.398</td></tr></tbody></table>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>
## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Faisal Amlani:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology

(equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Niema M. Pahlevan:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

<table><caption>TABLE X. Benchmark 2: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 7 and the errors calculated in Table IV.</caption><thead><tr><th>Method</th><th>Δt</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>1.00 × 10<sup>-3</sup></td><td>5th</td><td>4th</td><td>0.35<sup>b</sup> (MATLAB Classic)</td></tr><tr><td>DCG<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>5th<sup>a</sup></td><td>2nd<sup>a</sup></td><td>0.31<sup>b</sup> (C++)</td></tr><tr><td>LCG<sup>a</sup></td><td>2.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>FEM<sup>a</sup></td><td>1.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>1st<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>3.27 × 10<sup>-3</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>4.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>1st<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

<table><caption>TABLE XI. Benchmark 3: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 8 and the errors calculated in Tables V and VI for the aorta and the iliac arteries, respectively.</caption><thead><tr><th>Method</th><th>Δt</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>5.00 × 10<sup>-4</sup></td><td>5th</td><td>4th</td><td>1.93<sup>b</sup> (MATLAB+Coder<sup>b</sup>)</td></tr><tr><td>DCG<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>3rd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>1.40<sup>b</sup> (C++)</td></tr><tr><td>FEM<sup>a</sup></td><td>2.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>1st<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>LCG<sup>a</sup></td><td>1.99 × 10<sup>-3</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>5.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>1st<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

<table><caption>TABLE XII. Benchmark 4: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 9 and the errors calculated in Tables VII–IX.</caption><thead><tr><th>Method</th><th>Δt</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>5 × 10<sup>-5</sup></td><td>5th</td><td>4th</td><td>38.56<sup>b</sup> (MATLAB+Coder<sup>b</sup>)</td></tr><tr><td>DCG<sup>a</sup></td><td>1 × 10<sup>-4</sup></td><td>3rd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>108.25<sup>b</sup> (C++)</td></tr><tr><td>FEM<sup>a</sup></td><td>1 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>1st<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>7.6 × 10<sup>-4</sup></td><td>3rd<sup>a</sup></td><td>3rd<sup>a</sup></td><td>Not available</td></tr><tr><td>LCG<sup>a</sup></td><td>1 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>2.5 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>1 × 10<sup>-3</sup></td><td>1st<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Remark IV.2.</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Remark IV.2.</footnote>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Remark IV.2.</footnote>
<table><caption>TABLE XIII. Benchmark 5: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 10.</caption><thead><tr><th>Method</th><th>Δt</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>2.5 × 10<sup>-5</sup></td><td>5th</td><td>4th</td><td>95.47<sup>b</sup><br>(MATLAB+Coder<sup>b</sup>)</td></tr><tr><td>DCG<sup>a</sup></td><td>1 × 10<sup>-4</sup><sup>a</sup></td><td>3rd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>168.52<sup>b</sup> (C++)</td></tr><tr><td>FEM<sup>a</sup></td><td>1 × 10<sup>-3</sup><sup>a</sup></td><td>2nd<sup>a</sup></td><td>1st<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>5.9 × 10<sup>-4</sup><sup>a</sup></td><td>3rd<sup>a</sup></td><td>2rd<sup>a</sup></td><td>Not available</td></tr><tr><td>LCG<sup>a</sup></td><td>1 × 10<sup>-5</sup><sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>5 × 10<sup>-4</sup><sup>a</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>5 × 10<sup>-4</sup><sup>a</sup></td><td>1st<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

## DATA AVAILABILITY

The data that support the findings of this study are available from the authors upon reasonable request.

## APPENDIX A: TABLES OF ERRORS RELATIVE TO 3D OR IN VIVO DATA FOR THE OPEN-LOOP BENCHMARK PROBLEMS

Tables IV–IX report the errors (defined in the introduction of Sec. IV C) corresponding to the figures in the benchmarks of Secs. IV C 2–IV C 4.

## APPENDIX B: TABLES OF NUMERICAL PARAMETERS EMPLOYED BY EACH METHOD FOR THE OPEN-LOOP BENCHMARK PROBLEMS

Tables X–XIII report the numerical parameters utilized by each solver for the benchmarks of Secs. IV C 2–IV C 5, respectively.

## APPENDIX C: AN ADDITIONAL BENCHMARK (UPPER THORACIC AORTA)

An additional benchmark¹⁴ in the same spirit of Sec. IV C 2 similarly concerns the modeling of a single segment ($N_{seg} = 1$), but here representing the *upper thoracic aorta*. Such a domain is of length $L = L_{seg} = 24.137$ cm and consists of uniform material parameters with a wall thickness $h = h_i(x) = 1.2$ mm and a Young's modulus $E = E_i(x) = 400$ kPa. Fluid parameters correspond to blood density $\rho = 1060$ kg/m³, blood viscosity $\mu = 4$ mPa s, and velocity profile order $\zeta = 9$. Area and pressure at diastole are taken to be $A_d(x) = 4.5239$ cm² and $P_d(x) = 9.466$ kPa, respectively. Initial area and velocity are taken to be $A(x, 0) = A_d(x)$ and $U(x, 0) = 0$ m/s, respectively. External pressure is $P_{ext}(x) = 0$ kPa, where the outlet is terminated by a three-element Windkessel with parameters $R_1 = 1.1752 \times 10^7$ Pa s/m³, $C = 1.0163 \times 10^{-8}$ m³/Pa, $R_1 = 1.1167 \times 10^8$ Pa s/m³, and $P_{out} = 0$ kPa. A flow profile $Q_{in}(t)$, given by a closed-form expression provided in the supplementary material of Boileau et al.,¹⁴ is prescribed at the inlet, and the simulation is advanced to a final time of $t_f = 20$ s. A total of $N = 13$ discretization points is employed in the FC solver.

Figure 15 presents the time evolution, for the last complete cardiac cycle (length $T_{per} = 0.955$ s), of FC-generated pressure and flow at the midpoint $x = x^*$ of the modeled carotid artery (overlaid with corresponding benchmark solutions¹⁴). Results demonstrate excellent agreement between all simulations. Only minimal differences can be observed between any two methods in the line plots. Indeed, when compared to reference 3D solutions of pressure $P^*(x^*, t)$, flow $Q^*(x^*, t)$, pressure difference $\Delta P^*(t)$, and vessel radius difference $\Delta r^*(t)$, the various error metrics in Table XIV provide a mutual validation of the FC solver with other approaches. The qualities of the FC-based approach in terms of order-of-accuracy (the highest), timestep (among the highest), and computational cost (faster than DCG) is additionally demonstrated in Table XV.

<figure><img src="image_35.png" /><figcaption>FIG. 15. Time evolution at the midpoint of the upper thoracic aorta for pressure $P(x^*, t)$ (left) and flow $Q(x^*, t)$ (right). FC-generated solutions are overlaid with those generated by other methods. Solid lines are plotted at a resolution of $\Delta t = 1$ ms, and markers are plotted every twenty timesteps (for easier visualization). <sup>a</sup>Boileau et al.¹⁴ <sup>b</sup>Without DCG shift (see Remark IV.1).</figcaption></figure>

<footnote><sup>a</sup>Boileau et al.¹⁴</footnote>

<footnote><sup>b</sup>Remark IV.2.</footnote>
<table><caption>TABLE XIV. Additional benchmark: a comparison of relative errors (%) with 3D reference data for FC and other benchmarked solvers of the results plotted in Fig. 15.</caption><thead><tr><th>Method</th><th>ε<sup>RMS</sup><sub>P</sub></th><th>ε<sup>RMS</sup><sub>Q</sub></th><th>ε<sup>RMS</sup><sub>Δr</sub></th><th>ε<sup>RMS</sup><sub>ΔP</sub></th><th>ε<sup>MAX</sup><sub>P</sub></th><th>ε<sup>MAX</sup><sub>Q</sub></th><th>ε<sup>MAX</sup><sub>Δr</sub></th><th>ε<sup>MAX</sup><sub>ΔP</sub></th><th>ε<sup>SYS</sup><sub>P</sub></th><th>ε<sup>SYS</sup><sub>Q</sub></th><th>ε<sup>SYS</sup><sub>Δr</sub></th><th>ε<sup>SYS</sup><sub>ΔP</sub></th><th>ε<sup>DIA</sup><sub>P</sub></th><th>ε<sup>DIA</sup><sub>Q</sub></th><th>ε<sup>DIA</sup><sub>Δr</sub></th><th>ε<sup>DIA</sup><sub>ΔP</sub></th></tr></thead><tbody><tr><td>FC<sup>b</sup></td><td>1.197</td><td>2.199</td><td>2.387</td><td>7.246</td><td>3.598</td><td>7.256</td><td>7.534</td><td>29.427</td><td>-0.310</td><td>-5.342</td><td>-1.921</td><td>-8.798</td><td>1.104</td><td>2.458</td><td>2.201</td><td>7.183</td></tr><tr><td>FC<sup>c</sup></td><td>1.190</td><td>2.527</td><td>2.354</td><td>7.713</td><td>3.882</td><td>9.316</td><td>8.135</td><td>31.961</td><td>-0.310</td><td>-5.342</td><td>-1.921</td><td>-8.798</td><td>1.104</td><td>2.458</td><td>2.201</td><td>7.183</td></tr><tr><td>DCG<sup>a</sup></td><td>1.088</td><td>2.546</td><td>2.316</td><td>7.699</td><td>3.544</td><td>9.134</td><td>7.799</td><td>31.827</td><td>-0.575</td><td>-5.558</td><td>-2.541</td><td>-8.916</td><td>0.849</td><td>2.751</td><td>1.854</td><td>7.322</td></tr><tr><td>LCG<sup>a</sup></td><td>1.208</td><td>2.191</td><td>2.405</td><td>7.224</td><td>3.552</td><td>7.200</td><td>7.674</td><td>29.308</td><td>-0.317</td><td>-5.356</td><td>-1.928</td><td>-8.890</td><td>1.123</td><td>2.667</td><td>2.237</td><td>7.240</td></tr><tr><td>FEM<sup>a</sup></td><td>1.085</td><td>2.225</td><td>2.329</td><td>7.225</td><td>3.265</td><td>7.039</td><td>7.324</td><td>29.154</td><td>-0.556</td><td>-5.625</td><td>-2.496</td><td>-9.289</td><td>0.877</td><td>3.086</td><td>1.893</td><td>7.599</td></tr><tr><td>FVM<sup>a</sup></td><td>1.139</td><td>2.168</td><td>2.444</td><td>7.138</td><td>3.185</td><td>7.066</td><td>7.202</td><td>29.132</td><td>-0.711</td><td>-5.292</td><td>-2.859</td><td>-8.527</td><td>0.990</td><td>3.478</td><td>2.046</td><td>6.367</td></tr><tr><td>McC<sup>a</sup></td><td>1.085</td><td>2.220</td><td>2.328</td><td>7.223</td><td>3.267</td><td>7.037</td><td>7.086</td><td>29.307</td><td>-0.552</td><td>-5.578</td><td>-2.487</td><td>-9.036</td><td>0.872</td><td>2.674</td><td>1.886</td><td>7.321</td></tr><tr><td>STM<sup>a</sup></td><td>1.119</td><td>2.275</td><td>2.896</td><td>7.520</td><td>3.299</td><td>7.396</td><td>8.037</td><td>31.193</td><td>-0.593</td><td>-5.562</td><td>-4.478</td><td>-9.700</td><td>0.848</td><td>2.948</td><td>0.000</td><td>7.470</td></tr></tbody></table>

<table><caption>TABLE XV. Additional benchmark: a comparison of numerical parameters utilized to achieve the results plotted in Fig. 15 and the errors calculated in Table XIV.</caption><thead><tr><th>Method</th><th>Δt</th><th>Spatial order</th><th>Temporal order</th><th>Computation time (s)</th></tr></thead><tbody><tr><td>FC</td><td>1.00 × 10<sup>-3</sup></td><td>5th</td><td>4th</td><td>0.75<sup>b</sup> (MATLAB Classic)</td></tr><tr><td>DCG<sup>a</sup></td><td>5.00 × 10<sup>-4</sup></td><td>5th<sup>a</sup></td><td>2nd<sup>a</sup></td><td>0.97<sup>b</sup> (C++)</td></tr><tr><td>LCG<sup>a</sup></td><td>8.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>FEM<sup>a</sup></td><td>1.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>1st<sup>a</sup></td><td>Not available</td></tr><tr><td>FVM<sup>a</sup></td><td>3.31 × 10<sup>-3</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>McC<sup>a</sup></td><td>5.00 × 10<sup>-4</sup></td><td>2nd<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr><tr><td>STM<sup>a</sup></td><td>1.00 × 10<sup>-3</sup></td><td>1st<sup>a</sup></td><td>2nd<sup>a</sup></td><td>Not available</td></tr></tbody></table>

REFERENCES

¹D. S. Berger, J. K. Li, and A. Noordergraaf, "Differential effects of wave reflections and peripheral resistance on aortic blood pressure: A model-based study," *Am. J. Physiol.* **266**, H1626–H1642 (1994).

²K. Azer and C. S. Peskin, "A one-dimensional model of blood flow in arteries with friction and convection based on the Womersley velocity profile," *Cardiovasc. Eng.* **7**, 51–73 (2007).

³F. N. van de Vosse and N. Stergiopulos, "Pulse wave propagation in the arterial tree," *Annu. Rev. Fluid Mech.* **43**, 467–499 (2011).

⁴T. J. Hughes and J. Lubliner, "On the one-dimensional theory of blood flow in the larger vessels," *Math. Biosci.* **18**, 161–170 (1973).

⁵D. S. Berger, K. A. Robinson, and S. G. Shroff, "Wave propagation in coupled left ventricle-arterial system," *Hypertension* **27**, 1079–1089 (1996).

⁶N. Xiao, J. Alastruey, and C. Alberto Figueroa, "A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models," *Numer. Methods Biomed. Eng.* **30**, 204–231 (2014).

⁷J. Mynard, M. Davidson, D. Penny, and D. J. Smolich, "A simple, versatile valve model for use in lumped parameter and one-dimensional cardiovascular models," *Numer. Methods Biomed. Eng.* **28**, 626–641 (2012).

⁸J. P. Mynard and J. J. Smolich, "One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation," *Ann. Biomed. Eng.* **43**, 1443–1460 (2015).

⁹J. J. Wang and K. H. Parker, "Wave propagation in a model of the arterial circulation," *J. Biomech.* **37**, 457–470 (2004).

¹⁰M. S. Olufsen, C. S. Peskin, W. Y. Kim, E. M. Pedersen, A. Nadim, and J. Larsen, "Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions," *Ann. Biomed. Eng.* **28**, 1281–1299 (2000).

¹¹M. Saito, Y. Ikenaga, M. Matsukawa, Y. Watanabe, T. Asada, and P.-Y. Lagrée, "One-dimensional model for propagation of a pressure wave in a model of the human arterial network: Comparison of theoretical and experimental results," *J. Biomech. Eng.* **133**, 121005 (2011).

¹²T. Du, D. Hu, and D. Cai, "A fast algorithm for the simulation of arterial pulse waves," *J. Comput. Phys.* **314**, 450–464 (2016).

¹³L. O. Müller and E. F. Toro, "Well-balanced high-order solver for blood flow in networks of vessels with variable properties," *Numer. Methods Biomed. Eng.* **29**, 1388–1411 (2013).

¹⁴E. Boileau, P. Nithiarasu, P. J. Blanco, L. O. Müller, F. E. Fossan, L. R. Hellevik, W. P. Donders, W. Huberts, M. Willemet, and J. Alastruey, "A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling," *Numer. Methods Biomed. Eng.* **31**, e02732 (2015).

¹⁵S. Sherwin, V. Franke, J. Peiró, and K. Parker, "One-dimensional modelling of a vascular network in space-time variables," *J. Eng. Math.* **47**, 217–250 (2003).

¹⁶J. P. Mynard and P. Nithiarasu, "A 1D arterial blood flow model incorporating ventricular pressure, aortic valve and regional coronary flow using the locally conservative Galerkin (LCG) method," *Commun. Numer. Methods Eng.* **24**, 367–417 (2008).

¹⁷B. N. Steele, J. P. Ku, T. J. R. Hughes, and C. A. Taylor, "In vivo validation of a one-dimensional finite-element method for predicting blood flow in cardiovascular bypass grafts," *IEEE Trans. Biomed. Eng.* **50**, 649–656 (2003).

¹⁸K. S. Matthys, J. Alastruey, J. Peiró, A. W. Khir, P. Segers, P. R. Verdonck, K. H. Parker, and S. J. Sherwin, "Pulse wave propagation in a model human arterial network: Assessment of 1-D numerical simulations against in vitro measurements," *J. Biomech.* **40**, 3476–3486 (2007).

¹⁹L. Formaggia, D. Lamponi, and A. Quarteroni, "One-dimensional models for blood flow in arteries," *J. Eng. Math.* **47**, 251–276 (2003).

²⁰L. Formaggia, F. Nobile, A. Quarteroni, and A. Veneziani, "Multiscale modeling of the circulatory system: A preliminary analysis," *Comput. Visual. Sci.* **2**, 75–83 (1999).

²¹R. Raghu and C. Taylor, "Verification of a one-dimensional finite element method for modeling blood flow in the cardiovascular system incorporating a viscoelastic wall model," *Finite Elem. Anal. Des.* **47**, 586–592 (2011).

²²R. Raghu, I. Vignon-Clementel, C. A. Figueroa, and C. A. Taylor, "Comparative study of viscoelastic arterial wall models in nonlinear one-dimensional finite element simulations of blood flow," *J. Biomech. Eng.* **133**, 081003 (2011).

²³R. Mullen and T. Belytschko, "Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation," *Int. J. Numer. Methods Eng.* **18**, 11–29 (1982).

<footnote><sup>a</sup>Boileau *et al.*¹⁴</footnote>

<footnote><sup>b</sup>Without time shift (see Remark IV.1).</footnote>

<footnote><sup>c</sup>With time shift (see Remark IV.1).</footnote>

<footnote><sup>a</sup>Boileau *et al.*¹⁴</footnote>

<footnote><sup>b</sup>Remark IV.2.</footnote>
<footnote n="24">²⁴F. Amlani and N. M. Pahlevan, "A stable high-order FC-based methodology for hemodynamic wave propagation," *J. Comput. Phys.* **405**, 109130 (2020).</footnote>

<footnote n="25">²⁵F. Amlani, H. Wei, and N. M. Pahlevan, "A Fourier-based methodology without numerical diffusion for conducting dye simulations and particle residence time calculations," *J. Comput. Phys.* **493**, 112472 (2023).</footnote>

<footnote n="26">²⁶F. Amlani and O. P. Bruno, "An FC-based spectral solver for elastodynamic problems in general three-dimensional domains," *J. Comput. Phys.* **307**, 333–354 (2016).</footnote>

<footnote n="27">²⁷J. Alastruey, P. H. Charlton, V. Bikia, B. Paliakaitė, B. Hametner, R. M. Bruno, M. P. Mulder, S. Vennin, S. Piskin, A. W. Khir *et al.*, "Arterial pulse wave modelling and analysis for vascular age studies: A review from VascAgeNet," *Am. J. Physiol.* **325**, H1 (2023).</footnote>

<footnote n="28">²⁸Y. Shi, P. Lawford, and R. Hose, "Review of zero-D and 1-D models of blood flow in the cardiovascular system," *Biomed. Eng. Online* **10**, 33 (2011).</footnote>

<footnote n="29">²⁹X. Wang, J.-M. Fullana, and P.-Y. Lagrée, "Verification and comparison of four numerical schemes for a 1D viscoelastic blood flow model," *Comput. Methods Biomech. Biomed. Eng.* **18**, 1704–1725 (2015).</footnote>

<footnote n="30">³⁰O. P. Bruno and M. Lyon, "High-order unconditionally stable FC-AD solvers for general smooth domains. I. Basic elements," *J. Comput. Phys.* **229**, 2009–2033 (2010).</footnote>

<footnote n="31">³¹N. Albin and O. P. Bruno, "A spectral FC solver for the compressible Navier-Stokes equations in general domains. I: Explicit time-stepping," *J. Comput. Phys.* **230**, 6248–6270 (2011).</footnote>

<footnote n="32">³²M. Fontana, O. P. Bruno, P. D. Mininni, and P. Dmitruk, "Fourier continuation method for incompressible fluids with boundaries," *Comput. Phys. Commun.* **256**, 107482 (2020).</footnote>

<footnote n="33">³³K. Shahbazi, J. S. Hesthaven, and X. Zhu, "Multi-dimensional hybrid Fourier continuation–WENO solvers for conservation laws," *J. Comput. Phys.* **253**, 209–225 (2013).</footnote>

<footnote n="34">³⁴F. Amlani, H. S. Bhat, W. J. Simons, A. Schubnel, C. Vigny, A. J. Rosakis, J. Efendi, A. E. Elbanna, P. Dubernet, and H. Z. Abidin, "Supershear shock front contribution to the tsunami from the 2018 Mw 7.5 Palu, Indonesia earthquake," *Geophys. J. Int.* **230**, 2089–2097 (2022).</footnote>

<footnote n="35">³⁵J. P. Mynard, D. J. Penny, and J. J. Smolich, "Scalability and in vivo validation of a multiscale numerical model of the left coronary circulation," *Am. J. Physiol.* **306**, H517–H528 (2014).</footnote>

<footnote n="36">³⁶P. Reymond, Y. Bohraus, F. Perren, F. Lazeyras, and N. Stergiopulos, "Validation of a patient-specific one-dimensional model of the systemic arterial tree," *Am. J. Physiol.* **301**, H1173–H1182 (2011).</footnote>

<footnote n="37">³⁷S. Pagoulatou and N. Stergiopulos, "Evolution of aortic pressure during normal ageing: A model-based study," *PLoS One* **12**, e0182173 (2017).</footnote>

<footnote n="38">³⁸A. Aghilinejad, F. Amlani, S. P. Mazandarani, K. S. King, and N. M. Pahlevan, "Mechanistic insights on age-related changes in heart-aorta-brain hemodynamic coupling using a pulse wave model of the entire circulatory system," *Am. J. Physiol.* **325**, H1193–H1209 (2023).</footnote>

<footnote n="39">³⁹W. Jin and J. Alastruey, "Arterial pulse wave propagation across stenoses and aneurysms: Assessment of one-dimensional simulations against three-dimensional simulations and in vitro measurements," *J. R. Soc. Interface* **18**, 20200881 (2021).</footnote>

<footnote n="40">⁴⁰O. Vardoulis, E. Coppens, B. Martin, P. Reymond, P. Tozzi, and N. Stergiopulos, "Impact of aortic grafts on arterial pressure: A computational fluid dynamics study," *Eur. J. Vasc. Endovasc. Surg.* **42**, 704–710 (2011).</footnote>

<footnote n="41">⁴¹V. Fuster, R. A. O'Rourke, R. Walsh, and P. Poole-Wilson, *Hurst's the Heart* (McGraw Hill Professional, 2007).</footnote>

<footnote n="42">⁴²K. D. Lau and C. A. Figueroa, "Simulation of short-term pressure regulation during the tilt test in a coupled 3D–0D closed-loop model of the circulation," *Biomech. Model. Mechanobiol.* **14**, 915–929 (2015).</footnote>

<footnote n="43">⁴³H. Wei, F. Amlani, and N. M. Pahlevan, "Direct 0D–3D coupling of a lattice Boltzmann methodology for fluid-structure aortic flow simulations," *Int. J. Numer. Methods Biomed. Eng.* **39**, e3683 (2023).</footnote>

<footnote n="44">⁴⁴J. M. Huttunen, L. Kärkkäinen, and H. Lindholm, "Pulse transit time estimation of aortic pulse wave velocity and blood pressure using machine learning and simulated training data," *PLoS Comput. Biol.* **15**, e1007259 (2019).</footnote>

<footnote n="45">⁴⁵J. Carson, M. Lewis, D. Rassi, and R. Van Loon, "A data-driven model to study utero-ovarian blood flow physiology during pregnancy," *Biomech. Model. Mechanobiol.* **18**, 1155–1176 (2019).</footnote>

<footnote n="46">⁴⁶M. Willemet, S. Vennin, and J. Alastruey, "Computational assessment of hemodynamics-based diagnostic tools using a database of virtual subjects: Application to three case studies," *J. Biomech.* **49**, 3908–3914 (2016).</footnote>

<footnote n="47">⁴⁷C. D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. De Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot *et al.*, "Nektar++: An open-source spectral/hp element framework," *Comput. Phys. Commun.* **192**, 205–219 (2015).</footnote>

<footnote n="48">⁴⁸D. Moxey, C. D. Cantwell, Y. Bao, A. Cassinelli, G. Castiglioni, S. Chun, E. Juda, E. Kazemi, K. Lackhove, J. Marcon *et al.*, "Nektar++: Enhancing the capability and application of high-fidelity spectral/hp element methods," *Comput. Phys. Commun.* **249**, 107110 (2020).</footnote>

<footnote n="49">⁴⁹J. P. Mynard, *Computer Modelling and Wave Intensity Analysis of Perinatal Cardiovascular Function and Dysfunction* (University of Melbourne, 2011).</footnote>

<footnote n="50">⁵⁰N. Smith, A. Pullan, and P. J. Hunter, "An anatomically based model of transient coronary blood flow in the heart," *SIAM J. Appl. Math.* **62**, 990–1018 (2002).</footnote>

<footnote n="51">⁵¹T. W. Secomb, "Blood flow in the microcirculation," *Annu. Rev. Fluid Mech.* **49**, 443–461 (2017).</footnote>

<footnote n="52">⁵²S. Canić, "Blood flow through compliant vessels after endovascular repair: Wall deformations induced by the discontinuous wall properties," *Comput. Visual. Sci.* **4**, 147–155 (2002).</footnote>

<footnote n="53">⁵³J. Mynard, M. Davidson, D. Penny, and J. Smolich, "A numerical model of neonatal pulmonary atresia with intact ventricular septum and RV-dependent coronary flow," *Numer. Methods Biomed. Eng.* **26**, 843–861 (2010).</footnote>

<footnote n="54">⁵⁴W. W. Nichols, M. O'Rourke, E. R. Edelman, and C. Vlachopoulos, *McDonald's Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles* (CRC Press, 2022).</footnote>

<footnote n="55">⁵⁵J. Alastruey, K. H. Parker, and S. J. Sherwin, "Arterial pulse wave haemodynamics," in *11th International Conference on Pressure Surges* (Virtual PiE Led t/a BHR Group, 2012), pp. 401–443.</footnote>

<footnote n="56">⁵⁶K. W. Thompson, "Time-dependent boundary conditions for hyperbolic systems, ii," *J. Comput. Phys.* **89**, 439–461 (1990).</footnote>

<footnote n="57">⁵⁷S. J. Sherwin, L. Formaggia, J. Peiro, and V. Franke, "Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system," *Int. J. Numer. Methods Fluids* **43**, 673–700 (2003).</footnote>

<footnote n="58">⁵⁸O. P. Bruno, Y. Han, and M. M. Pohlman, "Accurate, high-order representation of complex three-dimensional surfaces via Fourier continuation analysis," *J. Comput. Phys.* **227**, 1094–1125 (2007).</footnote>

<footnote n="59">⁵⁹J. P. Boyd and J. R. Ong, "Exponentially-convergent strategies for defeating the Runge phenomenon for the approximation of non-periodic functions, part I: Single-interval schemes," *Commun. Comput. Phys.* **5**, 484–497 (2009).</footnote>

<footnote n="60">⁶⁰L. N. Trefethen and D. Bau, *Numerical Linear Algebra* (SIAM, 2022).</footnote>

<footnote n="61">⁶¹H. Musa, I. Saidu, and M. Y. Waziri, "A simplified derivation and analysis of fourth order Runge-Kutta method," *Int. J. Computer Appl.* **9**, 51–55 (2010).</footnote>

<footnote n="62">⁶²S. Abarbanel, D. Gottlieb, and M. H. Carpenter, "On the removal of boundary errors caused by Runge-Kutta integration of nonlinear partial differential equations," *SIAM J. Sci. Comput.* **17**, 777–782 (1996).</footnote>

<footnote n="63">⁶³M. H. Carpenter, D. Gottlieb, S. Abarbanel, and W. S. Don, "The theoretical accuracy of Runge-Kutta time discretizations for the initial boundary value problem: A study of the boundary error," *SIAM J. Sci. Comput.* **16**, 1241–1252 (1995).</footnote>

<footnote n="64">⁶⁴T. Elling, "GPU-accelerated Fourier-continuation solvers and physically exact computational boundary conditions for wave scattering problems," Ph.D. thesis (California Institute of Technology, 2012).</footnote>

<footnote n="65">⁶⁵P. J. Roache, "Code verification by the method of manufactured solutions," *J. Fluids Eng.* **124**, 4–10 (2002).</footnote>

<footnote n="66">⁶⁶C. Roy, C. Ober, and T. Smith, "Verification of a compressible CFD code using the method of manufactured solutions," in *32nd AIAA Fluid Dynamics Conference and Exhibit* (AIAA, 2002).</footnote>

<footnote n="67">⁶⁷J. M. Vedovoto, A. D. Silveira Neto, A. Mura, and L. F. Figueira da Silva, "Application of the method of manufactured solutions to the verification of a pressure-based finite-volume numerical scheme," *Computers Fluids* **51**, 85–99 (2011).</footnote>

<footnote n="68">⁶⁸S. Sankaranarayanan, N. Shankar, and H. Cheong, "Three-dimensional finite difference model for transport of conservative pollutants," *Ocean Eng.* **25**, 425–442 (1998).</footnote>
<footnote n="69"><sup>69</sup>M. B. Kalinowska and P. Rowiński, “Truncation errors of selected finite differ-
ence methods for two-dimensional advection-diffusion equation with mixed
derivatives,” *Acta Geophys.* **55**, 104–118 (2007).</footnote>

<footnote n="70">⁷⁰I. M. Babuska and S. A. Sauter, “Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?,” *SIAM Rev.* **42**, 451–484 (2000).</footnote>

<footnote n="71">⁷¹M. Lyon and O. P. Bruno, “High-order unconditionally stable FC-AD solvers for general smooth domains. II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations,” *J. Comput. Phys.* **229**, 3358–3381 (2010).</footnote>

<footnote n="72">⁷²N. Albin, O. P. Bruno, T. Y. Cheung, and R. O. Cleveland, “Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams,” *J. Acoust. Soc. Am.* **132**, 2371–2387 (2012).</footnote>

<footnote n="73">⁷³MathWorks, MATLAB Version 9.14.0.2674353 (R2023a) Update 7 (The Mathworks, Inc., Natick, MA, 2023).</footnote>

<footnote n="74">⁷⁴S. Sharkey, *A Guide to Interpretation of Hemodynamic Data in the Coronary Care Unit* (Lippincott-Raven Publishers, 1997).</footnote>

<footnote n="75">⁷⁵M. J. Kern, M. J. Lim, and J. A. Goldstein, *Hemodynamic Rounds: Interpretation of Cardiac Pathophysiology from Pressure Waveform Analysis* (John Wiley & Sons, 2018).</footnote>

<footnote n="76">⁷⁶C.-H. Chen, B. Fetics, E. Nevo, C. E. Rochitte, K.-R. Chiou, P.-A. Ding, M. Kawaguchi, and D. A. Kass, “Noninvasive single-beat determination of left ventricular end-systolic elastance in humans,” *J. Am. Coll. Cardiol.* **38**, 2028–2034 (2001).</footnote>

<footnote n="77">⁷⁷R. Klabunde, *Cardiovascular Physiology Concepts* (Lippincott Williams & Wilkins, 2011).</footnote>