2.1. Data dimensionality reduction

The processing of large-scale datasets directly poses a significant computational challenge, requiring the implementation of dimensionality reduction techniques to ensure manageability. Singular value decomposition (SVD) (Sirovich Reference Sirovich1987) is widely used in fluid dynamics to reduce data dimensionality and filter out noise. However, the intraventricular flow simulations considered in this work are high-dimensional, transient and highly nonlinear. They involve moving geometries subject to large temporal deformations and complex spatiotemporal patterns that appear. As a result, more advanced techniques are required to analyse and predict their behaviour effectively. In this study, we apply the high-order singular value decomposition (HOSVD) (Le Clainche & Vega Reference Le Clainche and Vega2017 a), a tensor-based extension of SVD that enables simultaneous decomposition along multiple data dimensions. This method improves compression and noise filtering by decoupling the modes in each dimension. As a result, it can distinguish and eliminate noise or redundancies independently across spatial, temporal and physical-magnitude components. Such enhanced denoising and compression capabilities make HOSVD particularly well-suited for building robust, low-rank representations from transient flow data (Groun et al. Reference Groun, Villalba-Orero, Casado-Martín, Lara-Pezzi, Valero, Le Clainche and Garicano-Mena2025), thereby supporting the construction of accurate and efficient predictive models.

The HOSVD is applied to the fifth-dimensional tensor $\mathbf{V}$ , which in turn is factorized as

(2.1) \begin{equation} \mathbf{V}_{j_1j_2j_3j_4k} \simeq \sum _{p_1=1}^{P_1}\sum _{p_2=1}^{P_2}\sum _{p_3=1}^{P_3}\sum _{p_4=1}^{P_4} \sum _{n=1}^{N} \mathbf{S}_{p_1p_2p_3p_4n} \mathbf{W}^{(v)}_{j_1p_1}\mathbf{W}^{(x)}_{j_2p_2} \mathbf{W}^{(y)}_{j_3p_3}\mathbf{W}^{(z)}_{j_4p_4}\mathbf{T}_{kn}, \end{equation}

where $\mathbf{S}_{p_1p_2p_3p_4n}$ is the fifth-dimensional core tensor, and the orthonormal matrices $\mathbf{W}^{(v)}$ , $\mathbf{W}^{(x)}$ , $\mathbf{W}^{(y)}$ , $\mathbf{W}^{(z)}$ and $\mathbf{T}$ contain the SVD modes for each dimension. A tuneable tolerance $ \varepsilon$ determines the number of retained modes $ N$ , ensuring that $ \sigma _{N +1}^{(j)}/\sigma _1^{(j)} \leq \varepsilon _1$ with $j = 1,\ldots ,5$ . This temporal matrix $\mathbf{T}$ containing the temporal coefficients associated with each mode (of dimensions $N \times K$ ) serves as input for the subsequent dynamic mode decomposition (DMD) computation.

2.2. Calculation of DMD modes

Dynamic mode decomposition (Schmid Reference Schmid2010) provides a ROM approach by expressing complex spatiotemporal data as a linear combination of $ M$ coherent structures, or modes, that evolve in time. The approximation of the dataset $ \mathbf{V}_{j_1}(x,y,z,t_k)$ (where $ j_1=1,\ldots ,3$ indexes the velocity components) is formulated as

(2.2) \begin{equation} \mathbf{V}_{j_1}(x,y,z,t_k) \approx \sum _{m=1}^{M} a_m \, \mathbf{u}_m(x,y,z) \, e^{(\delta _m + i\omega _m)t_k} \quad \text{for} \; k=1,\ldots ,K \, , \end{equation}

where $ \mathbf{u}_m$ denotes the spatial structure of each DMD mode, $ a_m$ is the mode amplitude, $ \delta _m$ is the temporal growth rate and $ \omega _m$ represents the oscillation frequency.

To obtain these modes, DMD assumes a linear mapping between consecutive data snapshots $ \{ \mathbf{v}_1, \mathbf{v}_2, \ldots , \mathbf{v}_K \}$ , posing the relation to be described as

(2.3) \begin{equation} \mathbf{V}_2^{K} \approx \mathbf{R} \, \mathbf{V}_1^{K-1} \, , \end{equation}

where $ \mathbf{R}$ is a linear operator that approximates the underlying dynamics, effectively capturing the temporal progression through the so-called Koopman framework (Mezić Reference Mezić2013; Otto & Rowley Reference Otto and Rowley2021).

Higher-order DMD (Le Clainche & Vega Reference Le Clainche and Vega2017 a) extends standard DMD by incorporating multiple time-delayed snapshots, enhancing robustness and accuracy in capturing complex flow structures. Furthermore, the mdHODMD algorithm preserves the tensorial formulation by applying the higher-order Koopman assumption to the temporal matrix $\mathbf{T}$ extracted in the previous step using HOSVD (see (2.1)) as

(2.4) \begin{equation} \mathbf{T}_{d+1}^K \approx \widehat {\mathbf{R}}_1 \mathbf{T}_1^{K-d}+ \widehat {\mathbf{R}}_2 \mathbf{T}_2^{K-d+1} + \cdots + \widehat {\mathbf{R}}_d \mathbf{T}_d^{K-1}. \end{equation}

This formulation links each snapshot to $d$ previous time steps, capturing intricate temporal dynamics beyond what standard DMD offers. The modified Koopman matrix $\widehat {\mathbf{R}}$ , which contains all the linear operators ( $ \hat {R_1}, \hat {R_2}, \ldots , \hat {R_d}$ ), is then computed to extract the DMD modes. Another tuneable tolerance is introduced for the HODMD; in this work, and following the guidelines discussed in Le Clainche & Vega (Reference Le Clainche and Vega2017 a) and Hetherington et al. (Reference Hetherington, Corrochano, Abadía-Heredia, Lazpita, Muñoz, Díaz, Maiora, López-Martín and Le Clainche2024) which is set equal to the one selected for the HOSVD, such that $ \varepsilon _1 = \varepsilon _2 = \varepsilon$ .

2.3. Prediction

Once the data is decomposed into DMD modes, future-state predictions are obtained by evaluating (2.2) at time instants beyond the training window $ t \gt t_K$ . Specifically, the model is constructed using data within the interval $ t \in [0,\, pT]$ , where $ T$ is the period of the cardiac cycle and $ p$ is the number of cycles employed for training. In this work, we will consider $ p \in \mathbb{R}$ and $p\leq 10$ . The resulting ROM is then used to predict the flow dynamics in a later interval, namely $ t \in [10T,\, 20T]$ , allowing us to assess the long-term predictive capabilities of the method.

The way modes are selected plays a crucial role in ensuring prediction accuracy and stability in the model. Depending on the calibration of the algorithm, spurious modes unrelated to the underlying physics may be captured, typically associated with noise or numerical artefacts. As a result, using all computed $ M$ DMD modes can degrade the quality of the extrapolated prediction. To address this difficulty, we adopt a filtering strategy in which only those modes with growth rates lower than a certain threshold are deemed physically relevant and retained for prediction. This threshold is set with a tuneable growth rate value $ \delta _{tune}$ . Thus, stable modes will correspond to $ -\delta _{tune} \lt \delta _m\lt 0$ with $\delta _{tune}\gt 0$ . Almost all of the selected modes correspond to frequencies associated with the dominant dynamics of the system.

Furthermore, when working with transient simulation data, certain modes, although physically relevant, may exhibit highly negative growth rates, leading to their rapid attenuation over time. This offers the potential to enhance the temporal stability of the model predictions; we also construct an alternative formulation in which we enforce $ \delta _m = 0$ , thereby eliminating exponential decay. This modification prevents excessive damping of key flow structures and enables more reliable long-term forecasts of the intraventricular dynamics.

To assess the accuracy of the ROM, we analyse the absolute error distribution over the entire dataset by computing the histogram of the normalized absolute error, defined as

(2.5) \begin{equation} E(\%) = \frac {|\textbf {v} - \textbf {v}_{\text{pred}}|}{\max |\textbf {v} - \textbf {v}_{\text{pred}}|} \times 100, \end{equation}

where $\textbf {v}$ denotes the reference simulation data and $\textbf {v}_{\text{pred}}$ is the reconstructed prediction. The normalization factor corresponds to the maximum absolute error computed independently for each velocity component, considering all spatial locations and time steps in the domain.

The number of bins in the histogram is determined using Sturges’s formula, applied to the spatial grid,

(2.6) \begin{equation} N_{\text{bins}} = 1 + \log _2(N_{\text{grid}}), \end{equation}

where $N_{\text{grid}} = N_x \times N_y \times N_z$ represents the total number of spatial points. In our case, this results in approximately 20 bins, allowing for an intuitive interpretation: the first bin represents the relative frequency of errors within the range 0 %–5 %, the second within 5 %–10 %, and so forth.

A key advantage of the ROM is its drastic reduction in computational cost. In this study, this increase in computational efficiency is quantified using the speed-up factor (SUF), defined as

(2.7) \begin{equation} \text{SUF} = \frac {t_{{\rm CFD}} \times N_{{\rm CPU}}}{t_{{\rm ROM}}} \end{equation}

where $t_{\text{CFD}}$ is the runtime of the full-order CFD simulation, $N_{\text{CPU}}$ is the number of processors used, and $t_{\text{ROM}}$ is the time required for the ROM prediction.

Finally, we note that the HODMD-based ROM involves several tuneable parameters that require careful calibration. These include the truncation tolerance $\varepsilon$ , which controls the number of modes retained; the window size $d$ , which determines the number of snapshots used in each HODMD window; and the threshold $\delta _{tune}$ applied to the growth rates to distinguish physically relevant modes from spurious ones. To calibrate the ROM, multiple combinations of $\varepsilon$ and $d$ are explored during training. While $\varepsilon$ is adjusted depending on the noise level, $d$ is typically chosen between 15 %–50 % of the total snapshots. Once the physical modes are identified, their growth rates and frequencies are analysed to set an appropriate threshold $\delta _{tune}$ , ensuring that only persistent, physically meaningful modes are retained for prediction. This procedure enhances the robustness and accuracy of the resulting ROM.

To calibrate the HODMD-based ROM, multiple combinations of the truncation tolerance $\varepsilon$ and the window size $d$ are explored. The tolerance $\varepsilon$ is adjusted depending on the noise present in the data, whereas $d$ is typically chosen to be between 10 %–50 % of the total snapshots (Le Clainche et al. Reference Le Clainche, Izbassarov, Rosti, Brandt and Tammisola2020; Hetherington et al. Reference Hetherington, Corrochano, Abadía-Heredia, Lazpita, Muñoz, Díaz, Maiora, López-Martín and Le Clainche2024). After identifying the physical modes, their growth rates and frequencies are analysed to define a threshold that separates persistent modes from spurious ones, ensuring the robustness of the ROM construction.

3.1. Database

To ensure that the idealized geometries represent physiologically realistic left ventricular conditions, particularly the end-systolic volume that corresponds to the minimum volume at the end of systole, we define a set of non-dimensional geometric parameters. These parameters are expressed in table 1 using the base radius $ a$ of the ventricular chamber as the reference length for each model.

Table 1. Non-dimensional geometric parameters for the two idealized ventricular models, using base radius $ a$ as reference length.

A schematic representation of these parameters and their spatial arrangement is provided in figure 1, illustrating the geometry and orientation of each feature. The figure also identifies the symmetry plane A–A’ (located at $y=0$ ) used throughout the analysis (highlighted in blue). Additionally, a probe line (L1), indicated in red, is defined to extract temporal data at specific spatial locations. For the Ideal 1 model, L1 is located at $ x/a = 0, \, y/a = 0$ , and it extends through all the vertical $z$ of the chamber. For the Ideal 2 model, L1 is positioned at $ x/a = 0, \, y/a = 0$ , and again it takes all the chamber height.

Figure 1. Representation of both idealized LV models: (a) Ideal 1 and (b) Ideal 2. Blue, symmetry plane A–A’ used for two-dimensional visualization; red, probe line L1 to represent temporal data.

Using the two idealized geometries of the LV, we conducted CFD simulations to investigate intraventricular flow dynamics. Simulations were carried out using the solver Ansys Fluent (Ansys, Inc. 2023). In both cases, blood was modelled as an incompressible Newtonian fluid with constant density $\rho = 1060\,\text{kg m}^{-3}$ and dynamic viscosity $\mu = 0.004\,\text{Pa}\, \text{s}$ .

To characterize the flow regime, we computed the Reynolds number ( $\textit{Re}$ ) based on the inlet tube diameter $D$ (non-dimensionalized as $D/a$ ) and the peak inlet velocity $V_{\text{max}}$ over the cardiac cycle. The resulting Reynolds numbers were approximately $\textit{Re} \approx 5500$ for the Ideal 1 model and $\textit{Re} \approx 4900$ for Ideal 2. Despite these values indicating transitional or potentially turbulent flow, we assumed a laminar regime, in line with previous studies (He et al. Reference He, Han, Zhang, Shah, Kaczorowski, Griffith and Wu2022; Tagliabue, Dedè & Quarteroni Reference Tagliabue, Dedè and Quarteroni2017), where coherent vortex structures, such as the main VR, were accurately captured under laminar assumptions.

A zero-velocity initial condition is set in the simulation pipeline for both models. Moreover, for both simulations, pressure-based boundary conditions were set at the inlet and outlet. These boundaries dynamically alternated between open and closed states to represent the cardiac cycle: during diastole, the inlet was open and the outlet treated as a wall; during systole, the roles were reversed. This boundary-switching was complemented by the dynamic motion of the ventricular wall, imposed by prescribing an interpolated sequence of meshes corresponding to each time instant. This approach induces pressure gradients within the chamber that drive blood in and out, thereby mimicking physiological behaviour.

The simulation framework and its assumptions were extensively validated in previous studies (Lazpita et al. Reference Lazpita, Garicano-Mena, Paniagua, Le Clainche and Valero2024 a,b, Reference Lazpita, Neidlin, Garicano-Mena and Clainche2025), including convergence analyses and comparisons with benchmark data. All numerical simulations were performed on the high-performance computing cluster Magerit, hosted at CeSViMa (Centro de Supercomputación y Visualización de Madrid). The computations were parallelized using 40 message passing interface (MPI) tasks within a single node.

For the rest of the manuscript, time is non-dimensionalized using the period $ T$ , such that $ t^* = t/T$ . This normalization allows for a consistent comparison of corresponding time instants across both simulation databases. It is also important to note that, for the first geometry, the diastolic phase takes place within the interval $ t^* \in [i,\, i + 0.67]$ , followed by the systolic phase in $ t^* \in [i + 0.67,\, i + 1]$ . In contrast, for the second geometry, diastole spans $ t^* \in [i,\, i + 0.8]$ , and systole occurs within $ t^* \in [i + 0.8,\, i + 1]$ , where $ i$ denotes the cycle index.

A key aspect of intraventricular hemodynamics is the formation and evolution of a VR during diastole, triggered by the inflow through the mitral valve (Nagargoje et al. Reference Nagargoje, Lazpita, Garicano-Mena and Clainche2025). This structure undergoes tilting and eventual breakdown due to instabilities, and its dynamics are considered a primary indicator of pumping efficiency. Accurately capturing the VR behaviour is, therefore, essential.

Figure 2 illustrates representative time instants for both geometries. The green isocontours show the Q-criterion at a level of 2500, highlighting vortex cores, while the colormap on the symmetry plane displays the out-of-plane vorticity in the range $[-100, 100] \: \text{s}^{-1}$ . Notably, the VR in Ideal 1 dissipates earlier than in Ideal 2; by $ t^* = 0.25$ , the ring in Ideal 1 is already breaking up, whereas Ideal 2 still maintains a coherent and smooth structure. Even at $ t^* = 0.40$ , the second model, though deformed, maintains a closed VR, while the first model’s ring has fully dissipated.

As shown in Lazpita et al. (Reference Lazpita, Neidlin, Garicano-Mena and Clainche2025), the vortex breakdown process is strongly influenced by ventricular geometry. In the Ideal 1 configuration, which features a narrower chamber, the vortex interacts with the ventricular wall at an earlier stage. This early interaction triggers the breakdown process before the VR can reach the apex. In contrast, the Ideal 2 model, with a wider chamber, allows the VR to fully develop, travel towards the apex, and break down as a combination of a reduction in velocity and wall interaction.

Figure 2. Representative snapshots of the intraventricular flow fields for Ideal 1 (a) and Ideal 2 (b). Green, Q-criterion isosurface at 2500; background, vorticity in the symmetry plane; range $[-100, 100]$ .

These contrasting behaviours can be attributed to two main factors:

1. (i) the narrower geometry of Ideal 1 causes the vortex to interact more strongly with the ventricular walls, preventing it from reaching the apex and accelerating its breakdown;

2. (ii) Ideal 2 exhibits a lower ejection fraction and volume change, resulting in reduced pressure gradients and slower flow dynamics.

These differences in VR dynamics will be key for assessing the ROM’s ability to capture flow physics across varying conditions.

Each simulation database consists of 20 cardiac cycles, with flow fields sampled at a time interval of $ \Delta t^* = 0.05$ : therefore, these are 20 snapshots per cycle and a total of $ K = 400$ snapshots. For training the ROM, we use data from the first 10 cycles, which include the initial transient phase. The remaining 10 cycles are reserved for validation, as they capture the dynamics once the system has reached temporal convergence. This distinction is further supported by the application of DMD, which reveals that growth rates stabilize only in the latter cycles.

Although it is common practice in the literature to simulate up to six cycles, we shall find that the flow remains in a transient state at that point (Canè et al. Reference Canè, Delcour, Luigi Redaelli, Segers and Degroote2022; Grünwald et al. Reference Grünwald, Korte, Wilmanns, Winkler, Linden, Herberg, Groß-Hardt, Steinseifer and Neidlin2022; Korte et al. Reference Korte, Rauwolf, Thiel, Mitrasch, Groschopp, Neidlin, Schmeißer, Braun-Dullaeus and Berg2023). This highlights the importance of longer simulations to ensure the reliability of the ROM and the physical relevance of the extracted modes.

The cases listed in table 2 are analogous for both LV geometries and thus we maintain the same denomination for consistency and brevity. Each case corresponds to a specific time window and number of cycles, which are used either for training or validation of the ROM. We also report the tuneable growth rate value selected for each case, which enables the identification of persistent modes that represent the main dynamics of the system, effectively discarding the spurious modes extracted by HODMD.

Table 2. Summary of simulation cases used to build and evaluate the ROM discussed in § 2, including the number of cardiac cycles $p$ , normalized time intervals $ t^*$ , snapshot count $K$ , their intended purpose and the selected tuneable growth rate parameter $ \delta _{tune}$ .

3.2. Database sensitivity assessment

Despite careful simulation set-up, the complex nature of moving-wall CFD simulations introduces inherent variability across cycles, even setting the boundary conditions perfectly periodic (Chnafa et al. Reference Chnafa, Mendez and Nicoud2014). While major flow structures, such as the primary VR, secondary vortices and outflow jets, are consistently reproduced, smaller-scale features might vary.

We quantify the variability across cardiac cycles by analysing the cycle-to-cycle variation of the total kinetic energy (TKE), defined at each non-dimensional time instant $ t^*$ as

(3.1) \begin{equation} \mathrm{TKE}(t^*) = \int _{V(t^*)} \frac {1}{2}\rho \mathbf{v}^2\, \text{d}V, \end{equation}

where $ V(t^*)$ is the time-dependent ventricular volume, $ \rho$ is the fluid density, and $ \mathbf{v}$ is the velocity vector. The TKE integrated over each cardiac cycle is computed as

(3.2) \begin{equation} \mathrm{TKE}_i = \int _{i-1}^{i} \mathrm{TKE}(t^*)\, \text{d}t^*, \quad i = 1, \ldots , 20, \end{equation}

where $ i$ , again, denotes the index of the cycle. To measure convergence and variability, we evaluate the relative root mean square error (RRMSE) of each cycle with respect to the last one,

(3.3) \begin{equation} \mathrm{RRMSE}_{\mathrm{TKE}}(\%) = \frac {|\mathrm{TKE}_{20} - \mathrm{TKE}_i|}{|\mathrm{TKE}_{20}|} \times 100, \end{equation}

where $ \mathrm{TKE}_{20}$ is the integrated kinetic energy for the 20th cycle, considered as the reference.

Figure 3 illustrates the cycle-to-cycle variability in TKE error, expressed as the RRMSE for both Ideal 1 and Ideal 2 geometries. Remember that the simulation was initialized from a zero-velocity condition. After the initial transient, the RRMSE remains below 4 % across all subsequent cycles, indicating good reproducibility of global flow features. However, local flow characteristics, such as instantaneous velocity fields, may still exhibit larger discrepancies due to intrinsic variability in the numerical solution (Chnafa et al. Reference Chnafa, Mendez and Nicoud2014; Canè et al. Reference Canè, Delcour, Luigi Redaelli, Segers and Degroote2022). As such, predictions obtained from ROMs should be interpreted within a baseline uncertainty margin of approximately 5 %–10 %.

Figure 3. Cycle-to-cycle RRMSE of TKE for Ideal 1 and Ideal 2, visualized using a grouped bar plot.

Figure 4. Spectrum of DMD modes for both idealized models: normalized amplitude versus frequency (a) and absolute value of growth rate versus frequency (b). The dashed line at $ \delta _{tune} = 5 \times 10^{-2}$ in (b) separates persistent physical modes from transient or spurious ones.

3.3. True frequency analysis

Figure 4 presents the DMD spectra obtained from a representative HODMD calibration applied to both Ideal 1 and Ideal 2 cases. Figure 4(a) shows the normalized amplitude of each mode with respect to frequency, while figure 4(b) displays the corresponding absolute values of the growth rates. Although only one calibration is depicted here for clarity, the modes shown have been previously identified as robust through an exhaustive sensitivity analysis involving the HODMD parameters $d$ and $\varepsilon$ , following the methodology described in Le Clainche et al. (Reference Le Clainche, Izbassarov, Rosti, Brandt and Tammisola2020) and Hetherington et al. (Reference Hetherington, Corrochano, Abadía-Heredia, Lazpita, Muñoz, Díaz, Maiora, López-Martín and Le Clainche2024).

The amplitude spectra in figure 4(a) reveal the dominant frequencies that appear consistently across both geometries. These include the fundamental frequency $\omega = 2\pi$ and its harmonics, such as $4\pi$ , $6\pi$ and up to approximately $14\pi$ , indicating strong periodicity and physical relevance of these modes. Lower-amplitude modes, scattered throughout the spectrum, are not persistent across calibrations and are interpreted as numerical noise or transient artefacts.

Figure 4(b) provides additional insight through the analysis of growth rates $|\delta _m|$ . Robust physical modes tend to have small or nearly zero-growth rates, indicating long-term persistence throughout the simulation. A threshold of $ \delta _{\text{tune}} = 5 \times 10^{-2}$ (shown as a dashed line) distinguishes these physical modes from the spurious ones, which exhibit larger decay or amplification rates and lack consistency. This threshold was established based on a qualitative shift observed in the spectral structure during the calibration study: below this value, the modes were consistently identified and physically meaningful, while above it, the spectrum became irregular and inconsistent.

These results confirm that the dominant dynamics of both Ideal 1 and Ideal 2 models are governed by a limited number of coherent structures oscillating at fundamental and harmonic frequencies, in agreement with the known periodicity of the cardiac cycle.

While the current study focuses on filtering spurious modes to construct a robust ROM, a detailed analysis of the retained modes and their physical interpretation is presented in a complementary work (Lazpita et al. Reference Lazpita, Neidlin, Garicano-Mena and Clainche2025), where the underlying flow dynamics in idealized LV models are extensively characterized. In that study, the leading modes were shown to capture the formation, tilting, and breakdown of the primary VR generated during early diastole, associated with the fundamental frequency $\omega = 2\pi$ , as well as the appearance of a secondary VR during the late diastolic inflow ( $\omega = 4\pi$ ). Higher-order modes revealed nonlinear interactions between these frequencies, describing the onset of vortex breakdown and the emergence of smaller-scale structures. This modal hierarchy thus provides a physically consistent picture of the intraventricular flow, linking coherent structures to characteristic temporal dynamics.