Synthetic Dataset Generation, Physical Consistency and Active Learning for Shape-Informed Surrogate Modeling of Hemodynamics

Introduction

Computational fluid dynamics (CFD), applied to blood flow modeling and simulation, can support non-invasive diagnosis and assessment of cardiovascular diseases, as well as treatment planning and patient monitoring (see, e.g., [1, 2]). The development of numerical solvers, together with the constantly growing potential of computing resources, allows tackling applications of extremely high complexity. However, especially in the case of aortic blood flow, where fluid dynamics is characterized by highly nonlinear regimes, strongly dependent on patient geometry, the physical and geometrical complexity still requires considerable computational resources for physiological patient-specific simulations. This aspect is particularly relevant for outer-loop applications concerning, e.g., inverse problems, optimization, or uncertainty quantification, in which a large number of (forward) simulations is required at each iteration step. Reduced-order models and scientific machine learning (ML) can address this challenge by considering approximate dynamics on lower-dimensional spaces. These techniques rely on existing simulation (training) data to define high-fidelity reduced representations of target applications, and their capabilities are strongly dependent on data availability.

Compared with other applications of model order reduction, biomedical applications are still constrained by limited data availability and substantial inter-patient variability in both anatomy and physiological conditions, concerning anatomical shape, cardiac output, peripheral pressure, and mechanical parameters, among others. Since full-order CFD simulations are computationally expensive, the predictive accuracy of inter-patient surrogate models is essential for clinically relevant applications.

This work addresses different challenges concerning inter-patient variability in the context of surrogate models of aortic blood flow. We study the approximation properties and the performance of ML surrogate models trained on synthetic data, i.e., data generated through high-fidelity numerical blood flow simulations (full-order models) on a dataset of variable patient anatomies sharing a common reference shape. In particular, we investigate the accuracy in terms of relevant biomarkers, namely, peak velocity, pressure, and wall shear stress (WSS), when predicting these quantities across shapes.

Handling the geometrical variability in the context of scientific ML requires tailored strategies for developing reduced spaces and deep learning architectures that can scale better during training and also generalize accurately to unseen geometries. Previous studies explored the use of the reduced basis method combined with domain decomposition strategies to attach archetypal subdomains [3], and combined with registration techniques to transport a reduced basis computed on a template shape onto patient geometries [4]. Recent works focusing on ML models have addressed the development of graph neural networks (GNNs) approximating the averaged quantities of interest (QoIs) from measurement data on the vessel centerline [5], the learning of equivalent one-dimensional (1D) models from measurement data, to be applied to variable three-dimensional (3D) anatomies [6], common GNN architectures applied to the 3D mesh [7], shape-informed GNNs, combined with registration strategies and trained on registered simulation data [8], and GNNs trained on image data of carotid arteries, which leverage first principles as soft constraints and equivariant layers [9]. In this work, we focus on graph- and transformer-based architectures [10], as they have proven effective for large-scale applications, and apply them to a database of physiological aortic geometries [11] generated from parametric deformations.

The accuracy of surrogate models and their capability of generalizing to different patients strongly depend on the variability of the available training data. This aspect is particularly important in the context of blood flow in large vessels, where geometrical details and patient-specific boundary conditions strongly influence the fluid dynamics. Flexible registration techniques and parameter-dependent shape databases, such as the one considered in this work, provide useful tools to increase the size of the training set with additional simulations. Typical methods to perform data augmentation, increasing the number of different anatomical shapes, include steady-state free precession (SSFP) magnetic resonance imaging (MRI) to directly acquire them from patients, statistical shape models [5], parametric deformations [11], and generative models [12]. Although the availability of large datasets keeps increasing, the cost of generating high-fidelity CFD solutions remains high. The exploration of the space of all possible anatomical shapes must thus be optimized to reduce the computational costs. In this regard, active learning (AL) [13] and reinforcement learning provide promising strategies to select the most informative samples to be added to the training dataset. We test different strategies for the selection of additional training samples, comparing label-dependent (supervised) techniques, i.e., techniques requiring the full CFD solution computed on the candidate geometries, with label-free (unsupervised) ones, which only rely on surrogate approximations and are therefore suitable for practical applications.

A further aspect investigated in this work concerns the introduction of physically consistent soft constraints in the training loss, besides purely data-dependent residuals computed from the solutions on the different shapes. The enforcement of these constraints has been tested extensively in the literature. Among others, physics-informed neural networks [14] and neural emulators [15] use ML models as ansatz spaces while minimizing equation-based loss terms or residuals of numerical schemes via automatic differentiation. We are interested instead in physics constraints used as regularizers for surrogate models, with the aim of modeling a whole joint distribution of input parameters and output solutions, where the parameters represent complex boundary conditions or different computational domains. The efficiency of such regularizers with respect to state-of-the-art ML methods is not often benchmarked, making it difficult to discern whether purely data-informed or physics-informed losses yield the highest gain in generalization error relative to computational cost.

This work provides the following main contributions.

(i) We present a detailed comparison of two state-of-the-art GNN and transformer-based architectures in terms of accuracy and computational cost. In particular, we show that the partial differential equation (PDE) residuals of the predicted solutions can be used as a posteriori error estimators.

(ii) Concerning AL strategies, our results show that physics-based strategies, often considered the most effective (e.g., [13]), are not necessarily the best choices for the considered applications.

(iii) Physical consistency soft constraints are more effective as a fine-tuning objective than as a regularizer applied throughout training.

(iv) We analyze how the considered surrogate models perform for out-of-distribution samples, based on shape databases published in different contexts, discussing how to improve the generalization capabilities of the surrogate models.

The rest of this paper is organized as follows: Section 2 summarizes the main aspects of the dataset generation and curation including the underlying modeling assumptions and their implications for the designed surrogate models. Section 3 then focuses on the validation of state-of-the-art deep learning architectures for surrogate models of aortic flows (Graph U-Nets and Transformer U-Nets), including the analysis of inter-patient AL methods, the application of physics-based soft constraints during training and in fine-tuning, and the evaluation of the models’ generalization capabilities. Section 4 discusses the conclusions, limitations, and possible future directions.

Dataset generation

This section describes the dataset generation, discussing the shape variability and parametrization of physiological, healthy aortas (Section 2.1), the relevant aspects of the vascular flow modeling (Section 2.2), and the details of the CFD solver (Section 2.3).

Geometric modeling

The shape database used in this work considers a geometric parametrization of healthy aortas based on the SynthAorta dataset [11], generated by mapping a single common reference centerline model. Equipped with a parametrized transformation and radii at each point of the centerline, a broad range of healthy physiological aortas can be captured, albeit in a simplified manner. Bośnjak et al. [11] proposed a set of hyperparameters defining healthy aortic geometries, which have a clear physiological meaning, are directly measurable, and are mutually independent. The stochastically sampled values for these hyperparameters are sourced from patient-specific data whenever possible. For the sources of the geometric hyperparameters, their definitions, and related distributions, we refer the reader to [11].

Based on the available collection of centerlines with corresponding radial data, smooth aortic surfaces of the SynthAorta dataset are reconstructed via convolution surfaces [18, 19]. Structured meshing methods based on convolution surfaces have been employed for tube-like structures [20], such as healthy aortas, but have also been extended to describe more general biological shapes [21], pathological vessel configurations such as aortic dissection [22], or vessel wall geometries [23].

Since each aorta in the SynthAorta parametrization is represented through a limited set of hyperparameters, it inevitably cannot reproduce the full geometric detail of any individual patient’s anatomy: local, high-fidelity features not captured by the chosen hyperparameters are necessarily lost. For aortic flow, however, Bošnjak et al. [24] have shown that such local geometric variations play only a secondary role for global hemodynamic QoIs, such as volumetric flow rates, pressures, and averaged WSS. Whenever local QoIs must be resolved with high accuracy, e.g., at elevated Reynolds numbers, patient-specific CFD remains the gold standard; global QoIs, in contrast, can be inferred reliably given a sufficiently rich training dataset. The methods presented in this work therefore necessarily trade some geometric fidelity to individual patients for the ability to generalize across the broader anatomical variability captured by the parametrization.

For the purposes herein, we consider cases 1 to 1000 of the SynthAorta dataset [11], wherein the hyperparameters are sampled with a Latin hypercube method, covering the entire parameter space. To showcase the geometric variability within the training and test data, Figure 1 shows eight samples.

Vascular flow modeling

Blood flow

To generate synthetic training data, we model the blood flow as an incompressible fluid. Namely, blood velocity and (kinematic) pressure are assumed to be described by the Navier–Stokes equations in a domain $\Omega\subset\mathbb{R}^{3}$, representing any aortic geometry in the dataset, and in a selected time interval $[0,T]$:

where $\nu$ denotes the kinematic viscosity, and the symmetric gradient is defined as $\nabla^{S}(\cdot) := \frac{1}{2}\left[\nabla(\cdot) + \nabla(\cdot)^{\mathsf{T}}\right]$. To account for the blood’s rheological behavior, a generalized Newtonian model is considered, where the apparent viscosity depends on the local shear rate, $\dot{\gamma}(x) := \sqrt{2\nabla^{S}u : \nabla^{S}u}$, via the Carreau model, i.e.,

with $\eta_{0}$, $\eta_{\infty}$, $\lambda$, and $n$ being rheological parameters. The physiological fluid parameters for Equation (3) are taken from [25]. In addition to the primary unknowns, i.e., the fluid velocity $u$ and pressure $p$, and their averages on the inlet and outlets, the WSS vector $\tau_{\mathrm{wall}} \in\mathbb{R}^{3}$, defined as

where $n$ is the unit outward normal vector and $\rho$ the fluid density, plays a critical role in various pathologies such as abdominal aortic aneurysms [16] or aortic dissection [17]. Defined on $\Gamma_{\mathrm{wall}}$ only, it is obtained in a postprocessing step. For the sake of brevity, other biomarkers relevant in hemodynamics are not discussed, noting that such QoIs depend on the fluid velocity and pressure in some way and might hence be inferred from the predicted primary unknowns.

Boundary conditions

The system of Equations (1)–(2) needs to be completed by suitable boundary conditions, which play a crucial role in recovering realistic and physiological flow fields. These boundary conditions depend on the considered patient geometry and on a set of parameters which, typically, need to be tuned according to the specific case. To handle this aspect, we consider an automatic parameter tuning procedure for each new geometry, with the goal of adapting the boundary conditions to maintain similar flow conditions across the dataset, corresponding to the physiological regime (see, e.g., Mills et al. [26]). In what follows, we consider the domain boundary to be decomposed as

i.e., inlet, vessel wall, and (four) outlet boundaries, respectively (see Figure 2, left). The definition of boundary conditions is then based on the following assumptions:

(i) We assume the heart rate and inflow profile’s temporal scaling to be equal for all considered cases, defining the spatial scaling according to the inlet geometry.

(ii) We neglect wall deformation, prescribing no-slip conditions at the vessel wall ($\Gamma_{\mathrm{wall}}$) for all cases.

(iii) The flow split is dictated by the outlet surface area, similar to [8, 27].

Physiological flow conditions are enforced by controlling target flow rates $Q_i$ and mean outlet pressures $\bar{p}_i$ on $\Gamma_i$ defined as

Inflow boundary condition At the inlet surface $\Gamma_0$, i.e., the cross section at the aortic root (see Figure 2), we prescribe a Dirichlet condition of the form

where $n$ stands for the unit outward normal vector at the inlet surface, $\eta_x(x)$ is a quasi-parabolic spatial profile, constructed by adopting the technique proposed in [28], and $\eta_t(t)$ is a time-dependent profile obtained from clinical measurements from a specific patient (see Figure 2, right). Finally, the scaling factor $\eta_Q$ is introduced to account for volumetric flow variation across geometries. Westerhof and Elzinga [29] found stroke volume,

per unit body mass to be independent of organism size. Assuming geometric similarity, organism mass is proportional to a volume scale, $\tilde{r}_0^3$, which implies $SV \propto\tilde{r}_0^3$. Consequently, preserving $Q_0/SV$, i.e., the heart rate $1/T_p$, requires $Q_0 \propto\tilde{r}_0^3$, i.e., a constant $Q_0/\tilde{r}_0^3$, in accordance with Murray’s cubic law [30].

Denoting quantities defined on the reference case analyzed in [31] with $(\widetilde{\cdot})$, we compute the inlet mean radius of the reference case from $\pi\tilde{r}_0^2 = \tilde{A}_0$ given the inlet area $\tilde{A}_0$. That is, we scale the inlet velocity by a constant factor $\eta_Q$ per geometry according to

Using this scaling, we obtain the inflow rates shown in Figure 3 (left). For demonstration purposes, Figure 3 (right) also contains a consistency check: inserting the numerically obtained flow rate into $\widetilde{Q}_0=\tilde{r}_0^3/r_0^3\,Q$ yields the reference flow rate from Bäumler et al. [31] itself.

Windkessel parameter tuning At the outlet, boundary conditions are imposed through lumped-parameter models, i.e., relating the mean outlet pressure to the outgoing flow. A popular choice consists of a 3-element Windkessel model on each outlet, in which the compliant pressure $p_{c,i}$ on the outlet $\Gamma_i$, to be enforced as the mean pressure $\bar{p}_i$ on that outlet, is computed as the solution of an ordinary differential equation of the form

depending on the proximal and distal resistances, $R_{p,i}$ and $R_{d,i}$, accounting for the flow resistance of proximal and distal vasculature, respectively, on a capacitance $C_i$, accounting for vessel elasticity, and on the distal (also venous or basal) pressure $p_{d,i}$. In Equations (8)–(9), the proximal pressure $p_{p,i}$ is an intermediate unknown, modeling the pressure at the compliant node between the proximal and distal resistances.

In this work, we consider the scheme by Grinberg and Karniadakis [27], neglecting the proximal resistances $R_{p,i}$ and yielding the 2-element (RC) model

The parameters of the Windkessel model are tuned in order to prescribe a systemic vascular pressure and a target flow split among the four outlets, according to the relation $Q_i/Q_j \approx R_{d,j}/R_{d,i}$, $(i,j=1,\ldots,4)$, which relies on a simplified zero-dimensional representation of the aortic domain [27]. We consider the flow split used in Katz et al. [32], i.e., 16.81% brachiocephalic trunk (BCT), 8.6% left common carotid artery (LCC), 8.21% left subclavian artery (LSA), and 66.38% thoracic aorta (TA). In particular, the resistance parameters are tuned for the first patient in the SynthAorta dataset (case 1), then scaled to other geometries in a second step. The distal pressures are set to $p_{d,i}=75$ mmHg. For case 1, we also fix the capacitance $C_4=12\times10^{-10}\ \mathrm{m^3/Pa}$ at the TA outlet. This procedure leads to the value $\widetilde{R}_{d,4}=2.5\times10^7\ \mathrm{Pa\,s/m^3}$ for the TA outlet, achieving physiological target pressure levels under resting conditions of $\approx130$ mmHg in systole and $\approx80$ mmHg in diastole, see, e.g., Mills et al. [26]. The resulting flow rates and pressures for case 1 are shown in Figure 4.

Since the inlet velocity is scaled across the database according to Equation (7), the proximal and distal resistances and the capacitance at each outlet are scaled accordingly to achieve similar pressure levels for the other cases. Moreover, to preserve the convergence properties, the capacitances $C_i$, $i=1,2,3,4$, are scaled to keep the product $R_{d,i}C_i$ constant on each outlet.

Based on the inflow scaling in Equation (7) and assuming that the flow split remains unchanged across the geometries, we can deduce scaling parameters by matching coefficients in

to obtain $p_{c,i}\equiv\widetilde{p}_{c,i}$, which for $Q=\frac{r_0^3}{\widetilde{r}_0^3}\widetilde{Q}$ from Equation (7) gives

Note that the resistance scaling is consistent with Westerhof and Elzinga [29]: since volume and mass scale similarly with organism size (assuming some constant average density of the organism), the chosen volume ratio in Equation (10) leads to a linear relationship between the resistance and organism size.

The resulting flow split and mean outlet pressure curves are given in Figure 5 for six selected cases from the SynthAorta dataset, showing that the chosen scaling achieves the desired flow distribution among outlets and that the pressure level is kept roughly constant, despite geometric variations. Setting a constant distal pressure $p_{d,i} = 75$ mmHg for all cases is vital for fast convergence to the periodic state and for achieving the desired pressure levels, as discussed in § 2.3.1.

Note, however, that the capacitance is scaled with the inverse factor of $R_{d,i}$, which can trigger numerical instabilities when the distal resistance is increased and the capacitance is decreased according to the proposed tuning procedure. This might lead to some oscillations for the cases with the smallest radii.

Computational fluid dynamics

The Navier–Stokes equations for incompressible flow of generalized Newtonian fluids, see Equations (1)–(2), are solved using an $L^2$-conforming discontinuous Galerkin approach based on a time-splitting scheme [33], in which the solution is split into a convective step with a fully explicit extrapolation of the nonlinear term, a pressure Poisson step, and a step projecting the velocity onto a divergence-free space. The formulation enforces weak mass conservation and energy stability by additional stabilization terms, which are critical for robust and accurate discretizations of turbulent flows on non-resolved scales. This strategy avoids solving linearized or nonlinear coupled velocity–pressure systems [34], typically leading to higher throughput in vascular flow. We formally achieve accuracy of order 2 in time, which is lower than that of consistent splitting schemes [35, 36], but sufficient given the spatial discretization accuracy within the SynthAorta dataset and practical performance considerations rooted in the problem parameters.

The solver is implemented in ExaDG [37], using the matrix-free framework provided by deal.II [38]. The preconditioning technique for the solution of the linear systems for the pressure Poisson equation and the momentum equation is based on a mixed-precision, hybrid multigrid method with matrix-free Chebyshev-accelerated Jacobi smoothers [39, 40]. Moreover, backflow stabilization via a velocity penalization method based on [41] is used on the Neumann boundaries to prevent instabilities for re-entrant flow (see also [42, 43]).

Convergence studies and grid selection

The numerical approach has to be verified in terms of spatiotemporal accuracy, while ensuring that the periodic state is reached up to the desired tolerance. The grids employed for CFD analysis, which feature trilinear mappings provided in the SynthAorta dataset, are depicted in Figure 1. We employ an adaptive backward differentiation formula of order 2 for time stepping based on the Courant–Friedrichs–Lewy (CFL) number, such that smaller grids automatically translate to smaller time steps. Second-order discontinuous polynomials for the velocity approximation and discontinuous linear polynomials for the pressure approximation are used. A grid sensitivity study on the whole set of geometries is not practically feasible. Hence, to validate the solver, we perform convergence studies for selected cases ($i = 1, 1000, 2000, 3000, 4000, 5000$) of the SynthAorta dataset, which are considered to be representative of the main features of the dataset. The geometric representation is trilinear, and the points inserted during refinement conform to the underlying geometry, see [11], such that the geometric approximation improves under grid refinement.

The observed quantities are the mean pressures and flow rates on the outlets, $\bar{p}_i$ and $Q_i$, $i = 1,2,3,4$, respectively, and the averaged WSS,

where $\tau_{\mathrm{wall}}$ is the WSS defined in Equation (4). For these quantities, we define the time-averaged error

for $\phi= Q,\bar{p},\overline{\tau}_{\mathrm{wall}}$, with respect to a reference solution computed on a fine grid containing $\approx10.21 \times10^{6}$ degrees of freedom (DoFs) and approximately $8 \times10^{4}$ time steps (uniform refinement level $l = 3$). We show the results of the convergence study for case 1, which yielded the largest error in the WSS, the most critical quantity. In particular, relative errors under uniform grid refinement are reported in Table 1. The mean outlet pressure shows relatively low errors already for grid level $l = 0$, while the flow rates are more prone to approximation errors. The relative difference in the integral mean WSS is the most sensitive QoI, with errors in the range of 10% even for a grid uniformly refined twice, and errors up to 60% in most geometries on levels 0 and 1. Based on the velocity and pressure errors, we select level $l = 1$, noting, however, that $l = 2$ would be required for reasonably accurate WSS. The corresponding transient profiles of the monitored quantities are depicted in Figure 6.

$l$	$\epsilon_{\bar{p}_1}$	$\epsilon_{\bar{p}_2}$	$\epsilon_{\bar{p}_3}$	$\epsilon_{\bar{p}_4}$	$\epsilon_{Q_1}$	$\epsilon_{Q_2}$	$\epsilon_{Q_3}$	$\epsilon_{Q_4}$	$\epsilon_{\overline{\tau}_{\mathrm{wall}}}$
0	0.0112	0.0056	0.0089	0.0082	0.4595	0.4823	0.4659	0.4682	0.2235
1	0.0019	0.0022	0.0019	0.0017	0.0572	0.0657	0.0641	0.0479	0.2716
2	0.0012	0.0014	0.0013	0.0007	0.0214	0.0266	0.0301	0.0086	0.1076

Periodicity In addition to the overall accuracy, we investigate the relative difference between the solutions over different cardiac cycles, i.e., within time intervals $[0.1, 0.78] + 0.78\,n$, $n = 0,1,2,3$. Results are presented only for the reference case 1 of the dataset. Namely, to investigate the convergence to the periodic state, we consider, for each refinement level $l$, as the reference solution the one computed on the same level $l$ but for the fourth cardiac cycle. Table 2 shows negligible errors with respect to periodicity compared to the spatiotemporal errors shown in Table 1, already for the first cycle (for $t \geq0.1$ s). Overall, we conclude that the first cardiac cycle is acceptable for all grid levels, since the solutions have mostly reached the periodic state, which is attributable to the specific Windkessel parameters chosen.

Cardiac cycle	$\epsilon_{\overline{p}_4}$ in %				$\epsilon_{Q_4}$ in %				$\epsilon_{\overline{\tau}_{\mathrm{wall}}}$ in %
	$l = 0$	$l = 1$	$l = 2$	$l = 3$	$l = 0$	$l = 1$	$l = 2$	$l = 3$	$l = 0$	$l = 1$	$l = 2$	$l = 3$
1	0.153	0.136	0.135	0.129	0.231	0.236	0.151	0.179	0.720	1.129	1.381	1.087
2	0.032	0.026	0.028	0.018	0.193	0.266	0.221	0.182	0.385	0.703	0.933	0.961

Uncertainty of QoIs

The inherent complexity of the vascular flow data generated from the SynthAorta dataset can be illustrated through the resulting uncertainty in the considered QoIs. Figure 7 depicts violin plots for cases 1 to 1000 over the cardiac cycle, demonstrating (i) the effectiveness of the algorithm used to tune the Windkessel parameters across geometries, leading to very little fluctuation from the target pressure levels despite vast geometric variation, and (ii) the impact of the geometric variation, explicitly through changes of the cross section and implicitly through the scaling of volumetric flow rate $Q_0$ via $Q_0 = \frac{r_0^3}{\widetilde{r}_0^3}\widetilde{Q}$, see Equation (7).

The variations in the WSS are directly related to the variations in the velocity: the violin plots of the flow rate show large variations in systole, while during the diastolic phase, variations are limited. The WSS shows similar variations over the cardiac cycle, but the widths of the violins are distinctly different. The flow rates show top-heavy, more uniform distributions, while the WSS features bottom-heavy distributions with thin but non-vanishing top regions. In diastole, the variation in the WSS seems more pronounced, although the mean absolute value is also smaller. In summary, peak systole shows the largest variations, while the correlation under geometric variation between velocities and flow rates, on the one hand, and WSS, on the other hand, is non-trivial.

Since the observed absolute values and their variations are largest at peak systole, we choose the point of maximum inflow at $t = 0.19$ s as the design point for our surrogates.

Dataset of full-order solutions

Maximizing the number and/or the fidelity of the CFD solutions considered for training and testing the ML models, while keeping the computational resources required to generate them tractable, requires handling the increasing amount of error within the data. In the considered cases, we fix the spatiotemporal resolution, as described in § 2.3.1, and we assume that the solution reaches a periodic state from 0.1 s on. In this setting, the solver tolerances and related numerical settings are chosen to obtain results as fast as possible while preserving overall acceptable accuracy and stability. This drive for efficiency leads to decreased robustness of the linear solvers. Hence, given a certain computational budget, a trade-off between accuracy/stability and number of samples has to be found.

Herein, we choose the route of postprocessing the simulation dataset to detect parameter configurations that deliver unphysical or diverging solutions, by (i) evaluating whether the final simulation time of 0.78 s was reached, and (ii) introducing a tolerance on the second derivative of the flow rate $\ddot{Q}_i$ and on the mean outlet pressure $\bar{p}_i$. A total of 15 of the first 1000 cases fail within the first time step due to bad element quality, while 2 cases are flagged as unphysical due to the variation in the pressure and flow rate, see Figure 8. These cases can be successfully simulated using lower grid resolution and/or different CFD solver settings (results omitted), indicating that the geometric parameters and the Windkessel tuning themselves are not problematic, but the CFD solver settings combined with the element quality are troublesome. In order not to introduce bias into the data, we restrict ourselves to these filtering criteria, although they are rather loose, so that more than 98% of the first 1000 cases remain in the dataset.

Additional filtering is performed during the processing of training and test data subsequently used in ML models, to remove cases which might not represent relevant physiological regimes. Specifically, we remove geometries in which the corresponding pressure solution yielded a maximum pressure drop at peak systole higher than a prescribed threshold of 5 mmHg. This threshold value is motivated by the distributions of filtered average and maximum values of the velocity, pressure drop, displacement, and WSS fields for the first 5000 cases of the SynthAorta dataset.

Adopting the filtering techniques described above, the resulting total number of accepted geometries is 4660 out of 5000 (93.2%). Out of these 4660 cases, $n_{\mathrm{train}} = 4560$ are used for training and $n_{\mathrm{test}} = 100$ for testing. The resulting distributions of the average and maximum values of the velocity, pressure drop, displacement, and WSS fields of the final training and test datasets are shown in Figure 9.

To reduce the dimension of the reduced space, the discontinuous second-order velocity $u$, the discontinuous first-order pressure $p$, and the continuous second-order nodal displacement field $d$ are interpolated onto continuous finite element spaces on uniformly refined grids. This step accounts for the higher-order contributions, since standard Graph U-Nets have no notion of higher-order elements. This leads to a scalar finite element space with 16125 nodes, so $N_u = N_d = 48375$ DoFs for vector fields (velocity and displacement), and $N_p = 16125$ DoFs for the pressure field. The WSS $\tau_{\mathrm{wall}}$ is obtained in postprocessing.

Reconstruction error on reduced-dimensional spaces

An effective way to explore the complexity of the snapshot dataset is to compute reduced-dimensional spaces via randomized singular value decomposition (rSVD) and investigate how the relative $\ell^2$-reconstruction error of the velocity and pressure drop decays with respect to its dimension. Knowing how the reconstruction error decays with respect to the reduced dimensions is crucial to select the most effective surrogate modeling approach. If the decay is fast or the threshold on the accuracy is reached with only a few rSVD modes, a linear subspace surrogate model might be sufficient (e.g., the reduced basis method or proper orthogonal decomposition combined with neural networks, POD-NN). Conversely, if the decay is slow, nonlinear surrogates and manifold learning approaches might be necessary and more effective. Given a set of training velocity and pressure snapshots $\{u_i\}_{i=1}^{n_{\mathrm{train}}}$ and $\{p_i\}_{i=1}^{n_{\mathrm{train}}}$, with $u_i \in\mathbb{R}^{N_u}$ and $p_i \in\mathbb{R}^{N_p}$, respectively, we define the snapshot matrices gathering the mean-centered snapshots as the column vectors:

where $u_{\mathrm{mean}}$ and $p_{\mathrm{mean}}$ denote the means computed over the training snapshots. We then compute their rank-$r$ rSVD, i.e.,

where $\Sigma_r = \operatorname{diag}(\sigma_1,\ldots,\sigma_r)$, $\Lambda_r = \operatorname{diag}(\lambda_1,\ldots,\lambda_r)$ contain the $r$ largest singular values in decreasing order, the $r$ columns of $\Phi_r = [\phi_1|\cdots|\phi_r]$ and $\Psi_r = [\psi_1|\cdots|\psi_r]$ are the leading left singular vectors (spatial modes), and the columns of $V_r = [v_1|\cdots|v_r]$ and $Q_r = [q_1|\cdots|q_r]$ are the corresponding right singular vectors. The columns of $\Phi_r$ and $\Psi_r$ define the bases of the (reduced) $r$-dimensional spaces that approximate the original set of snapshots. For both training and test snapshots, the reconstruction is defined by the orthogonal projection onto these spaces, using the same training means and bases, i.e.,

For the training snapshots, $i=1,\ldots,n_{\mathrm{train}}$, this coincides with the columns of $U_r$ and $P_r$ shifted by the means, i.e., $u_{r,i}=u_{\mathrm{mean}}+\sum_{j=1}^{r}\sigma_j(v_j)_i\phi_j$ and $p_{r,i}=p_{\mathrm{mean}}+\sum_{j=1}^{r}\lambda_j(q_j)_i\psi_j$. The relative reconstruction error, for an arbitrary training or test index $i$, is then computed as the $\ell^2$-norm of the difference between the original and reconstructed snapshot divided by the $\ell^2$-norm of the original centered snapshot, i.e.,

In the following study, we will consider the mean and maximum $\ell^2$-reconstruction errors on both the training and the test snapshot sets.

Figures 10 and 11 show the resulting relative reconstruction errors for increasing dimension $r$ of the reduced space for velocity and pressure, respectively. As expected, increasing $r$ yields a decrease in the reconstruction errors. Moreover, increasing $r$ from 500 to 1000 still yields a significant improvement, indicating the high complexity of the snapshot dataset. As the number of training samples increases, the test reconstruction error is observed to practically stabilize during the transition from 2048 to 4096 samples, albeit without reaching the training reconstruction error values. The training reconstruction error grows as the number of training samples increases to 4096. This behavior also indicates that the considered dataset is highly complex and cannot be easily represented by a low-dimensional space. Comparing the velocity and pressure reconstruction errors reveals that the velocity error is in general higher than the pressure error, indicating that pressure fields are more easily represented using a reduced basis.

Remark 2.2 On the training set, the decay of the reconstruction errors with increasing $r$ is closely related to the decay of the singular values $\sigma_j$ and $\lambda_j$: faster singular-value decay indicates that fewer modes are needed to capture a given fraction of the centered snapshot energy. However, the singular values alone do not determine the mean and maximum of the per-snapshot relative errors considered here.

Evaluating reconstruction errors on a held-out test set additionally assesses how well the reduced space computed from the training snapshots represents unseen snapshots drawn from the same underlying distribution. Applying this procedure to other datasets, e.g., for aortas affected by coarctation [8] or for simplified aortic arches [44], permits an empirical comparison of their linear approximation complexity, provided that sample sizes, sampling procedures, field representations, centering, and error norms and normalizations are comparable.

The relative $\ell^2$-reconstruction error can also be used to quantify the complexity of the parameter space for geometry generation. To this end, we consider the set of displacement fields $\{d_i\}_{i=1}^{n_{\mathrm{train}}}$ and $\{d_i\}_{i=1}^{n_{\mathrm{test}}}$, which map the reference geometry onto the other shapes, for the training and test sets, respectively, computing the rSVD and the corresponding relative reconstruction errors as described above. The resulting mean and maximum errors are shown in Figure 12. The mean reconstruction error for the displacement field stays almost constant when increasing the number of samples from 512 to 4096, which indicates that the geometric parameter space is reasonably well represented using 512 cases. However, the sensitivity of the Navier–Stokes equations to the domain variability causes the velocity snapshots to have a worse decay of the reconstruction error (Figure 10). Concerning the displacement, increasing $r$ reduces the reconstruction error, achieving a mean relative error of $\approx1 \times10^{-3}$. This error is lower than the best mean reconstruction error obtained for the velocity and comparable to the mean reconstruction error for pressure using the same number of modes. The spatial complexity of the velocity snapshot data is thus—as expected—much higher.

To further characterize the snapshot data, we employ the Uniform Manifold Approximation and Projection (UMAP) technique [45]. UMAP is a nonlinear dimensionality-reduction technique that maps high-dimensional data to a low-dimensional embedding while approximately preserving the neighborhood relationships between data points, revealing structures such as clusters, transitions, and outliers that are difficult to identify in the original space. Figure 13 shows the UMAP of the velocity snapshot data in $\mathbb{R}^2$, for the case $r=500$, identifying which cases yield the highest reconstruction errors. It can be observed that cases with high test $\ell^2$-reconstruction errors lie close to the boundaries of the UMAP embedding, while snapshots yielding errors in the 5%–20% range are associated with outliers in the flow pattern distribution, including mild stenoses (cases 1, 32) and mild coarctation (cases 38, 61, 75, 98, 100).

Inter-patient surrogate modeling

This section focuses on surrogate models built on the dataset of CFD solutions across the shape database introduced in Section 2. We consider two architectures for inter-patient surrogate modeling: (i) Graph U-Nets, relying on message passing on the hierarchy of reference graphs/meshes (level 1, level 0, and the

reference geometry’s centerline), and (ii) Transformer U-Nets, relying on convolutional and cross-attention layers on structured grids. These methods will be introduced in § 3.1. One of the most critical aspects is the high dimensionality of the space of possible aortic shapes, which can only be marginally approximated by the parametrization used to generate the training set [11].

For this purpose, we benchmark several AL strategies used to select the most informative samples during training (§ 3.2). Thereafter, we discuss possible approaches to improve the physical consistency of the model relying on the automatic differentiation of the associated PDE residual during training or within a fine-tuning step (§ 3.3). Finally, we assess the out-of-distribution generalization of the surrogates on different databases, considering a dataset of synthetic aortic coarctation patients [5, 12] (generated via statistical shape models) and the Vascular Model Repository (VMR) [46] dataset in § 3.4. Every computational study in this section is performed on a single NVIDIA A100-SXM4-40GB.

Inference with Graph U-Nets and Transformer U-Nets

Dataset preprocessing

All geometries in the SynthAorta dataset are characterized by the same mesh topology by definition, i.e., retaining a one-to-one correspondence between nodes and elements. Consequently, the nodal displacement field for a geometry $i$ is defined as

where $x_j^{\mathrm{ref}}$ and $x_j^i$ denote the positions of node $j$ in the reference and in the considered geometry, respectively. The nodal coordinates of the case-specific geometry can therefore be expressed as $x_j = x_j^{\mathrm{ref}} + d_j$, where no registration is necessary.

Based on the deformed geometry, we subsequently compute additional geometric features for the network. In particular, we construct a five-dimensional centerline encoding for each node by projecting it onto the case-specific, deformed centerlines of the aortic arch and the three supra-aortic branches. The encoding comprises the distance to the centerline, the arc-length coordinate of the closest centerline point, and the three components of the unit vector pointing from the node toward that point.

Each preprocessed training and test sample stores, at each node, the target fields, namely, the velocity $\boldsymbol{u}$, pressure $p$, and WSS, along with the quantities used as network inputs and conditioning:

(i) the case-specific nodal coordinate vector $\boldsymbol{x} = \boldsymbol{x}^{\mathrm{ref}} + \boldsymbol{d}$,

• (ii) the velocity field, provided either on the full domain or only on its boundary, depending on the field to be predicted, and

• (iii) the five-dimensional geometrical centerline encoding described above.

All input and target fields are normalized component-wise to [0, 1] using minima and maxima computed exclusively from the training set.

Concerning the input velocity field, for the network predicting the velocity, we use only the boundary data, i.e., the velocity restricted to the inlet, outlet, and wall nodes (zero-valued on all interior nodes), so that the model reconstructs the interior flow from the boundary values and the geometry. For the pressure and WSS networks, the full velocity field is provided—either the ground-truth field or the one predicted by the velocity network—so that $p$ and WSS are inferred from the complete velocity field and the geometric encoding. This choice is motivated by recent results suggesting that exploiting the correlation between velocity and pressure fields, rather than the correlation between geometry and pressure, can yield a better outcome in the context of surrogate modeling [8].

Architecture specifics

The Graph U-Net is composed of MeshGraphNet layers [47, 48] implemented in PhysicsNeMo [49], see Figure 14. The Transformer U-Net [10] substitutes graph convolutional layers with attention layers applied to a 3D grid obtained by padding the input vector with zeros. Convolutions and transpose convolutions are applied in a U-Net fashion to the 3D grid, first decreasing the size of the 3D encoded grid and then increasing it back to the original size, see Figure 15.

The variables considered in the network differ between the two architectures, motivated by the significantly higher computational cost for training GNN layers.

1. Graph U-Net: the boundary velocity (3 components), the centerline encoding (5 components), and the deformed position (3 components) are concatenated into an 11-dimensional nodal feature vector, mapped to a higher-dimensional embedding through a random Fourier feature (RFF) encoding followed by a feed-forward network (FFN), and fed as node features to

2. Transformer U-Net: the boundary velocity (3 components), the centerline encoding (5 components), and the deformed position (3 components) are split into two separate streams that are encoded independently: a geometric stream gathering the centerline encoding and the deformed position (8 components), which provides the input tokens processed by the network, and a boundary stream

gathering the boundary velocity (3 components), which is injected as a conditioning signal through the cross-attention layers. Each stream is given its own RFF encoding and FFN projection.

In both cases, the per-node input features (see § 3.1.1) are first lifted to a 64-dimensional embedding through an RFF map (with 2 features over 16 logarithmically spaced scales) followed by a two-layer FFN; the same embedding dimension is used throughout the network and as conditioning width.

The Graph U-Net operates directly on the unstructured mesh, organized as a three-level hierarchy: the volumetric reference mesh at levels 1 and 0, and its centerline. The pooling and unpooling operators are parameter-free $k$-nearest-neighbor ($k = 8$) linear interpolators that transfer node and edge features between consecutive levels. Each level applies a residual block composed of three MeshGraphNet sub-networks, each performing two message-passing steps with two-layer multilayer perceptrons, sigmoid linear unit activation functions, and mean aggregation; the node feature width grows along the encoder following the sequence $\{32, 64, 128\}$ and is mirrored along the decoder as $\{64, 32, 16\}$.

The Transformer U-Net instead embeds the node features into a zero-padded $32^{3}$ voxel grid with a base width of 64 channels. A depth-three encoder halves the grid resolution at each stage while doubling the number of channels ($64 \mathbin{\to} 128 \mathbin{\to} 256$). This is followed by a 512-channel bottleneck composed of three transformer blocks, and a symmetric decoder restores the original resolution with skip connections. Each transformer block combines a cross-attention layer (4 heads), whose queries come from a grouped $3 \times3 \times3$ convolution on the volume and whose keys/values come from the boundary conditioning token, with a $1 \times1 \times1$ convolutional FFN expansion ($\times4$) and group normalization.

We train three different networks, one for each of the three target fields (velocity $\mathbf{u}$, pressure $p$, and WSS), with output dimension 3, 1, and 3, respectively. The main hyperparameters and the resulting number of trainable parameters are summarized in Table 3, while Table 4 details the required training time. The Transformer U-Net is faster per iteration despite its larger parameter count, owing to its dense voxel-grid representation compared to message passing on the unstructured mesh of the Graph U-Net.

Component / hyperparameter	Graph U-Net	Transformer U-Net
Spatial representation	3-level mesh hierarchy	padded $32^{3}$ voxel grid
Conditioning mechanism	fused node features	boundary cross-attention
Input embedding dimension	64	64
U-Net depth / levels	3	3
Feature width (enc. stages)	32, 64, 128	64, 128, 256
Bottleneck width	128	512
Core block	MeshGraphNet ($\times3$)	cross-attention + conv-FFN
Message passing / blocks per stage	2	2
Attention heads	-	4 (every block)
Parameters, $\mathbf{u}$ net	$4.36 \times10^{6}$	$22.17 \times10^{6}$
Parameters, $p$ net	$4.36 \times10^{6}$	$22.17 \times10^{6}$
Parameters, WSS net	$4.36 \times10^{6}$	$22.17 \times10^{6}$

	Graph U-Net	Transformer U-Net
Total training time [h]	38.05	24.80
Average time per iteration [s]	0.61	0.40

Table 5 reports the layer-by-layer number of trainable parameters of the two architectures, taking the velocity network ($n_{\mathrm{out}} = 3$) as reference. The pressure network ($n_{\mathrm{out}} = 1$) only differs in the final projection, yielding 4 362 778 parameters for the Graph U-Net and 22 173 121 for the Transformer U-Net, against 4 363 558 and 22 173 251 for the velocity and WSS networks. The input RFF encoder contributes 8 320 parameters per stream (the Graph U-Net uses a single fused stream, whereas the Transformer U-Net uses two, resulting in 16 640 parameters).

Graph U-Net	Parameters	Transformer U-Net	Parameters
Encoder (3 ConvBlocks + pooling)	1 523 520	Input embedding	4 160
Bottleneck block mid_1	1 186 560	Conditioning projection	4 160
Bottleneck block mid_2	1 186 560	Encoder (3 stages)	3 565 760
Decoder (3 ConvBlocks + unpooling)	457 632	Bottleneck projection	131 584
Final projection block	966	Bottleneck (3 transformer blocks)	12 602 880
RFF encoder	8 320	Decoder (3 stages)	5 847 744
Total ($u$ net)	4 363 558	Output norm + projection	323
		RFF encoders (×2)	16 640
		Total ($u$ net)	22 173 251

Inference results

We use the Graph and Transformer U-Nets to infer the different QoIs (velocity, pressure, and WSS), focusing on the results at peak systole, given that the high-fidelity model’s response is most complex here. To benchmark the performance of the surrogate model, we compare the inference results against the rSVD reconstruction error using $r = 1000$ modes in terms of relative $\ell^{2}$ error. Figure 16 displays the results for the velocity inference. The error in the velocity is marginally higher than the rSVD reconstruction error, while large values in the reconstruction error coincide with cases featuring a large loss. Figure 16 also shows the relative PDE residual for both architectures, which also aligns well with the cases showing the largest losses and reconstruction errors. From these results, we conclude that the test cases featuring large losses are those in which both the reduced approximation and the inference are challenging due to the specific set of geometric hyperparameters. Overall, the errors range from 1% to 20%, with a mean loss of ≈ 2.5%. It is also observed that, although the overall performance of the two methods is comparable, the Transformer U-Net slightly yet consistently outperforms the Graph U-Net over all cases. The cases leading to the maximum, median, and minimum relative $\ell^{2}$ error for the velocity field predicted by the Graph U-Net are depicted in Figure 19, which shows that geometric hyperparameters leading to a narrowing of the descending aorta correspond to the cases with the largest errors.

The errors for pressure inference are reported in Figure 17. One can draw the following conclusions: (i) the prediction errors using the Graph U-Net and the Transformer U-Net are similar, (ii) the cases with increased inference error show large rSVD reconstruction errors, (iii) the two considered architectures have a very similar performance across the dataset, and (iv) compared with the velocity inference, the maximum test error is greatly increased to ≈ 100%, while the mean error lies around 2%.

The WSS U-Net results are given in Figure 18. In this case as well, we conclude that test cases performing poorly (high WSS losses) are characterized by high reconstruction errors. The WSS errors show similar trends to the velocity errors, and overall error levels lie around 2%, while the maximum observed error lies at ≈ 10.5%. In the same figure, we also compare the results with those obtained by replacing the input velocity field taken from the respective velocity U-Net with the full-order solution.

This yields a further improvement in the results. On the other hand, if the full-order solution is available, an approximation of the WSS can be readily obtained directly by postprocessing the field, yielding the appropriate standard numerical approximation, which serves as reference within the current work.

From these results, we conclude that Graph and Transformer U-Nets are valuable frameworks for surrogate modeling of aortic flow across a parametrized set of shapes. However, the accuracy of the results strongly depends on the type of geometric parametrization. Increasing the numerical resolution might be beneficial for the overall framework, removing noise from the snapshot data, but it is currently not computationally feasible. The rSVD reconstruction errors and high inference errors in all the quantities are aligned, indicating that certain parameter sets lead to challenging configurations, where the nonlinear physics cannot be captured as easily.

Inter-patient active learning

One of the crucial challenges in inter-patient surrogate modeling is sampling efficiency due to the high dimensionality of the space of possible aortic shapes. AL is specifically devised to address this problem by extending the training set with the most informative samples. In this section, we analyze the performance of different AL strategies for the considered problem.

Active learning strategies

Suppose that an initial training dataset and a pool $\mathcal{P}$ of candidate geometries not included in it are available. For each geometry $i$, we denote by $\hat{u}_i$ the network prediction (Graph U-Net or Transformer U-Net) and by $(u_i, p_i)$ the corresponding CFD solution. The AL strategies are based on ranking each candidate $i \in\mathcal{P}$ according to a certain acquisition score

and then selecting, for a given $N = |\mathcal{S}|$, the batch $\mathcal{S} \subset\mathcal{P}$ achieving the highest score. We benchmark several strategies, considering scoring functions based on worst error, residual, Fisher information, divergence, active subspaces, query by committee, geometry, and random sampling. Note that obtaining $(u_i, p_i)$ by evaluating the CFD model and retraining the network with the newly selected snapshot data are both costly operations. Therefore, selection criteria based solely on the network prediction $\hat{u}_i$ can significantly reduce the computational cost. The considered acquisition scores are detailed below.

Worst error In the first approach, one selects the geometries with the highest errors, which are the least well represented in the current training set. In this case, the acquisition score depends only on the relative error of the predicted velocity with respect to the CFD solution:

Residual Assuming that the full-order solution is accurate enough, the discrete residual associated with the spatial discretization might be close to optimal and hence considered a target. Based on this consideration, the idea of the second approach is that the relative error in the momentum residual, i.e.,

might be a useful selection criterion. For a given solution $(u,p)$, the residual $r(u,p)$ is assembled on the corresponding domain, i.e., the mapped geometry, and the score is defined as the associated momentum equation relative error:

Fisher information This approach considers the amount of information of each training sample with respect to the model parameters. To this end, we consider the normalized energy of the loss gradient. Let $L_i = \lVert\hat{u}_i-u_i\rVert$ denote the prediction error for sample $i$. The information score is then defined as

where the sum extends over all parameter tensors $\theta_l$ of the neural network (NN). The numerator measures the sensitivity of the loss to perturbations of the corresponding parameters, while the normalization by $\lVert\theta_l\rVert^2$ removes the dependence on the parameter scale. Consequently, samples that induce large parameter sensitivities receive a higher information score, indicating that they carry more information for updating the model during training.

Divergence This approach is similar to residual scoring and considers the magnitude of the discrete divergence of the predicted velocity, i.e.,

This score measures the physical conformity of the network output based solely on the predicted velocity and, therefore, entails a low computational cost.

Active subspaces The goal of this approach is to select the best candidates considering the geometric variations that have the strongest influence on the prediction error. The identification of those geometries relies on an active subspace analysis [50]. Each geometry is represented in a lower-dimensional space by a parameter vector $x_i\in\mathbb{R}^{20}$, corresponding to the leading rSVD coefficients of the nodal displacement field. Based on this dataset, a Gaussian-process surrogate $a(i)=\hat{e}(x_i)$ is trained to approximate the velocity prediction error as a function of these geometric parameters. The active subspace is derived from the covariance matrix of the surrogate gradients,

where $M$ denotes the number of sampled geometries. The dominant eigenvector $\mathbf{w}_1$ identifies the direction in the parameter space along which the predicted error varies most strongly. Each candidate geometry is subsequently projected onto this one-dimensional active coordinate, $z_i = \mathbf{w}_1^{T}\mathbf{x}_i$. The projected samples are partitioned using $k$-means clustering, with the number of clusters selected according to the silhouette criterion. Finally, the desired number of samples $N$ is obtained through uniform sampling of all clusters, ensuring that the selected geometries cover the dominant error-sensitive directions in the parameter space.

Query by committee In this approach, the selection is based on an estimate of the epistemic uncertainty of the surrogate model. We employ a query-by-committee strategy based on Monte Carlo DropConnect. For each candidate geometry $i$, $S$ stochastic forward passes are performed, yielding predictions $\widehat{\mathbf{u}}_i^{(s)}$, $s = 1,\ldots,S$, corresponding to different random realizations of the network weights. The acquisition score is defined as the mean predictive variance over all mesh nodes, i.e.,

where $V_i$ denotes the set of mesh nodes of sample $i$, and $\operatorname{Var}_s$ denotes the sample variance over the $S$ stochastic forward passes. A large predictive variance indicates strong disagreement among the committee members and, consequently, high model uncertainty. Samples with the largest acquisition scores are therefore considered the most informative and are preferentially selected for inclusion in the training set.

Geometry To increase geometric diversity within the training set, in this method, candidate geometries are ranked according to their distance from the geometries already contained in the current training set $\mathcal{T}$. Each candidate geometry $i$ is represented by its nodal displacement vector $\mathbf{d}_i$ with respect to the reference configuration. The acquisition score is defined as the mean Euclidean distance to all training samples,

High acquisition scores correspond to geometries distinct from those already present in the training set. Therefore, selecting samples with the largest scores promotes broad coverage of the geometric parameter space and reduces redundancy among the selected geometries.

Random As a baseline strategy, additional geometries are randomly selected from the pool of available samples $\mathcal{P}$, i.e., $\mathcal{S}\sim\mathcal{U}(\mathcal{P})$, such that every candidate geometry has the same probability of being selected. This strategy serves as a reference for evaluating the effectiveness of the other acquisition methods.

Numerical results

To evaluate the different strategies, we define the total improvement as the difference of the relative errors for the velocity field before and after the AL step:

where the superscripts before and after denote the results without and with the additional sample set $\mathcal{S}$. Here, only the Graph U-Net is considered. A positive value of the total improvement denotes an error reduction. The average is taken over the $n=n_{\mathrm{test}}=100$ held-out geometries introduced in § 2.4 (and also over the $n=n_{\mathrm{train}}$ training cases in Figure 20), i.e., over the same set used to report the inference results in § 3.1: no separate validation split is employed, and, thus, the terms validation and test are used interchangeably in the following. Both the initial training sets of size $n_{\mathrm{train}}^{\mathrm{before}}$ and the candidate pool $\mathcal{P}$ are selected from the remaining $n_{\mathrm{train}}=4560$ training geometries, ensuring that no candidate is ever part of the evaluation set.

Notice that the first three strategies (worst error, residual, Fisher information) are label-dependent, i.e., they require the calculation of the corresponding CFD solution, and will act in the following study as an oracle reference. The label-free ones (divergence, active subspaces, query by committee, geometry, random) only require an approximation with the surrogate, and are thus usable in deployment.

Warm-up regime In the first set of tests, the AL strategies are applied to small training sets, $n_{\mathrm{train}}^{\mathrm{before}} = \{100, 200, 400\}$, with a batch of size $N = n_{\mathrm{train}}^{\mathrm{before}}$ that doubles the set, i.e., $n_{\mathrm{train}}^{\mathrm{after}} = \{200, 400, 800\}$. Figure 20 reports the improvement in training and validation compared to the model trained without

the additional batch. Error- and residual-driven acquisition (worst error, residual, divergence) degrades validation at the smallest batches ($\epsilon_u$ improvement $-0.036$ for worst error and $-0.059$ for residual at $N = 100$), as it oversamples geometric outliers, and recovers to marginally positive values only for the largest batch size $N = 400$. The diversity/uncertainty strategies (geometry, Fisher information, active subspaces) yield small but consistently positive gains, yet none of the AL strategies beats random selection on validation. For the smallest training set and batch size ($N = 100$), Fisher information and random are comparable ($\approx+0.01$), and for $N = 400$ all strategies collapse to parity ($+0.005$ to $+0.008$). The results therefore suggest that the additional effort associated with the AL strategies does not pay off in a scenario with low data availability.

Fine-tuning regime In the second set of tests, we start from a data-rich model, $n_{\mathrm{train}}^{\mathrm{before}} = 2000$, and add three consecutive batches of size $N = 100$ (i.e., $n_{\mathrm{train}}^{\mathrm{after}} = 2000 \rightarrow2100 \rightarrow2200 \rightarrow2300$). Figure 21 reports the validation total improvement after each round. Gains are two orders of magnitude smaller ($\sim10^{-4}$) than in the warm-up regime, reflecting saturation, but targeting now clearly outperforms random selection (the weakest positive baseline, $\sim10^{-4}$). Residual is the only strategy that consistently worsens the model (negative total improvement at every round), while query by committee decays to nearly zero. Worst error dominates the first round ($4.1 \times10^{-4}$) but fades out in the following stages, whereas Fisher information and active subspaces give the most consistent and growing gains, the latter becoming best by the third round ($6.0 \times10^{-4}$). Thus, in the data-rich regime, gradient- and subspace-based acquisition is preferable.

These results indicate that the benefit of AL in the application at hand depends heavily on the specific regime. During the data-scarce warm-up phase, sophisticated acquisition strategies do not outperform random sampling, whereas in the data-rich fine-tuning phase, gradient- and subspace-based methods provide the largest improvements.

Physical consistency as a soft constraint

The NN surrogate models analyzed in the previous sections are trained to minimize the data misfit, but do not explicitly enforce the governing equations. To incorporate physical consistency during optimization, we augment the training objective with a differentiable residual-based penalty obtained by evaluating the weak form of the incompressible Navier–Stokes equations. The NN predictions are interpreted as finite element DoFs, allowing the residual and its gradients with respect to the predicted fields to be assembled efficiently via automatic differentiation. Within this section, we consider only the Graph U-Net for the sake of brevity, noting that no substantial differences in the results are expected when using the Transformer U-Net.

Velocity-dependent residual

For the physical consistency studies, we restrict ourselves to the surrogate model that predicts the velocity field. To compute the residual for a given solution $u$, we consider the velocity-dependent part of the weak formulation in the interior domain, i.e.,

where $V$ and $Q$ denote the function spaces for velocity and pressure, respectively. In particular, we consider $V \subset H^{1}_{0}(\Omega)$, hence restricting to test functions vanishing on the boundary.

The contributions to the residual of the pressure and time-derivative terms are not computed as a function of $u$, but are collected in a precomputed right-hand side $b$ assembled from the full-order velocity and pressure fields $(u_h^{\star},p_h^{\star})$, i.e.,

and

Combining Equations (22), (23), and (24), the physical consistency penalty is defined as the squared norm of the assembled residual. Thus

where $v_h$ and $q_h$ are the standard discrete test functions from the finite element space living on $\Omega_h$. Note that the physical consistency penalty vanishes when the predicted velocity $u_h$ coincides with the full-order velocity $u_h^\star$. Altogether, the loss function takes the form

where $n_{\mathrm{batch}}$ is the batch size, $\theta$ collects the trainable network parameters, $\hat{u}_{h,i}(\theta)$ denotes the predicted velocity for the $i$-th training geometry, interpreted as a finite element function, and $(u_{h,i}^\star,p_{h,i}^\star)$ is the corresponding full-order solution. The first term is the relative $\ell^2$ misfit $\varepsilon_{u,i}$ defined in Equation (12), here evaluated for the surrogate prediction rather than for the rSVD reconstruction, while the second term is the physical consistency penalty in Equation (25), normalized by the precomputed right-hand side. Both contributions are therefore dimensionless and of comparable magnitude. The weight $\beta\ge0$ balances the deviation from the data against physical consistency, with $\beta= 0$ recovering the purely data-driven reference case. We employ $\beta= 1$; for other possible approaches, see [51].

Remark 3.1 The right-hand side $b$ in Equation (24) is assembled from the full-order velocity and pressure fields, so that the penalty in Equation (25) is, in the terminology of § 3.2.2, label-dependent: it can be evaluated only on samples for which the CFD solution is available, and it vanishes exactly at the full-order solution. This is a consequence of the design of the surrogate models rather than of the formulation itself. The networks introduced in § 3.1 predict the velocity field at a single time instant (peak systole), so that neither the temporal derivative $\partial_t u$ nor the pressure $p$ can be recovered from the prediction alone; both are therefore supplied by the reference solution through $b$.

The surrogate models introduced in § 3.1 predict nodal velocity and pressure values on the shared reference grid. To compute the penalty term, we interpret these fields as the DoFs of a piecewise linear finite element function, i.e., the coefficients of the finite element fields with respect to basis functions with support points on the mesh vertices. In practice, the reference mesh is first deformed onto the particular patient-specific geometry using the nodal displacements predicted alongside the flow fields, and the resulting coefficient vector is passed as the evaluation point to an assembly routine. Specifically, for each cell, we compute the interpolation of $u_h$ and $\nabla u_h$ at the quadrature points from the nodal values via the finite element basis and accumulate the discrete weak residual defined in Equation (25). The element residual is registered with respect to the local nodal velocity values and differentiated by forward-mode automatic differentiation via Sacado [52]. The same pass returns both the global residual vector and its Jacobian, such that the squared residual norm enters the loss as a soft penalty for physical consistency, while the transpose Jacobian-vector product propagates its gradient back into the network. The implementation also takes into account batched predictions from the NN surrogate models.

Verification of convergence in a single geometry

As a preliminary step, we validate the physical consistency loss for a fixed single geometry:

For this convergence study, instead of an NN surrogate model for the velocity, we solve for the velocity field $u_h$ by minimizing the residual in Equation (27). Figure 22 reports the results for different numbers of training epochs and final learning rates $lr$ within the cosine annealing scheduler of the SOAP optimizer [53], showing that the relative $\ell^2$ error and the PDE residual converge. The accuracy could be improved by multi-stage training [54] or with other second-order optimizers. The optimization converges up to an accuracy of 2% in terms of $\varepsilon_u$, which validates the correctness of the implementation of the physical consistency loss and its gradient.

$lr_{\mathrm{last}}$	time in min	time/iter. (avg.) in s	epochs
$10^{-3}$	55	0.662	5000
$10^{-4}$	98	0.592	10000
$10^{-5}$	88	0.529	10000

Effectiveness in training and fine-tuning

The physical consistency loss can either be enforced throughout training or introduced only after a purely data-driven model has converged. We therefore investigate its effectiveness in both scenarios, considering the Graph U-NetU trained from scratch in low- and high-data regimes ($n_{\mathrm{train}}=50$ and $n_{\mathrm{train}}=500$), as well as the fine-tuning of the best pretrained model.

Physical consistency during training We compare the accuracy of the Graph U-NetU trained including the physical consistency term with the baseline model trained only with the data term. In the physically consistent variant, the loss to be minimized has two terms, the relative $\ell^2$ error $\epsilon_u$ and the relative residual $r(u_h)/\lVert b\rVert^2$, as in Equation (26).

The results are reported in Figure 23. The physical consistency loss yields a mild but consistent gain only in the low-data regime, while in the high-data regime the data term already dominates and the two variants are essentially indistinguishable. In terms of wall-clock time, enforcing the physical consistency loss increases the per-iteration cost by roughly $3$–$5\times$ over the purely data-driven baseline. Although the evaluation of the finite element residual inevitably incurs additional computational cost, it should be noted that the timings reported here are representative only of the present implementation, which has not yet been optimized for assembly and automatic differentiation.

Physical consistency during fine-tuning Next, we test the soft physics constraint in the loss of the Graph U-NetU model as a fine-tuning procedure. Starting from the best data-driven Graph U-Net (whose weights are frozen), a lightweight feed-forward correction network is trained with the physical consistency loss to predict a nodal velocity correction that is added to the frozen prediction on the interior nodes, while the boundary and wall velocities are kept at their prescribed values. Thus

	$n_{\mathrm{train}}$	time in h	time/iter. (avg.) in s
U-NetU	50	2.03	1.46
U-NetU-PC	50	5.60	4.04
U-NetU	500	8.44	0.61
U-NetU-PC	500	40.40	2.91

For the fine-tuning studies, the full training dataset with $n_{\mathrm{train}} = 4560$ is employed.

As in the two main architectures, the per-node geometric input (8 components: the 5-dimensional centerline encoding and the 3-dimensional deformed position) is first mapped by an RFF encoding (with 2 features over 16 scales with standard deviation $\sigma= 1$) to a 64-dimensional embedding and refined by a two-layer FFN (64$\rightarrow$64 map with sigmoid linear unit activation function). The correction network is a four-layer multilayer perceptron with three hidden layers of width 512 and sigmoid linear unit activation functions (64$\rightarrow$512$\rightarrow$512$\rightarrow$512$\rightarrow$3). The whole module is trained with the SOAP optimizer and a cosine annealing scheduler with an initial learning rate of $10^{-2}$ and a final learning rate of $10^{-4}$.

The results (Figure 24) demonstrate that the physical consistency loss is considerably more effective for fine-tuning than as a regularizer applied throughout training. Although evaluating the finite element residual increases the per-iteration cost by approximately a factor of 1.8 in fine-tuning, the physics-informed correction network reaches a lower validation error after only 10 fine-tuning epochs than the purely data-driven baseline achieves after 100 epochs, resulting in a lower overall wall-clock time. These findings suggest that the residual, once a surrogate has learned an accurate data-driven approximation, provides a complementary optimization signal that efficiently corrects remaining physical inconsistencies without requiring the entire network to be retrained.

	epochs	time in h	time/iter. (avg.) in s
U-NetU-PC	10	12.38	1.06
U-NetU	100	14.84	0.59

Out-of-distribution generalization

The surrogate models presented and tested in the previous sections have been trained and evaluated on geometries sampled from the SynthAorta parametrization [11], i.e., the training and test sets share the same underlying parametric distribution. To assess their robustness and practical applicability in more general contexts, we now evaluate the performance of the surrogate models on pools of patient-specific aortic geometries that lie outside this distribution. This experiment aims at quantifying the extent to which the learned geometric representation generalizes beyond the parametrized shape space and highlights the limitations of the current training dataset. For the sake of brevity, only results obtained via the fully trained Graph U-NetU model are presented.

To this end, we consider (i) a dataset of synthetic aortic shapes generated via statistical shape modeling from a set of aortic coarctation patients [5, 12] and (ii) aortic shapes obtained from the VMR [46]. The out-of-distribution geometries were registered to the current dataset’s template via the multi-level residual network-based large deformation diffeomorphic metric mapping (ResNet-LDDMM) registration method introduced in [8] to compute the corresponding displacement field from the reference shape.

The registration proceeds from a coarse resolution and progressively refines at finer resolutions. The ResNet-LDDMM framework allows for the computation of diffeomorphic transformations that preserve the topology of the geometries, ensuring that the registered shapes remain anatomically plausible. For both datasets, Figure 25 illustrates the registration process, showing the initial geometry ($t = 0.0$), the intermediate deformations ($t \in\{0.25, 0.5, 0.75\}$), and the final registered geometry ($t = 1.0$).

As a preliminary step, Figure 26 shows the embedding of the reduced velocity and displacement fields of the new geometries in the two-dimensional (2D) UMAP space of the SynthAorta dataset (analogous to what has been shown in Figure 13 for the SynthAorta dataset). Cases of coarctation are concentrated in a few sub-regions, whereas the VMR cases—particularly in the displacement embedding—are outliers away from the main manifold. This confirms that the new geometries are out-of-distribution with respect to the training set.

The results in Figure 27 underline the restricted validity of the parametrization used to describe inter-patient geometries. For the VMR dataset, the median relative $\ell^2$ error lies at 183%, whereas the coarctation dataset yields a median relative $\ell^2$ error of 69%. Both indicate that the parametrization underlying the SynthAorta dataset does not sufficiently cover the full physiological spectrum—let alone pathological configurations.

These findings should therefore be interpreted as a limitation of the available training data rather than of the proposed methodology itself. As larger collections of patient-specific anatomies become available, they can be registered to the common reference configuration and used to enrich the parametrized shape space through data augmentation, thereby extending the domain of validity of the surrogate models. The proposed framework is well suited for such incremental refinement, enabling continuous improvements in robustness and generalization without fundamental changes to the underlying methodology.

Conclusion

This work investigated different aspects of inter-patient surrogate modeling for the simulation of blood flow in the aorta. Building upon the SynthAorta [11] parametrization, we generated a large dataset of CFD solutions, which was then used to train graph- and transformer-based neural architectures for the approximation of relevant QoIs. Besides the predictive accuracy and the performance of the surrogate models, we tested several AL strategies for data augmentation, as well as the impact of including physical consistency constraints during the training stage. Geometric parametrization makes it possible to map functions defined on different geometries, such as numerical solutions, predictions, and residuals, onto a reference domain, thereby enabling the study of accuracy, AL, and physical consistency constraints in a common framework.

Several conclusions can be drawn from the numerical results. First, the parametric definition of the dataset allows defining a common reference geometry and displacement-field representations of all considered shapes. These representations provide an effective basis for inter-patient surrogate modeling, which does not require an additional geometric registration of shapes and automatically preserves point-to-point correspondence across patients. The surrogate achieves a median relative $\ell^2$ error of 2% for the velocity field: between 500 and 1000 rSVD modes would be needed to achieve a comparable reconstruction error, making the use of classical reduced basis methods too expensive. Second, the benefit of AL is strongly regime-dependent. Our results show that sophisticated acquisition strategies do not outperform random sampling during the warm-up phase—when only limited training data are available—whereas gradient- and subspace-based methods consistently improve performance when fine-tuning models that are already accurate. Third, incorporating the governing equations through a differentiable finite element residual is more effective as a fine-tuning objective than as a regularizer throughout training. Once the data-driven surrogate has converged, the residual provides a complementary optimization signal that efficiently corrects the remaining physically inconsistent errors at a lower overall computational cost than prolonged conventional training. Finally, the analysis of accuracy for out-of-distribution shapes demonstrates that the dominant limitation of the present framework is not the surrogate architecture itself, but the geometric support of the underlying parametrization and training dataset. Since these considerations apply to any dataset size—given that computational resources are always limited, problem sizes will be adjusted, and parametrizations extended—these results naturally lead to the question of how best to invest a given computational budget in extending snapshot data to minimize inference loss in geometrically parametrized aortic flow.

The primary contribution of this work is the demonstration of an integrated workflow for shape-informed surrogate modeling of vascular flow, combining synthetic dataset generation, geometry-aware neural architectures, differentiable physics constraints, and data acquisition strategies within a unified framework. Beyond the individual algorithmic components, the study provides practical insight into when additional information is most valuable during surrogate construction: increasing dataset diversity is more beneficial than targeted sample selection during early training, whereas physics-based optimization and informed sample selection become increasingly valuable once an accurate data-driven approximation has been established.

Nevertheless, the present study is subject to some limitations. The accuracy attainable with external cohorts depends strongly on the processing of the vascular surfaces, as the position of the planes truncating the inlet and the outlets, as well as the extent of the remaining vessel segments, condition the quality of the registration onto the common template. Moreover, the surrogates are trained at a single instant of the cardiac cycle and rely on rigid vessel walls and a single combination of boundary conditions per geometry. The presented framework is nevertheless readily extensible as larger collections of patient-specific anatomies become available, offering a scalable route towards more robust and clinically relevant surrogate models for cardiovascular flow simulations.

Future work will focus on validating the surrogates against clinical measurements—such as in vivo 4D-flow MRI acquisitions—representing a necessary step for their application in clinical practice.

Acknowledgements

This work received support from the European Union’s Horizon 2024 research and innovation program under the Marie Skłodowska-Curie project “VascularROM”, grant agreement no. 101209352. During preparation of this manuscript and underlying software development, the authors used Claude Code Opus 5.5 (Anthropic) and GPT-5.6/GPT-6 (OpenAI) to assist with style and language editing and coding. The authors take full responsibility for all content in the paper and the employed software.