A reduced order model formulation for left atrium flow: an atrial fibrillation case

A data-driven reduced order model (ROM) based on a proper orthogonal decomposition-radial basis function (POD-RBF) approach is adopted in this paper for the analysis of blood flow dynamics in a patient-specific case of atrial fibrillation (AF). The full order model (FOM) is represented by incompressible Navier–Stokes equations, discretized with a finite volume (FV) approach. Both the Newtonian and the Casson’s constitutive laws are employed. The aim is to build a computational tool able to efficiently and accurately reconstruct the patterns of relevant hemodynamics indices related to the stasis of the blood in a physical parametrization framework including the cardiac output in the Newtonian case and also the plasma viscosity and the hematocrit in the non-Newtonian one. Many FOM-ROM comparisons are shown to analyze the performance of our approach as regards errors and computational speed-up.


Introduction and motivation
Atrial Fibrillation (AF) is the most common type of cardiac arrhythmia, affecting around 1.5% of the population (more than 35 million people worldwide [9]).Its incidence seems to be correlated with age, since 8% of individuals above 80 years old are affected [35].An AF episode causes the Left Atrium (LA) to contract in an irregular and ineffective way, being typically triggered by irregular electrical impulses coming from the Pulmonary Vein (PV) roots.This abnormal contraction pattern seems to be related with stroke incidence due to thrombus formation within the Left Atrial Appendage (LAA) [69].
The LAA is a cavity which results from LA embryonic development, having a protruding and trabeculated morphology [2].When a patient develops AF, the LAA natural contractility reduces dramatically, making it prone to thrombi formation [36,60].Due to this, the LAA has received a lot of attention from both clinical and biomedical engineering fields.Although a large number of studies have attempted to delve deeper into the relationship between stroke risk and LAA morphology in AF patients, the underlying mechanisms are still not well understood.Some of these previous studies have suggested certain morphologies to be associated with a lower risk of stroke [20,48,70], while others could not find a correlation between the two phenomena [42,53].Other clinical studies associated stroke risk with certain atrial geometric parameters such as LAA volume [7,43], LAA ostium area [7,47,48], number of lobes [71], and LAA depth [7].Recent studies also drew attention to the importance of the atrial flow [47], and the PV morphology [55].Although everything seems to point to the existence of a mechanistic relationship between the thrombosis risk and the LA anatomy and flow, this one remains unknown at this time.
Recent advances in medical imaging have made it possible to apply Computational Fluid Dynamics (CFD) techniques to the study of LA flows.Early research on LA flow dynamics were part of whole left-heart simulations [18,64].More recent CFD studies have presented numerical analyses of the LA flow patterns [18,46,54], focused on investigating blood stasis, as it is considered a necessary thrombogenic factor.Some of them also studied specifically the LAA stasis in AF conditions [13,23,24,30,31,50].All of them have helped to provide insights into the AF phenomenon and calculate otherwise inaccessible parameters such as residence times and shear stress that are related to blood stasis.Other studies have showed that the residence time and flow patterns for flexible-wall and rigid-wall simulations are very similar in case of impaired atrial function, especially when both the reservoir and booster functions are decreased [23,31].An interesting approach is used by Dueñas-Pamplona et al. [22], who developed a morphing technique to study the risk of long-term stasis due to geometrical parameters.To reproduce the most critical AF case (in the absence of atrial contraction) the atrium was kept rigid, regardless of the atrial function at the time of medical imaging.Results show the enormous influence of cardiac output in the blood age indices, and the relatively minor role played by the PV orientation.To our knowledge, this is the only study that attempts to parametrize the LA flow problem.On the other hand, Gonzalo et al. [37] have drawn attention to the fact that the non-Newtonian blood rheology can impact the left atrial stasis in patient-specific simulations, submitting that hematocrit-dependent non-Newtonian blood rheology should be considered when calculating patient-specific blood stasis indices by CFD [37].
Reduced Order Models (ROMs) [10,11,12,39,58] have become increasingly important in hemodynamics applications due to the complex and multiscale nature of the cardiovascular system.ROM techniques allow for the creation of simplified models able to capture the essential features of the system and to significantly reduce the computational cost with respect to the standard CFD models based on classic discretization techniques, such Finite Volume (FV) or Finite Elements (FE) (hereinafter referred to as Full Order Model (FOM)).ROMs can also aid in the development of personalized medicine, as they enable the creation of patient-specific models that can help clinicians to better understand and treat cardiovascular diseases.
Classic projection-based ROMs were already adopted to speed up patient-specific cases or idealized cardiovascular benchmarks.A Proper Orthogonal Decomposition (POD)-Galerkin strategy is employed in combination with a FV full order solver in [15].The aim is to predict the pressure drop along an idealized vessel in a geometrical parameter setting.The same ROM approach but within a FE environment is adopted in [4,5,72] for the study of the blood flow patterns in patient-specific configurations of coronary artery bypass grafts where a physical parametrization involving the Reynolds number as well as an efficient centerlines-based geometrical parametrization are employed.
More recently, the combination of data-driven ROMs and FV method is becoming particularly appealing, due to the diffusion in the biomedical engineering community of commercial codes relying on FV schemes: see, e.g., [8,16,65].In [32] a POD with interpolation by Radial Basis Function (RBF) approach is adopted for the investigation of the hemodynamics in the aortic arch in presence of a left ventricular assist device.The authors report a speed up of O(10 6 ) associated with an error less than 15% which represents a promising result in such a direction.For this reason, in [6,61,62] a similar approach where the interpolation is carried out by feed-forward Neural Network instead of RBF is proposed for the analysis of the blood flow patterns in coronary artery bypass grafts.While in [6] only one physical parameter is introduced (the Reynolds number), in [61] a geometrical parametrization setting (with respect to the diameter of an isolated stenosis) is also considered.In such works, a speed up of O(10 5 ) and an average error below 5% are obtained.
Concerning the problem addressed in this paper, some recent works [56,59] use machine learningbased models to infer LAA blood stasis from LAA geometry stasis.Their goal is to elaborate the huge amount of information coming from the data in order to identify patients with AF at the highest risk of thrombus formation.Unlike these works, mainly based on data processing, our aim is to retain the physical problem and build a cooperation between ROM, CFD and data-driven techniques in order to build an efficient computational tool able to achieve faithful solutions.Following this research line and its encouraging results, the present work tries to extend the use of data-driven ROM approaches to the study of the blood flow in the LAA portion in a physical parametrization setting.This is of course a more complex case and a step forward with respect to vessel-like structures addressed in the previous works [6,32,34,61,62].Indeed, to the best of our knowledge, this work represents the first study about the application of ROM to the haemodynamics in the LA.For sake of completeness, we mention [28] where a ROM for the propagation of the electrical signal in the heart is analyzed.However, the authors do not address the fluid dynamics of the problem and furthermore use an idealized LA domain while a patient-specific one is used in this work, which represents an additional difficulty introduced in our research.
The paper is organized as follows.In Section 2 and 3, the FOM and the ROM adopted are explained in detail, along with the hemodynamics indices of interest and the techniques chosen for each algorithm.Then, Section 4 is reserved for error and efficiency analysis of our ROM approach.It also includes some FOM-ROM qualitative comparisons and clinical considerations on the patterns obtained.Finally, Section 5 is dedicated to draw conclusions and some possible extensions of this work.

The Full Order Model
The FOM employed in this paper is similar to the one proposed by Dueñas-Pamplona et al. in [24].It consists of parametrized Navier-Stokes equations for the blood flow in a patient specific domain Ω over a time interval of interest (t 0 , T ]: where v = v(x, t; µ) and p = p(x, t; µ) are the velocity and the pressure depending on the physical parameter vector µ, respectively.In addition, ρ = 1050 kg/m 3 is the blood density.Problem (1) is endowed with initial data v(x, t 0 ; µ) = 0 and suitable boundary conditions reported in Section 2.1.
T is the Cauchy stress tensor whose constitutive relation has the form: where the deviatoric component, T d , depends on the fluid model.It is known that, even if plasma is approximately Newtonian, whole blood could exhibit significant non-Newtonian features [17].Many works compare Newtonian and non-Newtonian models, showing that the Newtonian one is in general admissible.However, some differences can be found locally for some portions of the domain and/or during specific time instances of the cardiac cycle [41].For such a reason, in this work both models, Newtonian and non-Newtonian, are considered.For a Newtonian fluid, the tensor T d is: with µ = 0.0035 Pa•s the constant dynamic viscosity and D(u) = ∇u+∇u T 2 the strain rate tensor.On the other hand, the non-Newtonian behavior of the blood is handled with the Casson's model [21]: and where J 2 is the second invariant of D(v), η = η/(1 − H) 2.5 with η the plasma viscosity and H the hematocrit, and τ y = (a 1 + a 2 H) 3 where for blood a 1 = 0 and a 2 = 0.625 [29].For more details we refer to [26].We note that the plasma viscosity η and the hematocrit H will be treated as parameters in our ROM framework.Equations above are solved using the FV method implemented in OpenFOAM(R) 2204 [1].Direct numerical simulation of the blood flow is used here assuming that the flow is not turbulent.Second order schemes are adopted for space discretization and Euler implicit time scheme is used for time discretization.More details about the FV discretization of Navier-Stokes equations can be found in [32].
2.1.Patient-specific geometry and boundary conditions.In this research, the same patient geometry reported by Dueñas-Pamplona et al. in [22] is used as test configuration.The patient had a history of paroxysmal AF, but no previous stroke or transient ischemic attack, and underwent CT-imaging and Doppler transesophageal echocardiography at the Puerta de Hierro Hospital in Madrid.The images were obtained from the left atrial endocardial surface, including the LAA and the PVs, during normal sinus rhythm, using a SIEMENS Sensation 64 scanner with specific scanning parameters.The resulting DICOM images were manually segmented with the aid of the non-commercial code 3D-Slicer to extract the three-dimensional surface of the endocardium.The PVs were removed beyond the first bifurcation using Autodesk MeshMixer, and the LA endocardial surface geometry, including the LAA, was provided.In Figure 1 (left) the patient-specific geometry is reported.
Walls are assumed to be completely rigid due to the patient's condition.Regarding inflow and outflow boundary conditions, those reported in [22] are also used.They consist of Mitral Valve (MV) velocities during sinus rhythm after CT imaging.These blood velocities were used to establish the simulation boundary conditions, assuming each PV carries out the same flow-rate.In order to clarify the inflow and outflow boundaries, in Figure 1  2.2.Hemodynamics indices.Now we are going to introduce some relevant indices associated with the stasis of the blood flow [23].They represent fundamental quantities for the medical community because high residence times of the blood flow are related to thrombus formation.2.2.1.Age distribution.Following [63], the first and second moments of age, m 1 and m 2 , are computed as follows: where ϕ(x, t; µ) is the blood age distribution.In particular, m 1 represents the mean age of blood and describes the time needed for a given blood particle to reach another position in the computational domain, so large values of m 1 denote stasis of the blood.The second moment m 2 does not have any physical meaning.However we have decided to consider it as a further variable to be taken into account to show the versatility of our ROM approach.It should be noted that once m k is known, it would be possible to compute also ϕ (6): see, e.g., [40,63].For a laminar flow, we have for k = 1, 2: and initial data m k = 0, where n is the unit normal outward vector to the boundary, Γ i denotes the inflow boundaries (i.e. the PV inlets, see Figure 1), µ m k = 10 −10 kg/(m • s) is the mass diffusivity of the moment m k [23] and m 0 = 1 [63].
Another variable related to the blood age is the washout.It is computed by solving a scalar transport equation similar to (7) but without the source term and initial data φ = 1, where µ φ = 10 −10 kg/(m • s).The washout is the residual value of φ as the dynamics evolves over the cardiac cycle.It implies zones with high stasis and therefore high age.More details could be found in [40,63].Equations ( 7) and ( 8) are coupled with system (1) but, in order to reduce the computational cost, we adopt a segregated algorithm, so they are solved after the system (1).In other words, m k and φ are treated as passive scalars.Equations ( 7) and ( 8) are employed by using second order schemes for space discretization and Euler scheme for time discretization.

2.2.2.
Wall Shear Stress and Oscillating Shear Index.The Wall Shear Stress (WSS) can be defined as follows: The interest is on the Time Averaged Wall Shear Stress (TAWSS) representing the mean effect of the WSS on the entire cardiac cycle: Low values of TAWSS correspond to stasis regions.Another important quantity is the Oscillating Shear Index (OSI): It is related to the oscillations of the flow and ranges from 0, when the flow is unidirectional, to 0.5, when the direction of the flow is totally reversed.The OSI is a useful indicator in cardiovascular problems because it is correlated with the intimal thickness of the wall and the restenosis process [44].

The reduced order model
The basic assumption of ROM for a partial differential equations problem depending on time t and parameter vector µ is that any solution can be represented as a linear combination of a reduced number of global basis functions, that depend exclusively on space x, with the weights of the linear combination depending only on t and µ.For a generic variable Φ this is written as: where Φ rb is the reduced order approximation of Φ, L is the number of basis functions, the ℓ l are the basis functions and the α l are the weights of the linear combination (the so-called modal coefficients).
In this work, we are interested in the reconstruction of all the hemodynamics indices introduced in Section 2.2, so we have Φ = {m 1 , m 2 , φ, TAWSS, OSI}.It should be noted that TAWSS and OSI are steady-state variables.In this case, equation ( 12) becomes: However, for the sake of good order, in the following explanation we refer to a time and parameter dependent variable.We use the POD-RBF technique, which is divided into the following two phases: Offline: given a set of physical parameters values and time instances, the corresponding high fidelity solutions (the so-called snapshots) are computed and collected into a matrix.The POD algorithm is used to extract the reduced basis space from the snapshots matrix.Then, the snapshots are projected onto the POD space by obtaining the corresponding modal coefficients.Finally, a RBF interpolation algorithm is used to compute a map between the parameters and the modal coefficients.Online: given a set of new physical parameters values, the corresponding modal coefficients are computed via the RBF function and the approximated solution is recovered as a linear combination between these coefficients and the POD reduced basis (see equation ( 12)).The offline-online procedure is summarized in Algorithm 1.The implementation of the ROM relies on the Python library EZyRB [19].
The POD [3,38,45,66,68] is one of the most common techniques used to extract the essential information from the space generated by the solution manifold.Let N s = N t • N p be the dimension of this space: the goal of the POD is to construct a reduced basis space of dimension L ≪ N s which is optimal in the least-square sense.More precisely, the POD algorithm builds the reduced basis space which minimizes the quantity among all L-dimensional subspaces V rb spanned by the FOM solutions [39].
Let N h be the number of the mesh cells.We collect the FOM solutions into the snapshot matrix S Φ given by whose dimension is N h ×N s .Since S Φ is usually not squared, we introduce its rank R ≤ min{N h , N s }.
By applying the Singular Value Decomposition (SVD) to S Φ , we can rewrite it as Ns×Ns are orthogonal matrices whose columns are the left and right singular vectors, respectively, and Σ ∈ R N h ×Ns is a diagonal matrix with R non-zero real singular values As the rank R is typically large, we are now concerned with reducing the size of the problem.We rely on the Schmidt-Eckart-Young theorem [25], which states that the first L left singular vectors of S Φ are the POD bases of rank L, with L < R. Hence, the POD bases are the first L columns of the matrix L. In particular, the l-th column of L is the eigenvector associated with Cℓ l = σ 2 l ℓ l , where C = (S Φ ) T S Φ is the snapshot correlation matrix.Therefore, ℓ l is given by Then, the POD bases, also known as modes, are collected into the matrix It remains only to choose properly the value of L. A common choice is to define it as the smallest integer such that for a given threshold ε on the cumulative energy of the eigenvectors.Once defined the POD basis matrix B, the reduced solution Φ rb that approximates the truth one Φ is given by Φ rb (x, t i ; µ j ) = L l=1 α l (t i , µ j )ℓ l (x), for i = 1, . . ., N t and j = 1, . . ., N p and where α l (t i , µ j ) = B T Φ(t i , µ j ) l is the l-th modal coefficient.

Radial basis function interpolation.
The last step of the offline phase consists of the definition of a map F : T × K → R L from the parameters space to the modal coefficients one.To this end, we use the RBF interpolation [14,27,67].In general given a set of N s couples (x i , y i ) and a query point x, the RBF is defined as where w j are weights to be determined and P (x) is a polynomial required for stability reasons.For simplicity, we assume that P is of degree 1.By adding the conditions Ns j=1 w j = 0 and Ns j=1 w j x j = 0 the weights w j and the polynomial P are uniquely determined.As explained in Section 3.1, the projection of the snapshots onto the POD space gives us the modal coefficients α(t i , µ j ) = {α l (t i , µ j )} L l=1 corresponding to the samples (t i , µ j ), with i = 1, . . ., N t and j = 1, . . ., N p .Thus, in our case, we have that (x, y) ≡ ((t, µ), α(t, µ)).

Numerical results
In this section, we investigate the performance of our ROM model.For what concerns the validation of the FOM model in an artificial setup and the independence of the numerical results from the mesh for the patient-specific case of interest, we refer to Appendix A.
We consider a tetrahedral mesh inside the domain Ω and on its boundary ∂Ω, consisting of 542.560 cells and 79.708 cells respectively, i.e. the Grid 2 reported in Table 2.In Figure 2 we show a sketch of the mesh and we highlight the LAA portion on which we will focus our tests (the reason for this choice is discussed in Section 1).Since we are going to analyze the regime behaviour, i.e. downstream of the transient effects, we collect the FOM solutions corresponding to the fourth cardiac cycle whose period is 1.07 s.Therefore, the effective time interval is (t 0 , T ] = (4.33,5.4] s which, for sake of convenience, we report to (0, 1].We ran the FOM simulations with a time step ∆t = 0.01 s collecting N t = 107 time dependent snapshots.To unify the notation, where necessary, hereinafter we will use the superscript N to refer to the Newtonian model and the superscript C to refer to the Casson's one.In the Newtonian case we have one only parameter represented by the scaling factor f ranging in the interval [0.5, 1.5] (see Figure 1), i.e. µ N = f .On the other hand, in the non Newtonian case, we also consider the plasma viscosity η ∈ [1.5e − 05, 1.7e − 05] and the hematocrit H ∈ [35,50], i.e. µ C = {f, η, H}.Specifically, the discrete set of µ N consists of 20 points obtained by a uniform sampling procedure.Instead, for the Casson's model, we get a distribution of 30 points for µ C using Latin hypercube sampling [51].
For both the models, we split the initial database in a training set K train ⊂ K and in a validation set K test = K \ K train .All the snapshots belonging to the training set are stored in the matrix S Φ .The validation set is used to assess the accuracy of the ROM solution.In order to compare the results obtained with the two models, we choose similar test values for the scaling factor f , which is the only shared parameter.Therefore, we set where µ ∈ K test and ∥ • ∥ is the Frobenius norm.For the steady state variables, TAWSS and OSI, the relative error (18) becomes: 4.1.Choice of the number of modes.As explained in Section 3.1, the number of modes L is generally selected through the cumulative energy threshold ε in equation ( 16) affecting the accuracy of the ROM approximation.In Figure 3 we plot the modes number as well as the relative error e (see equations ( 18) and ( 19)), which is time average for m 1 , m 2 and φ, against the value of ε = {90%, 95%, 99%, 99.9%}, for both the Newtonian and Casson's models.As expected, the relative error decreases at increasing of modes number, i.e. for higher values of ε.However for ε ≥ 99% the relative error reaches a plateau for most of the variables.On the contrary, the number of modes increases.Such increase is particularly pronounced for m 1 , m 2 and φ (and appears more evident in the Newtonian case, although the error level exhibited by the two models is comparable).This could be due to the fact that TAWSS and OSI derive from an average in time over the cardiac cycle whilst the other indices are time-dependent and therefore characterized by a richer modal content.We fix the value of ε to 99% for the following tests, as it guarantees a good accuracy and at the same time it allows to monitor the computational cost.

ROM solutions.
Once set the cumulative energy threshold ε = 99%, we proceed with the ROM simulations.In Figure 4 we show the variation in time of the relative error e(t, µ) defined in equation ( 18) for m 1 , m 2 and φ associated with the parameters values in K N test and K C test , together with their mean values.Overall, we observe that the relative errors vary between 9% and 21% demonstrating a clinical relevance.More specifically, in the Newtonian case (left plots), the relative errors are quite similar to each other.They, along with the associated mean value, do not exhibit large oscillations during the cardiac cycle.For the Casson's model (right plots) we observe that the error increases immediately after the opening of the MV (happening around t = 0.4, see Figure 1), so unlike the Newtonian case the error is affected by wide oscillations.We also note the increase of the error during the first few time steps which might be due to the transient nature of the flow, as also noted in other configurations including academic benchmarks [33].
In  Newtonian (Casson) case.For each variable, the FOM simulations (top plots) are very similar to the ROM ones (bottom plots) as expected by the accuracy analysis based on the relative error carried out above.Furthermore, by looking at the patterns obtained, some speculations of clinical interest can be made.Both models reveal an higher momentum m 1 in the terminal region of the LAA.Also the washout φ and m 2 are characterized by a larger magnitude on the tip of the appendage.Therefore we can argue that a higher residence time at the end of the LAA is shown.Note that, although in the Newtonian case a slightly greater variation between the base of the appendage and its tip is shown for all the variables under consideration (compare Figures 5 and 6, Figures 7 and 8, Figures 9 and 10), it would seem that for this study case the introduction of Casson's model is not justified.Finally, in agreement with the trend of the relative errors, also the qualitative comparisons of the Casson's case reveal a worse reconstruction of the reduced solution at t = 0.For what concerns time independent indices, in Figure 11 and 12 we see that the TAWSS distribution is very well reconstructed by ROM.As expected, it shows greater values where the blood age    is lower and the computations provided by Newton and Casson's models are very similar.However, a more pronounced difference between the two rheologies can be found in the OSI distribution.However, deeper investigations should be conducted, also in direct contact with medical doctors, to establish whether the differences shown by Netwon and Casson's models may be substantial or negligible.In Table 1 we report the computational time taken by the FOM and the ROM.In this case we refer to the entire domain (and not only to the LAA region) in order to provide a fair comparison.All the simulations have been run on the SISSA HPC cluster Ulysses (200 TFLOPS, 2TB RAM, 7000 cores), the FOM ones in parallel using 16 processors while the ROM ones by using one processor only.Each FOM simulation takes roughly 3.75 hours in terms of wall time, or 60 hours in terms of total CPU time, while the online phase only needs a few seconds.Thus the speed-up is of the order of 10 5 and the ROM is able to practically work in a real-time way.In Table 1 we also report the estimation of the time required for the computation of the POD basis and the RBF interpolation for sake of completeness.Such results are very promising and could push toward the transfer of ROM techniques in hospitals and surgery rooms by means of the development of user-friendly digital platforms to be accessed with portable devices, such as a smartphone or a tablet.In this context, we have designed the ATLAS project that allows computations to be run from standard web browsers.Further details could be found in [34].

Conclusions and perspectives
A data-driven ROM based on POD-RBF technique is adopted in this work for the analysis of the blood flow in a patient-specific domain of LA when AF occurs.Such approach extracts a reduced basis space from a proper set of high fidelity solutions via POD and adopts RBF to compute the map between parameter space and reduced coefficients.The Newtonian and the Casson's models for the rheology of the blood are employed and compared.We consider a physical parametric framework involving the scaling factor of the cardiac output (both for Newtonian and Casson's case), the plasma viscosity and the hematocrit (for the Casson's case).
After an expensive offline phase, the POD-RBF approach demonstrated to be able to provide clinically relevant blood flow predictions for the problem at hand at a considerably lower computational cost.This indicates that in perspective such computational tool could be used in hospitals and surgery rooms to support the medical doctors.From a clinical point of view, the results revealed that, unless some differences between the Newtonian model and the Casson's one, the distribution of the mean blood age is higher on the tip of LAA.
In perspective, an improvement of considerable relevance could be the exploration of the performance of the POD-RBF approach in a geometrical parametric setting by extending to our case study what carried out in [61].We are also going to move towards non-linear ROM methods (see, e.g., [28,49,52]) in order to improve the accuracy as well as the efficiency of our approach.

Acknowledgements
José Sierra-Pallares wants to acknowledge "Movilidad de Investigadores e Investigadoras UVa -Banco Santander 2022" for funding his stay at SISSA Trieste and project number DPI2017-83911-R from Spanish Minitry of Science, Innovation and Universities.Jorge Dueñas-Pamplona wants to acknowledge the "Programa Propio -Universidad Politécnica de Madrid", and the "Ayuda Primeros Proyectos de Investigación ETSII-UPM".We also thank the "Programa de Excelencia para el Profesorado Universitario de la Comunidad de Madrid" for its financial support and the "CeSViMa UPM project" for its computational resources.The acknowledgments are addressed also to the support provided by European Union Funding for Research and Innovation -Horizon 2020 Program -in the framework of European Research Council Executive Agency: H2020 ERC CoG 2015 AROMA-CFD project 681447 "Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics" P.I.Professor Gianluigi Rozza.This work was also supported by the "Gruppo Nazionale per il Calcolo Scientifico" (GNCS -INdAM) and the European Union Funding for Research and Innovation -Horizon Europe Program -in the framework of European Research Council Executive Agency: ERC POC 2022 ARGOS project 101069319 "Advanced Reduced order modellinG: Online computational web server for complex parametric Systems" P.I.Professor Gianluigi Rozza.We also thank the PRIN NA FROM-PDEs project.
Appendix A. Supporting materials A.1.Model validation.Model validation is a critical step in the development of computational models and simulations, particularly in engineering and scientific research.Model validation involves comparing the model's predictions to experimental or observational data, verifying that the model accurately captures the behavior of the system under study.In this work, the FOM is validated against the experimental benchmark reported by [24].In that work, an in vitro model mimicking a real left atrium geometry working with a blood mimicking fluid is measured by Particle Image Velocimetry (PIV).Figure 13 shows the experimental geometry and the measurement plane employed.Figure 14 shows the experimental measurements of the mean velocity in the in vitro model (top of the figure) along with boundary conditions (bottom of the figure) for mass flow rates and the predictions of our FOM implementation.Mean velocity is well predicted by the FOM, which is in line with the results provided in the benchmark paper [24].In fact, current OpenFOAM implementation gives almost identical results to those obtained by the authors in that reference.A.2. Mesh convergence analysis.Mesh convergence analysis is a crucial step in verifying a CFD code, for ensuring that a CFD simulation is properly resolved and that the numerical results are reliable.The Grid Convergence Index (GCI) [57] is used here since it is a widely accepted method for assessing the mesh convergence of a CFD simulation.The GCI provides an estimate of the order of accuracy and the uncertainty associated with the numerical solution obtained from a given mesh.The GCI is calculated by comparing the solutions obtained from two or more grids with different mesh sizes.
Table 2 shows the result of a GCI study on four different grids for the patient-specific geometry at hand (see Figure 1), taking a point located at the geometric center of the ostium as sample point.Mean velocity is evaluated in such grid point for four different grid spacing, assuming in all cases a time step of 10 −3 s and running the simulation for 5 complete cycles.Results show the solution on Grid 2 is mesh independent, therefore this grid size is used for all the simulations in this study.
Table 2. Grid convergence study over 4 grids.v represents the mean velocity at a point located in the ostium of the LAA and vextrapolated its extrapolated value.N h is the number of grid elements, r the refinement ration between two successive grids.GCI is the grid convergence index in percent and its asymptotic value is provided by GCI asymptotic , where a value close to unity indicates a grid independent solution.The order achieved in the simulation is given by p [57]

Figure 1 .
Figure 1.Left: LA patient-specific geometry; the green surfaces are the PV inlets while the red region is the MV outlet.Right: boundary conditions for MV flow for different scaling factors ranging in [0.5, 1.5] (the dashed line corresponds to 1).

Figure 2 .
Figure 2. Sketch of the mesh.The red portion corresponds to LAA.

Figure 3 .
Figure 3. Variation of the number of modes and of the relative error with respect to the cumulative energy threshold ε for the Newtonian (first column) and non Newtonian (second column) case for all the variables involved.

Figures 5 -
12 we show some qualitative comparisons between the FOM and ROM solutions for the Newtonian and Casson's cases.The plots refer to the test parameters values µ N 3 and µ C 1 .Figures 5, 7 and 9(6, 8 and 10)  show the results obtained for the variables m 1 , m 2 and φ, respectively, for the Casson's case: φ

Figure 4 .
Figure 4. Variation in time of the relative error corresponding to the test parameters values for the Newtonian (first column) and non Newtonian (second column) case for m 1 , m 2 and φ.

Figure 5 .
Figure 5. Newtonian case: qualitative comparison between FOM (top) and ROM (bottom) solutions at different times for m 1 .

Figure 6 .
Figure 6.Casson's case: qualitative comparison between FOM (top) and ROM (bottom) solutions, different times for m 1 .

Figure 9 .
Figure 9. Newtonian case: qualitative comparison between FOM (top) and ROM (bottom) solutions at different times for φ.

Figure 10 .
Figure 10.Casson's case: qualitative comparison between FOM (top) and ROM (bottom) solutions at different times for φ.

Figure 11 .
Figure 11.Newtonian case: qualitative comparison between FOM and ROM solutions for TAWSS (panels A and B) and OSI (panels D and E).

Figure 12 .
Figure 12.Casson's case: qualitative comparison between FOM and ROM solutions for TAWSS (panels A and B) and OSI (panels D and E).

Figure 13 .
Figure 13.Description of validation experiment taken from [24].A: Isometric view of the in vitro LA with inflows Q P V, i , mitral valve outflow Q M V and LAA volume highlighted; B: Front view of the in vitro LA; C: Detailed view of the measurement section and plane at the in vitro LAA.

Table 1 .
CPU time taken by the offline and the online phases related to the whole computational domain. .