{'fraction_non_alphanumeric': 0.03218061159835929, 'fraction_numerical': 0.012388701770767332, 'mean_word_length': 4.428093040094126, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A scalable spectral Stokes solver for simulation of time-periodic flows in complex geometries

Simulation of unsteady creeping flows in complex geometries has traditionally required the use of a time-stepping procedure, which is typically costly and unscalable. To reduce the cost and allow for computations at much larger scales, we propose an alternative approach that is formulated based on the unsteady Stokes equation expressed in the time-spectral domain. This transformation results in a boundary value problem with an imaginary source term proportional to the computed mode that is discretized and solved in a complex-valued finite element solver using Bubnov-Galerkin formulation. This transformed spatio-spectral formulation presents several advantages over the traditional spatio-temporal techniques. Firstly, for cases with boundary conditions varying smoothly in time, it provides a significant saving in computational cost as it can resolve time-variation of the solution using a few modes rather than thousands of time steps. Secondly, in contrast to the traditional time integration scheme with a finite order of accuracy, this method exhibits a super convergence behavior versus the number of computed modes. Thirdly, in contrast to the stabilized finite element methods for fluid, no stabilization term is employed in our formulation, producing a solution that is consistent and more accurate. Fourthly, the proposed approach is embarrassingly parallelizable owing to the independence of the solution modes, thus enabling scalable calculations at a much larger number of processors. The comparison of the proposed technique against a standard stabilized finite element solver is performed using two- and three-dimensional canonical and complex geometries. The results show that the proposed method can produce more accurate results at 1% to 11% of the cost of the standard technique for the studied cases.


Introduction
Fast and cost-effective simulation of flow in complex geometries is highly desirable, at least for two sets of applications. The first is the time-critical problems where the results must be obtained in a given time window, thus imposing an upper bound on the time-to-solution. The second are those with a bound on the computational cost, where the number of CPU-hours spent on the simulation must be limited. While these two restrictions on time and cost are identical for sequential simulations, they may well diverge for parallel implementations depending on the scalability of the application.
A prominent area where all the above conditions and restrictions apply is cardiovascular modeling. The flow is unsteady and time-periodic, the geometry is often complex, the cost must remain low, and the available time for simulation is often limited. A specific example within this realm is patients with congenital heart disease who must undergo an operation right after birth [1]. Utilizing cardiovascular modeling as a part of patient-specific surgical planning will require the simulation to be completed between diagnosis (birth) and the time of operation (a few days after birth at most) [2]. Moreover, performing shape optimization to improve surgical anatomy relies on having a highly efficient computational framework since a formal optimization study typically requires tens to hundreds of simulations [3]. A similar constraint exists on the efficiency of the underlying computational method if such computations were to be performed outside of the academic setting and at an industrial scale, where access to high-performance computing resources for a large number of patients remains limited due to the cost considerations [4].
While there exists a broad range of numerical methods for simulating flow in complex geometries, here we consider stabilized finite element techniques to establish the relative accuracy and cost-efficiency of our method since these methods are widely adopted for cardiovascular blood flow simulation [5] and closest to the method proposed in this study. Within the realm of stabilized finite element methods (FEM), we particularly focus on the residual-based variational multiscale method (RBVMS), which is constructed on top of streamwise upwind Petrov-Galerkin (SUPG) technique [6,7,8].
As is the case with the majority of the time-dependent partial differential equations solvers, the application of the RBVMS technique to a time-periodic problem (e.g., blood flow in the circulatory system) requires discretization of the underlying governing equations in time. Thus, the cost of these methods is inherently dependent on the time step size, increasing as time step size decreases. While the maximum allowable time step size is limited by the stability considerations for explicit time integration schemes, it is often more relaxed for implicit schemes, where it is determined solely based on the desired accuracy [9,10]. The larger time step size for implicit schemes does not necessarily imply their lower cost, as the fewer time steps can be offset by the larger number of Newton-Raphson iteration per time step. Furthermore, the high cost associated with the time integration is exacerbated by the transient nature of the solution that requires simulation of many cycles to achieve cycle-to-cycle convergence. In addition to these setbacks, one well-known issue associated with the RBVMS formulation is its lack of consistency with regard to time step size that, for instance, leads to the time step size dependency of the final solution obtained from a steady-state problem.
A remedy for the aforementioned issue is to transfer the temporal problem to a frequency domain using Fourier transformation and subsequently solve the transformed boundary value problem mode-by-mode. Since the resulting spectral method bypasses the time integration, its accuracy and stability are not affected by the time integration scheme. We should note that the number of modes that must be incorporated into such calculation, at least in theory, must be equal to the number of time steps for all time scales in the problem to be captured. However, for cardiovascular flows, the majority of those modes can be neglected with minimal error, as we will show later in this article. In practice, a few modes (less than ten) suffices for sufficiently smooth and time-periodic boundary conditions, thus offering a significant reduction in overall cost in comparison to a traditional solver that requires thousands of time steps for simulation of a single cardiac cycle [11]. On top of this saving, the spectral technique can save on the computation required for cycleto-cycle convergence as it is based on a boundary value problem that does not include a transient solution. Additionally, for the linear Stokes equation that is the focus of this study, all the solution modes are orthogonal and independent, which present two major advantages. Firstly, capturing a new mode in the solution (either due to a modification in the boundary condition or stricter tolerance on the solution) can be done by solving one additional boundary value problem as there is no need to recompute previous modes. The minimal cost can be contrasted against the costly alternative of a traditional time integration scheme in which the entire computation must be repeated with a new time step size if higher accuracy is desired. Secondly, the spectral technique can be implemented as an embarrassingly parallelizable scheme with each group of processors computing the solution associated with only one of the many modes. Therefore, by formulating the original temporal problem as a set of orthogonal boundary value problems, one can significantly reduce the total computational cost and, at the same time, scale computations to a much larger number of processors.
The enumerated advantages of formulating a problem in the spectral domain have attracted researchers to investigate such schemes in the past. Perhaps the most relevant scheme to the present study that has come out of such efforts is the development of the Time Spectral Method, which has been primarily employed by the aerospace community for simulation of compressible flows over bluff bodies, such as blades in a compressor rotor and a pitching airfoil [12,13,14,15]. In this paper, we propose a new spatio-spectral technique named Complex-valued Stokes Solver (SCVS) that can be differentiated from the Time Spectral Method in a number of ways. First and foremost, the aforementioned Time Spectral Methods are formulated for the Euler's equations that are nonlinear and the SCVS for the Stokes' equations that are linear. This nonlinearity that creates a coupling between all modes is dealt with by establishing a relationship between solutions at all time steps and developing a pseudo-time-advancement technique that solves for all the modes (i.e., time steps) at once. In contrast, the SCVS exploits the linearity of the governing equations and is directly formulated in the spectral domain. As a result, the SCVS is implemented using complex arithmetic in contrast to the Time Spectral Method that is based on pure real arithmetic. Thirdly, we specifically introduce the SCVS for the simulation of creeping blood flow in the cardiovascular system. What distinguishes cardiovascular systems from prior applications is the geometrical complexity, as such models typically involve multiple inlets and outlets that may be coupled to lower-order models as boundary conditions [16]. Fourthly, in contrast to the aforementioned studies that employ finite volume technique for spatial discretization, we experiment with the finite element method in the SCVS as a widely adopted technique for cardiovascular modeling [17,18]. Lastly, in this study, we discuss in detail the convergence, accuracy, and computational efficiency of the SCVS as a function of element size, number of modeled modes, and the underlying linear solver tolerance.
The paper is organized as follows: we introduce the governing equations, namely Stokes equation with a complexvalued source term, in the spectral domain in Section 2.1. Second, we present the weak and discrete formulations of those equations based on finite element formulation in Sections 2.2 and 2.3, respectively. Then, we introduce the SCVS algorithm and elaborate on how to estimate the error in Sections 2.4 and 2.5. Next, we move on to Section 3 where we evaluate the introduced technique using four cases, where we compare the solutions obtained from the SCVS and a standard RBVMS solver. The convergence and accuracy of the SCVS are established for a canonical case where an analytical solution is available. The computational efficiency of the SCVS is demonstrated and contrasted against that of the RBVMS in Section 3.6. Finally, we have a look at possible research directions before drawing conclusions at the end in Sections 4 and 5 respectively.

Complex-valued Stokes solver formulation
In what follows, we first briefly describe the complex-valued Stokes equation as the governing equation for low Reynolds number flows, discuss its discrete formulation using FEM, and then present the SCVS algorithm in detail.

Governing equations
Consider creeping fluid flow in a domain Ω with boundary Γ = ∂Ω that is governed by the incompressible unsteady Stokes equations as where u(x, t) refers to the velocity of the fluid at location x and time t, p(x, t) is the pressure, ρ is the density, n is the boundary Γ outward normal vector, and µ is the viscosity. g and h refer to the imposed velocity and traction on the Dirichlet Γ g and Neumann Γ h boundaries where Γ = Γ g ∪ Γ h . Additionally, suppose that these boundary conditions are time-periodic with a period of T . To represent Eq. (1) in the spectral domain, we make use of Here the frequency ω i is defined in terms of period T as ω i = 2πi/T for i = 0, 1, . . . , N m with N m denoting the largest computed mode. For time-periodic flows, such as those encountered in cardiovascular modeling, T is the cardiac cycle, which is the time it takes for boundary conditions to repeat themselves. Similar to the state variables, Dirichlet and Neumann boundary conditions can be expressed in the spectral domain as g(x, t) = igi (x)eˆj tω i t and h(x, t) = ihi (x)eˆj tω i t , respectively. With these definitions, Eq. (1) and continuity equation can be written aŝ These equations resemble those of the steady Stokes equation with a nonzero source term. The only difference is that the first (source) term in Eq. (3) is complex-valued unless ω i = 0, for which the steady Stokes equations in the real domain are recovered.

Weak formulation
Since the solution to Eq. (3) at ω i and ω j for i j are independent, we only need to provide a solution procedure at a single frequency ω i . Thus, for the sake of notation brevity, we drop subscript i and denote our formulation in terms of variable ω below.
The weak form of Eq. (3) is stated as follows. Given the frequency ω, findũ ∈ S andp ∈ P such that for any w ∈ W and q ∈ Q holds. In Eq. (4), w and q are test functions for velocity and pressure, respectively, and are velocity and pressure solution and test function spaces. L 2 denotes the space of scalar-valued functions that are square-integrable on Ω, and (H 1 ) n sd denotes the space of vector-valued functions with square-integrable derivatives on Ω.

Discrete formulation
Denoting the finite-dimensional subspace of S, W, P, and Q by S h , W h , P h , and Q h , respectively, we attempt to solve the discrete form of Eq. (4) using Galerkin's approximation. Namely, we seekũ h ∈ S h andp h ∈ P h such that holds. In writing Eq. (6), we assumed Ω h = Ω, i.e., the computational domain after discretization remains unchanged. If the two differ, then one must perform the integrals over Ω h rather than Ω when computing B G and F G in Eq. (6). Galerkin's formulation of incompressible flow has a saddle-point nature, which produces a singular system if one were to adopt similar shape functions for velocity and pressure. Various techniques have been proposed to overcome this issue, ranging from mixed-element [19] and penalty techniques [20] to stabilized finite element methods [21]. In the present study, we adopt the mixed finite element method, which allows us to reduce the number of variables that influence the accuracy of the SCVS formulation, thus simplifying the measurement of its numerical properties. More specifically, we employ linear and quadratic shape functions for pressure and velocity, respectively, to satisfy inf-sup condition (also known as LBB condition) [22,23,24]. In 2D and 3D, we use linear triangular and tetrahedral elements with 3 and 4 nodal points, respectively, for pressure. For velocity, we use quadratic shape functions that produce 6 and 10 nodes in each of the aforementioned elements in 2D and 3D, respectively. Denoting these linear and quadratic shape functions at global node A by N A (x) and M A (x) , respectively, the test functions and unknowns are interpolated in space using where η, η g , andη refers to the velocity nodes, velocity nodes on the Dirichlet boundaries, and pressure nodes, respectively. U A , P A , W A , Q A in Eq. (7) are the velocity and pressure unknowns and their respective test functions. G A is the prescribed velocity defined on the Dirichlet boundaries after discretization such that where Π hg is an operator that projectsg to the finite-dimensional discrete space. Substituting for the variables appearing in Eq. (6) using Eq. (7) while ensuring the results holds for any W A and Q A produces the following system of linear equations where and K and D matrices and B vector are computed as The solution to this linear system is obtained using a brute-force approach. Namely, we treat the entire linear system as a single sparse matrix and solve it iteratively using the Generalized minimal residual method (GMRES) method [25] with a Jacobi preconditioner. All these methods are implemented in our in-house linear solver that employs an optimized data structure for scalable distributed parallel processing [26]. The GMRES algorithm used for this purpose is slightly modified from its real counterpart owing to K being complex-valued. The same linear solver and preconditioner are also used for the RBVMS to allow for a one-to-one comparison of computational cost between the two. This choice has been made due to the underlying structure of the tangent matrix, which is not symmetric for the RBVMS, as well as its relative simplicity and wider use by the finite element community. Provided that the left-hand side matrix in the case of SCVS is symmetric and has a block structure, one may adopt a linear solution and preconditioning strategy based on Schur complement to reduce the overall cost of solving this linear system in the future [27,28,29,30]. The use of multigrid [31] is also expected to significantly reduce the cost of linear solve, particularly on finer grids.

SCVS algorithm
In this section, we discuss the implementation of the SCVS algorithm. Overall, the procedure involves representing the boundary conditions in the spectral domain, then solving for flow at each mode, and reconstructing the solution in the real domain.
1. Take the Fourier transformation of the boundary conditions g(x, t) and h(x, t) as 2. Select the number of modes to truncate the above series based on the smoothness of the boundary conditions. Note that the number of required modes can be determined directly from the frequency content of the prescribed boundary conditions due to the linearity of the Stokes equations.
More specifically, one can select a tolerance ǫ m to obtain a priori estimate of the number of required modes N m + 1 such that eˆj ω i t are the truncated Neumann and Dirichlet boundary conditions using N m +1 Fourier modes and h 2 Even though we are selecting the first N m modes to truncate the series in this study, in general one may select a non-consecutive set of modes solely based on their amplitude to accelerate the rate of convergence. 3. Construct and solve the linear system in Eq. (9) for ω 0 , ω 1 , . . . , and ω N m given the boundary conditions computed from the previous step, i.e.,g i (x) andh i (x). 4. Given the solution U A and P A , reconstruct the solution in the spectral domain using Eq. (7) and then in time as The number of computed modes N m + 1, which was selected based on a priori error estimate, can be refined following the computation of each solution mode. This way, one can compute the difference between the solutions computed from Eq. (14) using N m − 1 and N m modes and decide whether to continue computing the next mode if that difference is larger than a prespecified tolerance.
This flexibility presents a significant advantage over the conventional spatio-temporal techniques where capturing an additional mode requires repeating the entire simulation.
Later in Section 3.5, we will show the rate at which a priori (Eq. (13)) and a posteriori (Eq. (15)) errors drop as N m increases are the same, implying that one can simply rely on Eq. (13) to select N m . In our experience, N m ≤ 7 is typically sufficient for cardiovascular applications where the boundary conditions are smooth and periodic in time.

Error estimation
To ensure consistency and accuracy of the proposed formulation, we will consider a set of canonical cases where an analytical solution to the governing equations is available. The overall error of a conventional scheme stems from spatial discretization error, the linear solver error, time integration error, and in the case of nonlinear equations, Newton-Raphson iteration error. If we consider the solution at a single mode, i.e., boundary conditions oscillating at a single frequency, then SCVS is only prone to the first two types of error enumerated above (discretization and linear solver). Thus, in what follows, we establish the overall accuracy of the SCVS algorithm as the mesh size or linear solver tolerance is varied. Also, for scenarios in which the boundary conditions are generic functions of time, we will study the convergence of the SCVS solution as a function of N m .
The overall error e(t) is defined based on the relative difference between the computed velocity u h (x, t) and the reference velocity u(x, t) over the domain Ω as where the L 2 (Ω) norm is defined as u 2 L 2 (Ω) =´Ω u · udΩ. As discussed earlier, this total error can be decomposed to the error due to spatial discretization e H and the error due to the linear solver e L . Denoting the best approximation of the solution in the finite dimensional solution space S h by Π h u, L 2 -norm of the interpolation error for the quadratic shape functions employed here can be estimated as [32] e where C 1 is a constant and h is the diameter of the smallest element-bounding circle in 2D and sphere in 3D.
The error related to the linear solver e L is produced by the approximate nature of the iterative solution procedure employed in solving Eq. (9). This error can be related to the linear solver tolerance ǫ L as in which X is the approximate solution obtained from the iterative linear solver. In Eq. (17), the constant C 2 depends on the inverse of tangent matrix A −1 from Eq. (10) and thus is independent of the boundary conditions. How the overall error e changes with e H and e L for a pipe flow will be investigated in Section 3.4.

Results
The solution procedure and four case studies selected for evaluating the SCVS algorithm are discussed in this section. To establish the relative accuracy and cost of the SCVS, we will compare it against a standard RBVMS solver. The details pertaining to the RBVMS algorithm are included in Appendix A. The four test cases considered in this study are (1) a 2D channel, (2) a 2D diverging nozzle, (3) a 3D pipe, (4) a patient-specific model of Glenn procedure with an anastomosis [33]. Since the analytical solution is available for the first and third cases, these two cases will be employed to validate our implementation and evaluate its numerical characteristics. While the boundary conditions are periodic and contain a single frequency for those two cases, they are generic functions of time for the case two and four, thus allowing us to study the convergence of solution with regard to N m and its relationship to e M . At the end of this section, the performance of the SCVS and the RBVMS are compared in terms of total CPU time and wall-clock time to demonstrate the performance of the SCVS.
The SCVS uses mixed quadratic-linear shape functions for velocity and pressure, respectively, whereas the RB-VMS uses linear shape functions for both. This incompatibility raises a concern that comparing the two is similar to comparing apples to oranges. To address this concern, we have implemented a Stokes solver with mixed shape functions (called MSS) that is identical to the SCVS except for the traditional discretization of the acceleration term in the time domain. The formulation of the MSS along with its performance metrics relative to the SCVS and the RBVMS are included in Appendix B, where we also explain why the SCVS is solely compared against the RBVMS for the remainder of this article.

Solution procedure
The 2D and 3D models are initially discretized using triangular and tetrahedral elements. For this purpose, we use a combination of Tetgen [34] and Simvascular [35] software. An in-house script was developed to generate quadratic elements from the linear triangular and tetrahedral elements through a node insertion process. For boundary nodes, special care was taken to ensure that inserted nodes are properly projected to be located on the curved boundaries. Since the SCVS simulations are performed using quadratic and the RBVMS using linear elements, we reduce the number of elements for the quadratic meshes such that the total number of degrees of freedom is roughly the same as the linear meshes, thus allowing for a one-to-one comparison between the two. A zero initial condition is used for all the RBVMS simulations. The generalized-α method is used for time integration of the RBVMS, which is an implicit second-order method [36]. As detailed in Appendix A, despite some similarities with the operator splitting techniques [37], velocity and pressure are solved together during each linear solve. The time step size for the RBVMS is selected to sufficiently resolve the time variation of boundary conditions. More specifically, we employ 2,000 time steps for cases with time-varying boundary conditions. For the steady simulations, the time step size is selected such that the diffusive Courant-Friedrichs-Lewy (CFL) number ∆tν/h 2 = 1. The time integration is continued for steady cases until the relative residual falls below 10 −6 . A total of five cycles are simulated to achieve cycle-to-cycle convergence for the RBVMS simulations. All the results shown below correspond to the last simulated cycle.
For all cases, the Reynolds number is Re < 10 −3 to ensure nonlinear effects are negligible in the case of RBVMS solver. The Reynolds number is defined based on the maximum flow rate through boundaries, the area of the boundary with the maximum flow rate, and the kinematic viscosity ν = µ/ρ. We only perform a single Newton-Raphson iteration in the RBVMS simulations, given that the governing equations are linear.
The GMRES method is adopted for solving the resulting linear system for both the SCVS and RBVMS using the same tolerance ǫ L = 10 −6 . We used our in-house linear solver for this purpose [30,29,26]. As discussed earlier, the linear system in the case of the SCVS is complex-valued, requiring us to make a necessary adjustment to the realvalued version of our GMRES solver to accommodate for a complex-valued linear system. Both the RBVMS and the SCVS solvers are also implemented in our in-house finite element solver, which is written in object-oriented Fortran and parallelized using the MPI library.

2D channel flow
A 2D channel with a length-to-height aspect ratio of 5 is considered as our first case study (Fig. 1). The geometry is discretized using 3,762 linear triangle elements for the RBVMS solver and 882 quadratic elements for the SCVS solver, producing 2,000 and 1,881 nodes and 6,000 and 4,262 degrees of freedom for two cases, respectively. The non-slip boundary condition is imposed for the top and bottom walls. Either a nonzero constant or a cosinusoidal Neumann boundary condition is imposed on the inlet, and zero traction is imposed at the outlet. These boundary conditions produce a solution that corresponds to the fully developed condition observed in an infinitely long channel (i.e., pressure gradient does not vary as a function of x). The amplitude of the wave, which otherwise is irrelevant due to the linearity of the Stokes equations, is selected such that the Reynolds number is sufficiently small for the RBVMS computations. The oscillation frequency ω is varied to simulate a wide range of conditions. More specifically, the Womersley number W = ωH 2 /ν with H denoting the half-channel height varies from 0 to 20π for the results shown below. An analytical solution is available for the oscillatory flow in the 2D channel, expressed in the spectral domain as a function of local height y and ω according to [38] where Λ 2 =ĵW,ũ x is the streamwise velocity in the spectral domain, L is the channel length, andh x is the traction amplitude at ω imposed at the inlet acting in the streamwise direction. This solution can also be expressed in time as The SCVS and the RBVMS simulation results are compared against the analytical solution for W = 2π, W = 10π, and W = 20π at two time points t = T/4 and t = T/2 in Fig. 2. In general, the error increases as the flow becomes more oscillatory (at higher W). This larger error can be attributed to the velocity profile developing sharper gradients near the walls. At a similar W and time point, the SCVS provides more accurate predictions in comparison to the RBVMS. This higher accuracy is despite the fact that a larger number of degrees of freedom were employed in the RBVMS simulations. The primary reason for this improved accuracy is the use of quadratic shape function for the SCVS. The time integration scheme and the stabilization terms also contribute to the larger errors in the RBVMS results.
A more condensed version of these results is provided in Table 1, where the relative error at time point T/4 and T/2 are computed using Eq. (15) for the RBVMS and the SCVS. In the steady case, since the analytical solution is a parabola and can be exactly represented by the quadratic shape functions, the SCVS error of e = 1.34 × 10 −5 is solely due to the linear solver. Repeating this computation with a smaller ǫ L confirms that e for the SCVS can be arbitrarily lowered without refining the grid. The steady solution for the RBVMS, on the other hand, shows a relatively larger error of 1.93% that is grid-dependent given that it utilizes linear shape functions. The same improved accuracy is observed at higher modes, where for instance, the RBVMS produces a solution with a 13.4% error at T/2 and W = 20π, whereas the SCVS error remains as low as 1.83%. In general, in comparison to the RBVMS method, the SCVS is an order of magnitude more accurate for a mesh with a similar number of nodes.   Fig. 1. The velocity is normalized using analytical u max . The results on the left and right columns are extracted at t = T/4 and t = T/2, respectively, and those on the first, second, and third row correspond to W = 2π, 10π, and 20π, respectively.

2D diverging nozzle
In the second example, we consider a 2D diverging nozzle with an expansion ratio of 2 (Fig. 3). This geometry is discretized using 2,400 linear triangle elements for the RBVMS solver and 600 quadratic elements for the SCVS solver, producing 1,281 nodes in both cases with 5,124 and 4,184 degrees of freedom, respectively. A time-periodic Neumann boundary condition is imposed on the inlet, and zero traction is imposed on the outlet (Fig. 3).  Even though the inlet Neumann boundary contains a wide range of frequencies, only the first few modes have a large enough amplitude to have a significant effect on the results. As shown in Fig. 4, h i follows a power-law trend versus ω i , experiencing a sharp drop that is approximately proportional ω −2 i . Such a sharp drop in the amplitude at higher frequencies justifies the truncation of the series at a relatively low wavenumber. Thus, to construct the SCVS solution, we select N m = 5 and solve for ω 1 , . . . , ω 5 while neglecting ω 0 given the time-average of the inlet boundary condition is zero. The Womersley number, i.e., W = ωH 2 /ν with H being the half-inlet height, for these simulated modes ranges from 0 to 10π. To show the convergence of SCVS solution with N m , we reconstructed solution at t = T using one (N m = 1), three (N m = 3), and five (N m = 5) modes. The results of these computations normalized by u * = |h x (ω = 1)|H/µ for N m = 1, 3, and 5 are shown in Fig. 5 along with the results obtained from the RBVMS. This figure shows the qualitative convergence of the SCVS at N m = 3. For a more quantitative comparison of two formulations, the dimensionless flow rate Q/Q * , with Q * = Hu * being the characteristic flow rate, is computed and shown in Fig. 6 for the RBVMS and the SVC with N m = 1, 3, and 5.
Taking the flow computed from the SCVS solution with N m = 10 as the reference, the relative error in the prediction of the RBVMS is 1.0%, whereas the relative error for the SCVS at N m = 1, 3, and 5 is 25%, 5.0%, and 1.5%, respectively. The fast rate of convergence of the SCVS with respect to N m confirms our earlier argument that only a few modes are needed to obtain a converged solution using the SCVS algorithm when boundary conditions vary smoothly in time.

3D Pipe flow
For the third case study, we consider an oscillatory laminar pipe flow. A pipe with an aspect ratio of L/R = 15 is considered with a cosinusoidal inlet and zero outlet Neumann boundary condition (Fig. 7). Similar to the 2D channel flow case in Section 3.2, these boundary conditions produce a solution that corresponds to the fully developed condition observed in an infinitely long pipe. The oscillation period is varied to simulate flow at eleven Womersley numbers W = ωR 2 /ν = 0, 8π, . . . , 80π with R denoting the radius of the pipe. As we discussed earlier, the domain is discretized in space using linear and mixed quadratic-linear tetrahedral elements for the RBVMS and the SCVS simulations, respectively. A wide range of meshes is constructed to investigate the convergence of both solvers ( Table 2). The number of elements in corresponding meshes is selected such that the number of nodes is roughly the same for both solvers. Overall, the element size normalized by the pipe radius varies from 0.0340 to 0.228 among simulated cases.  Similar to the first case considered above, an analytical solution is available for this case that permits us to establish the accuracy of each solver. For an oscillatory flow in a pipe, the solution in the spectral domain is expressed as [39] where Λ 2 = −ĵW and J 0 is the zero order Bessel function of the first kind. This solution can be expressed in time using u x (r, t) = real ũ x (r, ω)eˆj ωt .
For a qualitative evaluation of the SCVS and the RBVMS solutions, simulated u x (r, T/4) and u x (r, T/2) along with the analytical prediction of Eq. (21) at W = 8π, 40π, and 80π are shown in Fig. 8. These results are obtained using the finest grids, namely QM4 and LM4 in Table 2. Similar to what we observed earlier in Section 3.2 for the 2D channel flow, the SCVS provides a better approximation, particularly at larger Womersely numbers or at t = T/2, in which the solution exhibits sharper gradients. The error in the steady state solution (W = 0) obtained from the RBVMS and SCVS solver using the finest grids is e = 9.0 × 10 −3 and e = 1.4 × 10 −4 , respectively. The superior accuracy of the SCVS, which is achieved despite using fewer degrees of freedom, is primarily a result of utilizing a higher-order shape function. The numerical integration and stabilization terms are the secondary contributors to the larger error in the case of the RBVMS.   Fig. 7) predicted using the SCVS (solid black), the RBVMS (dashed red), and analytical solution (circles). The results on the left and right columns are extracted at t = T/4 and t = T/2, respectively, and those on the first, second, and third row correspond to W = 8π, 40π, and 80π, respectively. Computations are performed on QM4 and LM4 grids for the SCVS and the RBVMS, respectively. All the results are normalized using the maximum velocity from the analytical solution u max .
To study mesh convergence, the error e(T/2) defined in Eq. (15) is computed for the SCVS and the RBVMS as a function of h at W = 8π (Fig. 9). Third-order accuracy is observed for the SCVS, indicating that the measured error is dominated by the discretization error e H . This third order accuracy is in agreement with the estimate in Eq. (16), which also predicted e H ∝ h 3 . A similar relationship can be obtained for the RBVMS, which utilizes linear shape functions, as The second-order accuracy of RBVMS is also confirmed by the results shown in Fig. 9. To further analyze each method's accuracy for various flow conditions, we next study the overall error as a function of the Womersley number W. As we observed earlier in Figs. 2 and 8, the velocity profile develops sharper gradients near the wall at higher modes, thereby producing a larger error at higher W. That qualitative observation is confirmed quantitatively in Fig. 10 that shows a direct power-law relationship between e and W. This relationship indicates that resolving the velocity field at higher modes requires improved spatial resolution. It also shows that the quality of the solution degrades faster for quadratic elements than the linear elements as W increases. In what follows, we explain these observations in more detail and analytically predict the exponents that appear in the power-law relationship between e and W for the RBVMS and SCVS. Since the linear solver tolerance ǫ L is selected to be sufficiently small in this case and e L is negligible in comparison to e H , we can explain e ∝ W 3/2 for the SCVS and e ∝ W for the RBVMS by further analyzing the behavior of e H for each solver. To utilize Eqs. (16) and (22) for this purpose, we first need to establish how u L 2 , u H 2 , and u H 3 vary with W. Considering only cases in which W 0 and the fact that u x = u x (r, t), from Eq. (20) we can write where , andz 1 denotes complex conjugate of z 1 . Neglecting the dependence of these norms on the integrals of z 1 , z 2 , and z 3 , it is straightforward to show that indicating u L 2 , u H 2 , and u H 3 should vary proportional to W −1 , W 0 , and W 1/2 . Comparing these approximate predictions against reference quantities (obtained from 1D numerical integration of Eq. (23)) show that all exponents are over-predicted by 0.2 (Fig. 11). The ratio of these norms, however, is correctly predicted using Eq. (24) as indicating that the RBVMS and the SCVS error should grow proportional to W and W 3/2 at sufficiently small ǫ L owing to Eqs. (16) and (22), respectively. These analytical predictions are in agreement with our numerical results that were shown earlier in Fig. (10). Next, we investigate the effect of linear solver tolerance ǫ L on the overall accuracy of the SCVS, utilizing QM4 mesh and the solution corresponding to W = 8π evaluated at t = T/2. The left and right tail of the results shown in Fig. 12 indicate that the overall error e is independent of e L for sufficiently small e L and proportional to ǫ L when e L is sufficiently large. This observation confirms that (as expected) the overall error is dominated by the discretization error e H and that of the linear solver e L at small and large ǫ L , respectively. The transition between dominance of e H and e L occurs at approximately ǫ L = 10 −6 , which is a value specific to this case study and a function h among others. Note that this tolerance is the optimal tolerance with regard to the computational cost, as further decrease in ǫ L is not met with an improved overall solution accuracy. Lastly, note that the slope of 1 observed on the right tail of Fig. 12 is in agreement with Eq. (17), which predicts e L ∝ ǫ L . If we take the main objective of designing a computational method to be making a viable prediction at a specified accuracy with the lowest cost possible, then it is desirable to compare the SCVS and the RBVMS methods on an error-versus-cost plot. Thus, we have extracted the cost (in terms of CPU time) and error (in terms of e(T/2)) for all the simulations performed in this section using various grids and W. The results are shown in Fig. 13. Comparing the symbols with the same color, i.e., for a flow at a given W, the SCVS always provides more accurate predictions at a lower cost. Focusing on W = 80π cases (blue symbols), the SCVS yields a similar accuracy at the coarsest grid (solid triangle) to that of the RBVMS at the finest grid (hollow circle) while reducing the cost by three orders of magnitude. If one holds the computational cost relatively similar (solid blue circle versus hollow blue triangle, for instance), the SCVS offers over an order of magnitude improvement in accuracy relative to the RBVMS.  Table 2. The black, red, and blue colors correspond to W = 8π, 40π, and 80π, respectively. These results show that the SCVS, compared to the RBVMS, always provides a higher accuracy at a lower cost.

3D patient-specific Glenn
For the last test case, we will consider a complex patient-specific geometry acquired from a patient undergoing Glenn operation [33]. As shown in Fig. 14, the superior vena cava (SVC) is anastomosed to the left pulmonary artery (LPA) and right pulmonary artery (RPA) in this operation. The geometry is discretized using 988,747 linear tetrahedral elements for the RBVMS solver and 132,066 quadratic tetrahedral elements for the SCVS. This discretization results in 163,791 and 183,708 nodes and 655,164 and 574,401 degrees of freedom for the linear and quadratic meshes, respectively. The boundary conditions are selected based on physiologic data and a Windkessel model [3,40] for the RPA and LPA faces (Fig. 14). The wall is considered a non-slip boundary. In this example, twelve modes (N m = 11) are simulated in total for the SCVS. The case with N m = 11 is used as the reference to measure error for the remaining cases, given that the solution accuracy of the SCVS was superior to that of the RBVMS based on the cases discussed earlier. The Womersley number W = ωD 2 h /ν is based on the hydraulic diameter of the SVC (D h ), which is defined such that πD 2 h /4 is equal to the area of the SVC boundary. This definition, which is the square of what is normally defined in the literature as Womersley number, leads to W = 0, 146π, 293π, . . . , 1611π for simulated ω 0 , ω 1 , ω 2 , . . . , ω 11 . The results of these computations are normalized by u * = |h SVC (ω = 0)|D h /µ and Q * = πD 2 h u * /4. A closer examination of the solutions associated with each mode shows that their spatial distribution widely varies as W increases (Fig. 15). While the velocity peaks at the center of the vessels for the steady solution (Fig. 15-(a)), it develops two peaks closer to the walls at a larger W. This behavior resembles the trend that we observed earlier in the canonical geometries (e.g., cylinder in Section 3.5). This behavior is a result of a change in the relative importance of various terms that appear in the Stokes equation. While only the viscous and pressure terms are active at W = 0, three terms (acceleration, pressure, and viscous) balance each other for W > 0. As W increases, the relative importance of the acceleration term increases, leading to a solution that exhibits peaks near the walls rather than at the center of the geometry.
The SCVS solution in the time domain is reconstructed using Eq. (14) and qualitatively compared against that of the RBVMS (Fig. 16). This qualitative comparison shows that a few modes (in this case 6, including the steady solution) are sufficient for resolving the solution in this complex geometry. For a more quantitative comparison between the SCVS and the RBMVS solutions, the predicted flow through three branches is shown in Fig. 17. This figure confirms that the difference between the two solvers' solution reduces as more modes are employed. Using the case with N m = 11 as a reference, the relative errors of flux computed from the RBVMS are 0.65%, 0.38%, and 0.34% on the SVC, LPA, and RPA, respectively. The relative errors in the SCVS solution with N m = 7 are 0.14%, 0.19%, and 0.12% for the SVC, LPA, and RPA, respectively. Note that, due to the error in the RBVMS solution (as well as discretization error in the SCVS), the two solutions do not converge even at very large N m . Nevertheless, the SCVS does converge to a solution as N m → ∞. The rate of this convergence depends on the smoothness of the boundary conditions. The relationship between the smoothness of the boundary condition time variation and the convergence rate of the SCVS with respect to N m is demonstrated in Fig. 18. e M from Eq. (13), which measures the truncation error in the imposed boundary conditions associated with the finite N m , follows the same trend as the overall error e. As more modes are included in the solution, the boundary conditions approach their temporal reference profile, thereby improving the overall accuracy of the SCVS solution. Note that the two trends slightly diverge on the right tail of this plot. This divergence is an artifact of taking the case with N m = 11 as the reference when measuring e , whereas e M is measured against the reference solution. If the exact solution were available and employed as the reference solution, then e for very large N m would have converged to the larger of discretization e H and linear solver e L errors rather than going to zero.

Computational costs
The primary reason for the development of the SCVS algorithm was to obtain a more computationally efficient solution procedure. We discussed some of the performance metrics of the SCVS algorithm in the discussion pertaining to Fig. 13. In this section, we provide a more comprehensive account of the performance of these two solvers in terms of computational costs as well as wall-clock-time for all cases studied above.
For cases with multiple grids, we considered the finest grid in computing performance values. Also, we consider N m = 5 for the SCVS as it proved sufficient for sufficiently smooth boundary conditions. Thus, the computational cost t C for the SCVS is computed as the cumulative cost of simulating the first 6 modes. Its wall-clock-time t W , on the other hand, is taken as the maximum of the wall-clock-time of those 6 simulated modes since they are embarrassingly parallelizable. The simulation results presented in Sections 3.2 and 3.4 for the canonical 2D channel and 3D pipe flow cases were for individual modes. Nevertheless, we consider the sum of the cost of the first six modes for those cases as well. This choice was made to represent a situation in which the flow in these models is simulated with an arbitrary boundary condition rather than one with a unimodal oscillation. Furthermore, note that both t C and t W for the RBVMS scale linearly with the total number of time steps 1 . The results reported here are what is considered typical, namely simulating five cycles with each including 2,000 time steps for a total of 10,000 time steps. Finally, the performance metrics are in general implementation-and machine-dependent. However, given that both the RBVMS and the SCVS solvers are implemented by our group, written and compiled using the same language and compiler, linked against the same libraries, and ran on the same machine, the figures in Tables 3 provide an apple-to-apple comparison of the two formulations. Table 3: Comparison of the computational performance of the SCVS and RBVMS solvers. N p , t C , t W denote the number of processors, computational costs in CPU-hours, and the simulation wall-clock-time in hours, respectively. The last two columns are t C and t W for the SCVS relative to those of the RBVMS. All the SCVS results correspond to N m = 5. For cases with multiple grids, we adopted the finest grid in computing these performance figures. The relative figures are in percent.

Case
SCVS RBVMS Relative (%) As was hypothesized earlier, the SCVS reduces the overall computational cost significantly when it is compared against the standard RBVMS formulation. This reduction in cost that ranges from less than 1% for 2D cases to 11% for the 3D Glenn geometry is primarily achieved by reducing the number of linear solves. While this number is N m +1 = 6 for the SCVS, it is equal to the total number of time steps for the RBVMS, i.e., 10,000. The performance gap between the SCVS and RBVMS is not as large as 1667 = 10, 000/6 as a single SCVS linear solver is much more expensive than that of the RBVMS. On average, the SCVS requires 30,000 GMRES iterations, whereas 500 iterations suffice for solving the linear system obtained for the RBMVS formulation. The reason for the large number of iterations for the SCVS is mostly attributed to the high condition number of A matrix that has a zero diagonal sub-block (Eq. (10)). As discussed in Section 4, reducing the cost of solving this linear solver presents an opportunity for significantly reducing the cost of the SCVS algorithm in the future.
The performance gap between the SCVS and RBVMS widens when wall-clock time is concerned as t W for the SCVS is less than 4% of that of the RBVMS for all cases considered here ( Table 3). As we discussed earlier, the modes in the SCVS algorithm are independent and can be solved concurrently. Therefore, the number of processors utilized in the SCVS computation can far exceed that of the RBVMS without loss of parallel efficiency. Differently put, the wall-clock time for the SCVS is roughly independent of the number of computed modes as long as there are sufficient parallel computational resources available. The cumulative effect of this added dimension for parallelization and lower overall computational cost is a formulation that allows for a much faster turn-around time for a given problem.

Future work
The biggest hurdle to be overcome in the future is the extension of the SCVS formulation to high Reynolds number flows. Such an extension will not be trivial, as one can not exploit the linearity of the governing equations to solve for velocity and pressure at different modes independently.
Solving for all the modes at the same time will be one way to deal with this coupling. However, this brute force approach can become quickly unaffordable as the number of modes increases. This limitation, which has been observed in the past for the Time Spectral Method discussed in Section 1, must be overcome if the extension of SCVS to the Navier-Stokes equation were to find widespread use in the future.
The efficiency of the SCVS algorithm can be significantly improved in the future by making minor adjustments in its formulation or its underlying linear solver. Relaxing the incompressibility constraint by using a penalty method or including stabilization terms in its formulation is expected to reduce the condition of the underlying linear system significantly. Experimenting with more efficient linear solver algorithms (e.g., bi-partitioned method [30]) as well as preconditioners [29] can lead to additional improvements in the overall performance of the SCVS algorithm in the future.
Depending on the hardware architecture, reformulating the problem so that it avoids complex arithmetic may or may not reduce the cost of this scheme as such reformulation leaves the number of real floating-point operations unchanged.

Conclusions
In this paper, we proposed the SCVS as an alternative approach for fast simulation of time-periodic flow at low Reynolds numbers in complex geometries. Starting from the unsteady Stokes equation, we showed that it could be expressed as a steady Stokes equation with an imaginary source term in the time-spectral domain. The resulting equation can then be discretized and solved as a boundary value problem at a few selected modes, avoiding the use of costly and unscalable time integration schemes. As a proof of concept, we showed how this boundary value problem could be solved using Galerkin's formulation with mixed elements. We later employed this formulation for simulating flow in a variety of 2D and 3D geometries. To provide a point of reference, all these simulations were also performed using the standard RBVMS formulation.
For cases with an analytical solution available, the SCVS showed about an order of magnitude improvement in accuracy relative to the RBVMS at a similar number of degrees of freedom. This difference was attributed to the use of quadratic shape functions for the SCVS and to a lesser extent to the lack of stabilization terms and time integrator in our formulation. The variation of the overall error in the SCVS solution as a function of grid size, linear solver tolerance, and the mode number was measured through our numerical test cases and was shown to follow analyticallyderived power-law formulas. For cases with boundary conditions changing arbitrarily in time, we showed the overall error follows the same trend as the truncation error associated with simulating a finite number of modes. For smooth time-varying boundary conditions, such as those encountered in cardiovascular flows, we showed that using as few as 6 modes is sufficient for reducing the overall error to O(10 −3 ).
The SCVS led to significant improvement in performance in comparison to the RBVMS. The total computational cost of the SCVS for the cases considered in this study varied between 0.74% and 11% of that of the RBVMS, whereas its wall-clock time was consistently lower than 4% of that of the RBVMS. Further improvement in performance through the use of stabilized or penalty formulation as well as the extension to higher Reynolds numbers remain to be explored in the future.

Appendix A. RBVMS formulation
A brief account of the RBVMS formulation employed in this study is provided below. A more detailed description of this algorithm is provided in reference [6,7,30].
The weak formulation of the RBVMS is stated as follows. Find u ∈ S and p ∈ P, such that for all w ∈ W and q ∈ Q B G (w, q; u, p) + B S (w, q; u, p) = F (w, q) , In this equation, B G contains Galerkin's term whereas B S are the stabilization added to allow for equal-order velocity and pressure functions and prevent convective instability associated with Galerkin's method. Other parameters appearing in (A.2) are defined as in which c 1 = 1 and c 2 = 3 are the model constants, ξ ∈ R n sd × R n sd is covariant tensor obtained from a mapping between the physical and parent domains, and ∆t is the time step size. The resulting equations are discretized using triangular and tetrahedral elements in 2D and 3D. To relateu to u and integrate the resulting equations in time, the second order generalized-α method is adopted [36]. This algorithm, which is implicit in time, requires linearization of the nonlinear terms in Eq. (A.2). Thus, we employ the Newton-Raphson method for linearization of the resulting equations and a predictor-corrector procedure to converge on the solution at the next time step. The spectral radius of the infinite time step, which appears in the generalized-α time integration scheme, is set to 0.2 in this study. We verified that this parameter has no discernible influence on the results reported in this study. The working of the generalized-α time-integration scheme [36] is explained in greater detail below.
First, we predict unknowns at the next time step asu n+1 = γ−1 γu n , u n+1 = u n , and p n+1 = p n . Next, the unknowns are computed at an intermediate time, namelyu n+α m , u n+α f and p n+α f aṡ where, R M and R C are the momentum and continuity residuals and ∆u and ∆p are correction to our initial prediction of unknowns at the next time step.
The linear system in Eq. (A.5) is solved iteratively using the GMRES technique. Then corrections are applied usingu n+1 ←u n+1 + ∆u, u n+1 ← u n+1 + γ∆t∆u, Following these corrections, the Newton-Raphson iteration loop will be terminated if the norm of the residual is smaller than a specified tolerance. Otherwise, new intermediate variables are computed via Eq. (A.4) and a new iteration is performed. As mentioned before, for the linear Stokes equation, only one iteration is enough. This entire process is repeated at each time step.
Overall, the RBVMS algorithm requires three nested iteration loops: the outer loop for time-stepping, the first inner loop for Newton-Raphson iterations, and the innermost iterative loop for the linear solver. For the problems investigated in this study, however, the nonlinear terms in Eq. (A.2) can be neglected, resulting in two nested loops for the linear solver and time stepping.
Unlike operator splitting techniques that time advance velocity first and then enforce divergence-free condition by solving a Poisson equation, the RBVMS solves for both at the same time during a single linear solve as implied by Eq. (A.5). Despite this difference, the sub-block matrix that relates pressure to continuity (L in Eq. (A.5) that arises from u p · ∇q stabilization term in Eq. (A.2)) is identical to the discrete Laplacian operator that appears in the operator splitting technique. Note that this sub-block was zero in the linear system that arises from the SCVS formulation due to the lack of any stabilization term in our formulation (c.f., Eq. (10)). As a consequence of this zero sub-block, the tangent matrix is ill-conditioned, resulting in a lower rate of convergence for the SCVS in comparison to the RBVMS.

Appendix B. MSS formulation
To allow for a more direct one-to-one comparison against the SCVS technique, we have implemented unsteady Stokes equations in the time domain using mixed shape functions. The weak formulation of this solver (called MSS) is identical to the SCVS except for the traditional treatment of the acceleration term. Namely, we use a conventional time integration scheme (i.e., the second-order implicit generalized-α method explained in Appendix A) for this term. The weak form of the MSS states: find u ∈ S and p ∈ P, such that for all w ∈ W and q ∈ Q B G (w, q; u, p) = F (w, q) , (B.1) where B G =ˆΩ ρw ·u + ∇w : (−pI + µ∇ s u) + q∇ · u dΩ F =ˆΓ h w · hdΓ.

(B.2)
Note that the test function and solution spaces are identical to those of the SCVS except for being defined in real rather than complex domain. Equation (B.2) is discretized using triangular and tetrahedral elements in 2D and 3D that were used for the SCVS. We use the same time integration scheme (i.e., generalized-α) and linear solver as the RBVMS (i.e., GMRES). The MSS also requires two nested loops: an outer loop for time-stepping, and an inner loop for the linear solver.   (15)) and the cost of MSS with the SCVS and RBVMS solvers for W = 8π and t = T/2. The details of the meshes used for these computations are reported in Table 2, where QM1 is employed for the SCVS and MSS and LM1 was adopted for the RBVMS. The parameters of the time integration scheme are kept the same as the RBVMS, which were reported in Section 3.1.
The huge cost reduction from the MSS to SCVS can be explained by the fact that the SCVS solver performs just one linear solve, whereas the MSS solver performs 10,000 linear solves (each time step requiring one linear solve). However, the performance gap between the MSS and SCVS is 1522/1.25 = 1,220 that is less than 10,000 as the average number of GMRES iterations for the MSS is less than that of the SCVS (i.e., 2,336.7 versus 7,943). Having pure real operations as oppose to complex operations in the GMRES algorithm also reduces the performance gap between the two.
To explain the lower performance of the MSS in comparison to the RBVMS, we must analyze the structure of the linear system for these two solvers. Since there is no product of p and q in Eq. (B.2), the discrete form of the continuity equation will not depend on the pressure. This lack of dependence is in contrast to the RBVMS that through u p · ∇q term in Eq. (A.2) produces a non-zero L block in the tangent matrix in Eq. (A.6). This non-zero block reduces the condition number of the tangent matrix, leading to a much faster convergence rate of the GMRES algorithm for the RBVMS (493.97 for the RBVMS versus 2,336.7 for the MSS).
In terms of error, the SCVS that does not suffer from the time integration error and employs quadratic shape functions for velocity is the most accurate formulation. The SCVS is followed by the MSS solution that contains time integration error and the RBVMS that on top of that error uses linear elements.
The comparison of the SCVS against the MSS is kept limited to this appendix. For the remainder of this article, the SCVS is only compared against the RBVMS for three primary reasons. Firstly, as confirmed by the results shown in Table B.4, in general, the RBVMS performs better than the MSS, hence if the SCVS outperforms the RBVMS, it will also outperform the MSS. Secondly, the RBMVS is a more widely adopted formulation than the MSS, thus using it as a benchmark can provide a better point of reference. Thirdly, the cost and convergence properties of the MSS significantly suffer for larger meshes, making the simulation of some of the cases reported in Sections 3.4 and 3.5 nearly impractical. For these reasons, we benchmark the SCVS against the RBVMS for all cases discussed under Section 3.