Breaking the Kolmogorov Barrier in Model Reduction of Fluid Flows

Turbulence modeling has been always a challenge, given the degree of underlying spatial and temporal complexity. In this paper, we propose the use of a partitioned reduced order modeling (ROM) approach for efficient and effective approximation of turbulent flows. A piecewise linear subspace is tailored to capture the fine flow details in addition to the larger scales. We test the partitioned ROM for a decaying two-dimensional (2D) turbulent flow, known as 2D Kraichnan turbulence. The flow is initiated using an array of random vortices, corresponding to an arbitrary energy spectrum. We show that partitioning produces more accurate and stable results than standard ROM based on a global application of modal decomposition techniques. We also demonstrate the predictive capability of partitioned ROM through an energy spectrum analysis, where the recovered energy spectrum significantly converges to the full order model’s statistics with increased partitioning. Although the proposed approach incurs increased memory requirements to store the local basis functions for each partition, we emphasize that it permits the construction of more compact ROMs (i.e., of smaller dimension) with comparable accuracy, which in turn significantly reduces the online computational burden. Therefore, we consider that partitioning acts as a converter which reduces the cost of online deployment at the expense of offline and memory costs. Finally, we investigate the application of closure modeling to account for the effects of modal truncation on ROM dynamics. We illustrate that closure techniques can help to stabilize the results in the inertial range, but over-stabilization might take place in the dissipative range.

where the left and right set of basis functions match. To mitigate the instabilities encountered in Galerkin projection-based ROMs, Carlberg et al. [82] showed that least-squares Petrov-Galerkin can indeed build more accurate ROMs while Galerkin ROM becomes unstable for long time intervals. Alternatively, closure and stabilization models represent a feasible choice in the accuracy/stability-cost trade-off for Galerkin ROMs. Closure models have been popular in turbulence community, and found their way into ROM for fluid flows as well. Therefore, we also investigate the effect of a linear mode dependent eddy viscosity type closure model to stabilize the solution of ROM in case of 2D Kraichnan turbulence.

Two-Dimensional Kraichnan Turbulence
The two-dimensional (2D) turbulence cannot be realized in practice or experiments, but only in numerical simulations. Furthermore, 2D turbulence behaves differently from realistic 3D turbulence. In 3D turbulence, energy is transferred forward, from large scales to smaller scales, via vortex stretching. Conversely, in 2D turbulence, most of the forcing and dissipation energy will be transferred from smaller scales to larger scales and follows the Kraichnan-Batchelor-Leith (KBL) theory [83][84][85]. Despite the apparent simplicity in dealing with two spatial dimensions, 2D turbulence is very rich in its dynamics due to the inverse energy and forward enstrophy cascading mechanisms, which makes it a feasible testbed for a study of ROM in turbulence. Nonetheless, readers should be cautious while trying to extend the investigated framework into more realistic 3D turbulence cases.
The decaying homogeneous isotropic 2D Kraichnan turbulence follows the 2D Navier-Stokes equations, which can be written in vorticity-streamfunction formulation as below, where ω is the vorticity and ψ is the streamfunction. Re is Reynolds number relating the inertial and viscous effects. J(ω, ψ) and ∇ 2 ω are the Jacobian and Laplacian operators, respectively, which can be defined as The vorticity-streamfunction formulation enforces the incompressibility condition, where the vorticity and streamfunction fields are related by the following Poisson equation, According to the KBL theory for 2D turbulence, the inertial range in the energy spectrum is proportional to k −3 in the inviscid limit. In our numerical experiments, the initial energy spectrum in Fourier space is given by where k = k 2 x + k 2 y and k p is the wavenumber at which the maximum value of initial energy spectrum occurs. During the time evolution process, due to the nonlinear interactions, this spectrum quickly approaches toward k −3 spectrum as we will show later when we present our results. The magnitude of the vorticity Fourier coefficients is related to the energy spectrum as Then, the initial vorticity distribution (in Fourier space) is obtained by introducing a random phase. More details regarding derivation of the initial vorticity distribution from an assumed energy spectrum can be found in [86]. In the present study, we use a spatial computational domain of (x, y) ∈ [0, 2π] × [0, 2π] and a time domain of t ∈ [0, 4]. Periodic boundary conditions are applied in both xand y-directions. A spatial grid of 1024 2 and a timestep of ∆t = 0.001 are used. A fourth-order accurate Arakawa scheme [87] is used for spatial discretization and a third-order total variation diminishing Runge-Kutta scheme (TVD-RK3) [88] is adopted for temporal discretization. The energy spectrum for different values of Reynolds number is given in Figure 1 at t = 2 and t = 4. We can observe that the energy spectrum in the inertial range flattens toward the classical k −3 scaling limit as Re increases, in agreement with the KBL theory for two-dimensional turbulence. It should be noted that the same initial vorticity distribution is enforced for all numerical experiments, by fixing the seed of the random numbers generator. In order to verify the validity of the used spatial resolution, Figure 2 shows the energy spectrum for Re = 16,000 (left) and Re = 32,000 (right) using two spatial resolutions (i.e., 1024 2 and 2048 2 ). A good match between the energy spectra corresponding to the two spatial resolutions can be seen, even for Re = 32,000. Nonetheless, we will consider Re = 16,000 with a spatial grid of 1024 2 for the rest of the present study. In addition, results obtained from solving Equation (2) using the aforementioned numerical schemes and resolution will be denoted as the full order model (FOM) results.
(a) (b) Figure 1. Evolution of energy spectra in 2D decaying Kraichnan turbulence for Re ∈ {2000, 4000, 8000, 16,000, 32,000 } at t = 2 (a) and t = 4 (b). Notice the convergence of energy spectrum in the inertial range to the k −3 scaling with increasing Reynolds number.
(a) (b) Figure 2. A mesh independence test for Re = 16,000 (a) and Re = 32,000 (b) using two spatial resolutions. An acceptable match is obtained from both resolutions.

Reduced Order Modeling
In order to build a Galerkin projection-based reduced order model (ROM) of a physical system, the system's state, u(x, t), is assumed to live in a reduced trial subspace. Then, the high-dimensional operators are projected onto the same subspace, constraining the residual to be orthogonal to that subspace. Therefore, generating effective reduced subspaces that sufficiently contain the solution manifold is vital for the successful construction of Galerkin ROMs. One of the very early-developed and well-known approaches to build reduced subspaces is Fourier analysis. However, it assumes universal basis functions (or modes) which have no specific relation to the system in hand, except for boundary conditions. More efficiently, snapshot-based model reduction techniques can tailor a reduced subspace that best fits the problem by extracting the underlying coherent structures that control the major dynamical evolution we are interested in. Among those snapshot-based techniques, proper orthogonal decomposition (POD) is the most popular and well-established approach extracting the modes which contribute the most to the total system's energy [4,89]. The basic implementation of POD is presented in the following subsection.

Proper Orthogonal Decomposition
As mentioned before, proper orthogonal decomposition (POD) is a snapshot-based approach. Therefore, a collection of system's realizations (called snapshots) are required in advance. Assuming that a number of snapshots (N s ) of the 2D vorticity field, denoted as ω(x, y, t n ), are stored at consecutive times t n for n = 1, 2, . . . , N s , the time-averaged field, called "base flow", can be computed as The mean-subtracted snapshots, also called anomaly or fluctuation fields, are then computed as the difference between the instantaneous fields and the mean field, ω (x, y, t n ) = ω(x, y, t n ) −ω(x, y).
This mean-subtraction is common in ROM community, and it guarantees that ROM solution would satisfy the same boundary conditions as a full order model [90]. This anomaly field procedure is also interpreted as a mapping of snapshot data to its origin (i.e., centering), and will be used for constructing the POD basis functions. The POD approach basically relies on an eigenvalue decomposition of a snapshot correlation (or covariance) matrix. Since the number of snapshots (O(10 2 )) is usually much less than the number of grid points (O(10 6 )), we follow the method of snapshots [54] for efficient implementation of the POD algorithm. An N s × N s temporal correlation data matrix A = [a ij ] is computed from the inner product of temporal anomaly fields as where the angle-parenthesis denotes the inner product in Hilbert space defined as where Ω represent the physical (spatial) domain. Then, an eigenvalue decomposition of A is performed as where Λ is a diagonal matrix whose entries are the eigenvalues λ k of A, and V is a matrix whose columns v k are the corresponding eigenvectors. A is a Gramian matrix, a property which can be utilized for efficient solution of the eigenvalue problem (e.g., using Jacobi transformations [91]).
The eigenvalues λ k represent the amount of information contained in each mode; therefore, they should be arranged in a descending order (i.e., λ 1 ≥ λ 2 ≥ · · · ≥ λ N ), for proper selection of the POD modes. The POD modes φ k are finally computed as where v n k is the n-th component of the eigenvector v k . It can be easily verified that the POD modes are orthonormal (i.e., φ i ; φ j = δ ij ), which simplifies the derivation of ROM.

Galerkin Projection
With the POD algorithm implemented, a set of POD basis functions are obtained {φ k (x, y)} R k=1 , where R is the number of POD modes retained in ROM (i.e., the size of the reduced subspace). The solution is therefore approximated as a superposition of the contributions of these selected modes as follows: where α k (t) are the so-called temporal modal coefficients. Since the vorticity and streamfunction are related by the kinematic relationship given by Equation (5), the basis functions (θ k (x, y)) and mean field (ψ(x, y)) corresponding to the streamfunction can be obtained from those of the vorticity as follows: which might result in a set of basis functions for the streamfunction that are not necessarily orthonormal. Moreover, the reduced order approximation of streamfunction can share the same temporal coefficients α k (t), Those approximations are then substituted in Equation (2), as follows: Then, we use the definition of Jacobian operator and expand the above expression as Note that the first term ∂ ∂tω (x, y) is zero. In addition, the Jacobian and Laplacian operators are spatial functions; therefore, the following can be written (after change of indices) Then, we can collect the constant, linear, and nonlinear terms (in time) as follows: Finally, an inner product with an arbitrary basis function φ k (x, y) is performed, where we make use of the orthonormality of the vorticity basis functions. Therefore, a reduced order model is obtained, which can be simply represented by the following set of ODEs for α k (t), where the predetermined model coefficients can be computed by the following numerical integration (offline computing) To complete the dynamical system given by Equation (18), the initial conditions for α k (t) may be obtained by the following projection of the initial vorticity field onto the respective basis functions, where ω(x, y, t 0 ) is the vorticity field specified at initial time t 0 . We note that the approach described here is called tensorial ROM, since the model predetermined terms (especially those corresponding to nonlinear terms) are computed in advance as tensors during an offline stage, while the online solution of Equation (18) scales with R 3 . Alternatively, the model terms can be calculated online while minimizing the cost of computing the nonlinear terms through approximating approaches like discrete empirical interpolation method (DEIM) [92,93], gappy POD [60,94], and missing point estimation (MPE) [95,96]. An overview of these techniques can be found in [97,98].

Partitioned ROM
A representative solution subspace is a key factor in the construction of ROMs of physical systems and POD has been successfully offering a systematic approach for building such a subspace. However, POD in principle depends on a global approximation of the snapshot data. This has several consequences, a few of which are mentioned here. First, local temporary excursions in state space which contain small amounts of energy can be overlooked by POD because their contribution to the total energy may be negligible. These excursions can, however, be of interest and have significant impact on the dynamical evolution (see [99], for example). Second, the global nature of POD approach can result in overall deformation of the obtained modes for systems with fast variations in state. As a result, the constructed POD modes are smoothed out and cannot capture any dominant structure at all. This in turn causes an energy distribution across a large number of modes; therefore, the truncated modes contain a significant amount of system's energy, which causes inaccuracies and instabilities in the solution unless the size of the ROM is increased.
An alternative approach which preserves the optimality of POD, while giving care to the system's local characteristics, is called principal interval decomposition (PID) [69][70][71][72][73]. In PID, the time domain is divided into a few partitions, each characterizing a specific stage in the system's dynamics and evolution. Then, the same POD algorithm is applied within each partition to generate a set of basis functions that best fit the respective partition locally. In a sense, the PID approach can be viewed as decomposing the global POD subspace into a few of locally optimal subspaces. As a result, accurate partitioned ROMs with smaller sizes can be dedicated to each individual sub-interval. More details about the properties of PID can be found in [73]. This partitioning idea can also be performed in the physical domain, state space, or parameter space [64,[100][101][102].
In brief, the time domain T = [t 0 , t f ] is divided into a number P of non-overlapping partitions (every two consecutive partitions share only the interface). Adaptive partitioning and clustering techniques can be used to tailor partitions with arbitrary widths if no a priori information is available about the dynamics of the systems. However, for simplicity and clarity of presentation, we assume equidistant decomposition of the time domain. In other words, the total number of snapshots N s is divided into N s /P sets. After that, local mean fieldsω (l) (x, y) andψ (l) (x, y) are computed and the POD algorithm described in Section 3.1 is applied to each group of snapshots to generate local basis k (x, y)} R k=1 for l = 1, 2, . . . , P. Then, the Galerkin ROM equations are obtained for each interval as where the predetermined model coefficients can be computed partition-wise by the following numerical integrations: and the initial conditions for α (1) k (t) are obtained by the following projection: It remains to enforce the interface constraints to enable solution transition from one partition to the next, if no access to the true field at the interface instant is available. We do so by imposing the following condition at the interface [103], where t l represents the time instant at the interface between the interval l and l + 1. Here, ω (l) (x, y, t l ) is the reconstructed vorticity field at the interface using the left interval (l), while ω (l+1) (x, y, t l ) is reconstructed at the same instant, but using the right interval (l + 1). Those can be rewritten as This reduces to performing the following computation at the interface to re-initiate our solver at the first timestep of the new time interval,

Results
As indicated in Section 2, the full order model is solved over a square domain with a side length of 2π using a spatial grid of 1024 2 for t ∈ [0, 4] with a timestep of 0.001. These were shown to be sufficient resolutions to solve the problem at a Reynolds number of 16, 000. To build the reduced order subspace using the POD and PID algorithms, described in Sections 3.1 and 4, 800 snapshots of vorticity field were collected (i.e., every five timesteps). It was also shown that the sampling rate has minor effect on POD-based modal decomposition, as long as the Nyquist-Shannon sampling criterion is met [104]. Another factor we consider in our choice of sampling rate is to have adequate number of snapshots in each sub-interval (after partitioning the time domain).
The partitioned ROM approach is used to solve Equation (21) using different numbers of partitions. In particular, we present the results for P ∈ {1, 4, 8}. Note that, for P = 1, the approach reduces to the standard Galerkin-POD ROM with a single interval covering the whole time domain. In Figure 3, the energy spectrum at t = 4 is used to illustrate the effect of increasing the number of partitions on ROM accuracy and convergence to the FOM solution. This is shown for P = 1 (top), P = 4 (middle), and P = 8 (bottom), along with a compensated (for zooming-in) view of the energy spectrum at the inertial range on the right. We can easily observe that increasing the number of partitions improves the convergence of ROM energy spectrum to the true one for all values of R.
Although increasing the number of retained modes generally improves the ROM accuracy, this does not help a lot in the case with P = 1. This is shown in Figure 3, where we still see a large deviation of ROM results in the inertial range even with 32 modes. This might be attributed to the smoothing-out of constructed basis functions in such a way that information is distributed across a large number of modes. Consequently, the effect of modal truncation on system's dynamics is no more negligible. On the other hand, a partitioned ROM with P = 8 enables us to build accurate ROMs with relatively compact sizes (e.g., R = 16). To provide more insights about the scales resolved with partitioned ROM, we plot the evolution of the temporal coefficients of the first and eighth modes in Figure 4. We can notice that a ROM with global POD modes (i.e., P = 1) can adequately capture the evolution of the first mode's contribution (with increasing deviations toward the end of the time domain). On the other hand, a phase shift and amplification of the temporal evolution of the 8 th modal coefficient are spotted after a small initial period. From a physical point of view, the first modes correspond to the large dominating scales, while the last ones correspond to the small (or fine) scales. This indicates that a global POD application provides ROMs that moderately capture the large scales, but can fail to capture the finer details. Moreover, we can deduce from Figure 4a that the dynamics at different time instants is different (where all the snapshots are projected onto the same set of basis functions). This justifies the use of paritioned ROM to tailor unique bases for each sub-interval. On the other hand, with partitioned ROMs (i.e., P = 4 and P = 8), the evolution of both α 1 and α 8 is tracked accurately, implying that the presented partitioning approach improves the ROM's capabilities to represent the fine scales. Moreover, the accuracy of α 1 prediction is improved compared to P = 1. Although the evolution of temporal coefficients in different sub-intervals looks similar, this is because the fields in each region are projected onto their respective bases. Therefore, once reconstructed, the dynamics in the full order space will be different. In addition, we can note the "jump" at the interfaces (e.g., at t ∈ {1, 2, 3, 4} for P = 4) corresponding to updating the solution manifold. The true modal coefficients (denoted as FOM) is obtained from projecting the vorticity field at different times on the corresponding basis function as α (l) More visualizations of ROM's predictability are provided in Figures 5 and 6, showing the vorticity (in color) and streamfunction (in height) at the final time (i.e., t = 4), compared to the FOM results. In Figure 5, we illustrate that ROM with P = 1 fails to correctly predict the final fields at different values of R. At R = 16, even the large scales are not captured accurately (recall the larger deviation of α 1 at final time from Figure 4). Conversely, both the large fine scales are predicted reliably using the partitioned ROM approach as shown in Figure 6 with R = 16.
In order to obtain more insights about the scales resolved by partitioned ROM, compared to standard POD-based ROM, we compare PID ROM results for P = 8 and R = 16, with standard ROM (i.e., P = 1) and increased number of modes. The root-mean-squares errors between Galerkin ROMs and FOM are shown in Table 1. We find that at least around R = 128 modes are needed for standard Galerkin ROM to achieve similar accuracy as partitioned ROM, which suggests that the relevant scales represented in those global 128 modes are almost encapsulated in 16 local modes if eight partitions are used. The computational gain of the PID approach can be easily seen from this analysis since online computational cost scales with O(R 3 ).

Kolmogorov n-Width and Relative Information Content
The Kologorov n-width provides a measure of the system's reducibility, and in POD/PID context, it can be considered as a measure of how well a linear superposition of POD modes might represent the underlying dynamics. Similarly, a relative information content (RIC) of the snapshot matrix can be used as a simple quantitative metric to understand the Kolmogorov n-width of the system. The following RIC formula [105] is used to compute the percentage modal energy, In other words, RIC(k) represents the fraction of information (variance) of the snapshot matrix that can be recovered using a specific number (k) of basis functions. Figure 7 shows the growth of RIC values with increasing k for different values of P. We can observe that, in case of P = 1, the increase of the fraction of preserved information is relatively slow and more modes are required to exceed a 90% RIC. On the other hand, with the use of a partitioned ROM approach, the slope of RIC curve increases significantly. For example, less than 10 modes are enough to hit the 90% RIC value. In Kolmogorov n-width context, a faster increase in RIC corresponds to a faster decrease of the Kolmogorov n-width. Therefore, we consider that partitioned ROM can, in turn, help to break the Kolmogorov barrier and increase the system's reducibility. We note that, in Figure 7a, we only show the RIC in the first sub-interval (i.e., l = 1) for the sake of clarity. In Figure 7b, we plot the RIC curves for P = 4 in the four sub-intervals, where they almost match each other. An interesting observation from Figure 7b is that the slopes of RIC in the second and third sub-intervals are the lowest, compared to those in the first and last sub-intervals. This might be justified by the fact that flow is just initiated in the first sub-interval with relatively smooth dynamics. After that, most of the system's vortical evolution and energy exchange take place within the following sub-intervals, decaying to quasi steady-state in the last sub-interval. For fast evolution and strong dynamics, more modes are required to capture sufficient amount of information. Alternatively, constructing local modes (through further partitioning) can help to accurately capture those dynamics. Since we use fixed number of modes in each sub-interval, this might correspond to different RIC values. That is the ROM approximation in some regions might be more or less accurate than the approximation in other regions. This might be crucial for flow problems with stronger and more complicated dynamics, where that difference in RIC values can cause inconsistencies of approximations. Alternatively, varying values of R can be used to satisfy pre-determined level of RIC in each partition.

Closure Modeling
It turns out that the partitioning approach helps with building more accurate and stable ROM due to several effects. First, it helps to construct more compact and localized subspaces which capture fine scales efficiently. Second, the energy (information) is concentrated in fewer modes. As a result of these two effects, modal truncation works effectively as the discarded modes have negligible effects on the system's dynamics. However, as shown in previous discussions, a large number of partitions might be needed to provide acceptable accuracy. Other than the increased storage requirement to store basis function for each sub-interval, a projection error is accumulated due to interface treatment. If the true field data are available at the interface, it can be projected onto the respective basis functions to obtain the corresponding modal coefficients at the beginning of the sub-interval. Since this is not generally available, the update of ROM solution manifold is performed using Equation (27), where the obtained reduced order approximation of the field from the previous partition is projected onto the basis functions of the following partition. Since that approximation depends on the accuracy of the ROM in the current partition, and lives in a reduced subspace different than that of the following subspace, discontinuities can easily occur at the interface. In Figure 8, we show the modal coefficients at the interfaces for P = 4 obtained from ROMs with different R dimensions. Moreover, the mismatch between initial conditions at each interface increases in subsequent partitions, due to the successive interface treatment, resulting in amplification of the error. However, we highlight that the projection error generated at the interface is bounded and overall accuracy gain due to the localization in the PID approach surpasses that error. Conversely, for ROM with global POD application without any partitioning, the truncated modes might include substantial amount of energy and/or have strong interactions with the retained modes. Hence, inaccuracies and instabilities can easily take place in the produced ROM. Therefore, a compromise between accuracy gain from localization and projection error from interfaces should be considered while selecting the number of intervals. To enhance the stability of ROMs, closure models have been effectively incorporated (inspired from large eddy simulation (LES) studies). Here, we investigate the effect of closure modeling to stabilize and correct ROMs based on global application of POD (i.e., P = 1). We utilize a linear mode length eddy viscosity closure based on Rempfer's model [72,106,107], where the 2D Navier-Stokes equation (given in Equation (2)) is re-written as where κ is a free-parameter controlling the amount of stabilization added to the model, R is the number of modes, and k is the mode index (i.e., k = 1, 2, . . . , R). More recently, several ideas from a plethora of predictive turbulence models have been used in ROMs in the spirit of closure modeling to account for the effects of the truncated scales [108][109][110][111][112][113][114][115]. In the present study, we test the effect of the linear mode dependent eddy viscosity closure model for κ ∈ {0, 4, 8}, where κ = 0 corresponds to the standard Galerkin ROM without stabilization.
In Figure 9, we present the temporal evolution of α 1 and α 8 for ROM with 16 modes and a single global interval (i.e., P = 1). We find that the closure significantly improves the predictive capabilities of the ROM, and with increasing κ (i.e., adding more stabilization), results converge to the true values of the modal coefficients. In addition, the final vorticity and streamfunction fields at t = 4 are shown in Figure 10. We can notice tremendous improvement of predicted results with closure modeling. With close investigation of surface topology, we can observe that Figure 10d captures the correct large scales as Figure 10a. However, we can see that smaller and finer details are overlooked and smoothed-out. This is to be expected given the fact that closure models try to compensate the effect of the truncated modes on the dynamics of the retained modes. In other words, it helps to improve our predictions of the dynamics of these retained modes, but does not add the contribution of the truncated modes into our reduced order approximation. Therefore, our ROM prediction will always be restricted to the scales represented by the retained modes. In case of global POD modes, only the large scales are captured and the fine scales are smoothed-out, which justifies Figure 10d. Finally, the energy spectrum at final time is given in Figure 11, where we can see the influence of closure modeling in improving ROM results. A very nice observation from this figure is that stabilization helps a lot in the inertial regime of energy spectrum, but it causes over-stabilization in the dissipative regime. The root-mean-squares errors (RMSE) between Galerkin ROM approximation (with P = 1 and R = 16), in the presence of stabilization, and FOM solution are shown in Table 2, compared to the base case without closure. It is clear that the addition of stabilization to Galerkin ROM helps to improve the predictions and reduces the approximation error. However, compared to Table 1, we can deduce that the accuracy gain due to partitioning is much more superior than that of closure. In addition, the RMSE is sensitive to the value of the stabilization parameter. Therefore, a dynamic selection of this parameter might be needed to guarantee sufficient correction.
Moreover, it should be noted that there might exist situations where the gain from only closure is minimum, even in case of "optimal" closure [73]. Hence, a combination of both partitioning and stabilization seems to be a feasible choice for building accurate, stable, and compact ROMs for turbulent flows. Indeed, Ahmed and San [72] demonstrated that an eddy viscosity based closure in conjunction with PID yields significant improvements in accuracy for capturing strong moving discontinuities with a negligibly small computational overhead. In addition, Table 3 shows the effect of stabilization through closure modeling in the partitioned ROM (with P = and R = 16). It is clear that the closure slightly improves the RMSE of vorticity field approximation. In addition, we can see that further increase in κ can actually increase the RMSE due to over stabilization. It can be also observed that RMSE in streamfunction predictions in case of closure modeling with PID is generally larger than that in case of PID without closure. This can be attributed to the nature of streamfunction being inverse of the Laplacian operator of vorticity, making it a smoother field by construction. Therefore, closure with PID might cause over-stabilization and excessive smoothing of the streamfunction fields predictions. This again suggests the development of a dynamical strategy to estimate optimal values of eddy viscosity coefficient. Table 2. Root-mean-squares error at time t = 4 between the standard Galerkin ROM (P = 1, R = 16) and FOM in order to demonstrate the sensitivity with respect to the additional stabilization through κ at Re = 16,000.

Computational Cost
In addition to the accuracy gains through the partitioned ROM approach, presented in previous subsections, a comparison in terms of computational cost is required for fair assessment. Here, we consider two phases of ROM computations: offline and online. By offline, we refer to the evaluation of the vectors, matrices, and tensors corresponding to the constant, linear, and nonlinear predetermined model coefficients (i.e., B, L, and R). Meanwhile, we denote the actual prediction and solution of Equation (21) as the online stage. Table 4 provides the CPU times in the offline and online stages for different number of partitions and modes. From this table, we can see that increasing the number of modes increases the CPU times for both the offline (due to computation of larger vectors, matrices, and tensors) and online (due to solving an increased number of equations) stages. Interestingly, we can observe that the computational cost of the online stage scales cubically with increasing R. This conforms to the quadratic nonlinearity in the investigated system, making the online cost an order of O(R 3 ).
On the other hand, increasing the number of partitions has a negligible effect on the online CPU time. Therefore, the partitioned ROM approach provides more accurate results with similar online costs as standard ROM with the same number of modes. However, increasing P (with fixed R) results in an increase in the offline CPU time. This corresponds to the computation of B, L, and R for each individual partition. That is why the offline CPU time scales linearly with increasing P. Nonetheless, it was shown in [73] that the CPU times for basis construction are largely reduced in the PID framework due to the solution of smaller eigenvalue problems rather than solving a large one. Finally, it is noteworthy that the partitioned ROM approach requires more memory to store the local basis functions and model predetermined coefficients in each partition. Therefore, we can consider that partitioned ROM converts the offline and storage costs into a speed-up of the online deployment.

Conclusions and Future Works
We utilize a partitioning approach in the current study in order to build a compact reduced order models of 2D turbulent flows. We demonstrate the framework using the two-dimensional decaying Kraichnan turbulence. It is shown that this partitioned ROM helps to break the Kolmogorov n-width barrier, in the context of a relative information content criterion. An energy spectrum analysis is performed to emphasize the capabilities of the presented framework to sufficiently capture both the large and fine scales and recover the energy spectra within the inertial and dissipative ranges. The partitioned ROM produces significantly more accurate predictions than standard ROM without partitioning. A computational CPU time analysis reveals that the partitioned ROM minimizes the online cost at the expense of offline cost and increased memory to build reduced order approximation with acceptable accuracy. In addition, we investigate the use of a linear mode dependent eddy viscosity to stabilize ROM prediction. We find that closure improves the results in the inertial kinetic energy regime, corresponding to the large scales. However, it presents over-stabilization in the dissipative regime, which smoothens out the field finer details. This might advocate a strategy to develop optimal (or near optimal) estimates for the eddy viscosity coefficients to be posed. Therefore, we suggest the use of partitioning along with closure as an enabler for more efficient ROMs turbulent flow systems. Incorporating PID with closure would provide the flexibility of constructing ROMs with less modes (and same number of partitions), or with less partitions (and same number of modes). The former would provide higher speed-ups during online deployment, while the latter would reduce the offline cost.
To this end, the partitioned ROM has been shown to improve the predictions for the investigated 2D Kraichnan turbulence case. However, scaling this framework into a 3D setting is still required to fully assess the validity of the framework to address realistic turbulent flows with strong anisotropy. In addition, the PID has been shown to capture finer details than standard global POD. Application of this approach onto flow situations where global basis fails drastically and results in completely unstable and inaccurate ROMs can be an interesting extension to the current study. Also, the adopted PID approach assumes evenly-spaced partitions with an equal number of modes. Despite the resulting simplicity of implementation, this fixed partitioning might be less meaningful for more challenging flow problems where dynamics change rapidly. Therefore, an adaptive partitioning with varying number of modes can be another extension to generalize the present work. The assumption of non-overlapping partitions might also be relaxed to allow overlapping sub-intervals with an indicator for smooth transition between partitions based on the actual flow conditions. Finally, the development of a strategy to dynamically update the basis functions using data assimilation techniques would enable the extrapolation of ROM in time while reflecting local phenomena.