3.1. Modal decomposition

First, we consider models generated using POD modes. POD produces economic reduced-order models, but has the well-known shortcoming of mixing together fluid motions at different temporal/spatial scales (Mendez, Balabane & Buchlin Reference Mendez, Balabane and Buchlin2019). Second, we consider models generated from modes oscillating at a single frequency obtained from a procedure that is equivalent to a DFT of the velocity snapshots. For practical convenience, we obtain the two distinct sets of modes using the same computational technique, based on the approach proposed by Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016) which only operates on the temporal correlation matrix. Briefly, the method considers the temporal correlation matrix ${\boldsymbol{\mathsf{R}}} \in \mathbb {R}^{N_t, N_t}$, with entries

(3.1)\begin{equation} {\mathsf{R}}_{ij} = \frac{1}{N_t}( \boldsymbol{u}^{\prime}(t_i, \boldsymbol{x}), \boldsymbol{u}^{\prime}(t_j, \boldsymbol{x}) ), \end{equation}

and then defines a filtered correlation matrix ${\boldsymbol{\mathsf{S}}}$, with elements

(3.2)\begin{equation} {\mathsf{S}}_{ij} = \sum_{k = -N_f}^{k = N_f} g_k {\mathsf{R}}_{i+k,j+k}, \end{equation}

given by the application of the filter coefficient vector $\boldsymbol {g}$ along the diagonals of the correlation matrix. An ordered set of temporal coefficients $\boldsymbol {a}_i = [a_i(t_1), \ldots , a_i(t_{N_t})]$ and associated mode energies $\lambda _i$ is then obtained from the eigendecomposition of ${\boldsymbol{\mathsf{S}}}$,

(3.3)\begin{equation} {\boldsymbol{\mathsf{S}}} \boldsymbol{a}_i = \lambda_i \boldsymbol{a}_i, \end{equation}

so that $\lambda _i\delta _{ij} = \boldsymbol {a}_i^{\top }\boldsymbol {\cdot }\boldsymbol {a}_j$. As discussed in Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016), when the filter is extended over the entire dataset and in the limit of number of samples tending to infinity, the filtered correlation matrix converges to a Toeplitz, circulant matrix. Then, its eigenvalues trace the power spectral density of the underlying dataset. On the other hand, the eigenvectors $\boldsymbol {a}_i$ corresponds to the Fourier basis. This procedure generates conjugate pairs of modal structures with same energy oscillating at a single frequency. These can be viewed as a set of modal oscillators exhibiting periodic fluctuations (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017) and tracing fluid motion at on a two-dimensional subspace. In practice, for a finite-length dataset, we filter the temporal correlation matrix assuming periodicity using a box-car filter, as suggested in Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016). Hereafter, we will refer to the modal structures identified by this procedure as DFT modes.

One important consideration is that, unlike dynamic mode decomposition (see Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Schmid Reference Schmid2010; Chen, Tu & Rowley Reference Chen, Tu and Rowley2012), DFT lacks the ability to discern and identify dominant frequency components. Instead, a number of modes equal to the number of snapshots utilised is produced, oscillating in conjugate pairs at specific frequencies determined by the sampling period ${\rm \Delta} t$ and observation time $T$ (Mendez et al. Reference Mendez, Balabane and Buchlin2019). This property, picket fencing, results in frequencies that are integer multiples of the fundamental frequency $f_1 = T^{-1}$, up to the Nyquist component $f_{Nyq} = (2{\rm \Delta} t)^{-1}$. In addition, unlike for POD, as the length of the dataset is increased, the number of energy-relevant modes increases and low-frequency modes with little dynamical importance appear. The approach we use here is to divide the dataset into five partition of thirty time units, covering an average of 20 cycles of the dominant oscillatory component, and providing sufficient frequency resolution to distinguish small scale spectral features. In addition, two possible ways of sorting pairs of modal structures are possible, i.e. by energy content (using the eigenvalues $\lambda _i$) or by frequency. Models obtained with the two sorting schemes will be referred to as $\mathrm {DFT}_e$ and $\mathrm {DFT}_f$, respectively.

We now focus on the characteristics of the modal structures obtained by these two methods. We denote the normalised cumulative sum of the eigenvalues $\lambda _i$ of the (filtered) correlation matrix as

(3.4)\begin{equation} e(n) =\left. { \sum_{i=1}^{n} \lambda_i }\right/{ \sum_{i=1}^{N_t} \lambda_i}, \end{equation}

describing the fraction of the fluctuation kinetic energy captured by the first $n$ elements of the expansion (2.2).

This quantity is shown in figure 2(a) for the POD and for the two possible DFT sorting schemes. As expected, a larger energy is captured by the POD basis. For the DFT decomposition, the energy-based sorting is more efficient at data compression, although the difference vanishes for large $n$, since for low-energy modes the two sorting schemes are equivalent. The modal energies associated with the $\mathrm {DFT}_f$ modes are shown in figure 2(b) as a function of the modal index $i$. The distribution is characterised by a continuous component, with modal energy decaying with frequency, and a discrete component, with a fundamental peak for the pair of modes (31, 32) and its first few harmonics. The peak, at a non-dimensional frequency $f = 0.7$, is physically originated from the high-energy structures transported along the shear layer by the rotation of the main vortical structure.

Figure 2.

(a) Cumulative sum of the first $100$ eigenvalues of ${\mathsf{S}}_{ij}$ for the three decompositions considered. (b) Distribution of the modal energies of $\mathrm {DFT}$ modes sorted by frequency.

This can be observed in figures 3(a) and 3(b), showing the vorticity field $\omega$ of the DFT mode pair (31, 32).This pair of modes describes a vorticity perturbation having the form of a wave travelling along the edge of the main vortex. Hence, the spatial structure of the two modes is shifted in the direction of the shear layer by half-wave. Travelling-wave structures in cavity flows have already been observed in simulation by Poliashenko & Aidun (Reference Poliashenko and Aidun1995); Auteri et al. (Reference Auteri, Parolini and Quartapelle2002) and characterised by global stability analysis and Koopman analysis by Boppana & Gajjar (Reference Boppana and Gajjar2010) and Arbabi & Mezić (Reference Arbabi and Mezić2017), respectively. The two leading POD modes, reported in figure 3(c,d), have the same energy and capture the same travelling-wave pattern described by the leading DFT mode pair.

Figure 3.

Vorticity field of the most energetic pair of DFT modes, panels (a,b), and of the first two POD modes, panels (c,d).

3.2. Energy analysis

To provide a more robust foundation to understand the sparsification results reported in §§ 3.3 and 3.6, we first focus on the analysis of the average energy interactions.The structure of the interaction tensor ${\boldsymbol{\mathsf{N}}}$ for a large POD-based model with $N=75$, reconstructing more than $99\,\%$ of the fluctuation kinetic energy, is reported in figure 4, showing the magnitude of the interactions for three slices for $i=1, 10$ and $75$, in panels (a), (b) and (c), respectively. All entries of the tensor ${\boldsymbol{\mathsf{N}}}$ are generally non-zero, although the strength of the interactions varies across several orders of magnitude. This is a combined result of the projection coefficients tensor ${\boldsymbol{\mathsf{Q}}}$ (shown later), whose entries are typically non-zero, and of the complex spectral structure of the temporal coefficients $a_{i}(t)$. The most important feature of figure 4 is that interactions are highly organised and there exists a subset of interactions that are more active. Specifically, for any mode $i$, triadic interactions can be classified as illustrated in panel (a) in four different categories by introducing a cutoff modal index $n$. The subset of interactions denoted as $\mathrm {LL}$ corresponds to nonlinear energy transfer involving pairs of low-index modes, $\mathrm {HL}$ and $\mathrm {LH}$ denote interactions involving high–low/low–high index modes, while $\mathrm {HH}$ denotes the subset of interactions involving pairs of high-index modes. We observe that the areas corresponding to $\mathrm {LL}$ and $\mathrm {HL}/\mathrm {LH}$ are the most active. If we map low/high modal indices to large/small scales, this result is in agreement with the picture of energy transfer between scales in homogeneous isotropic two-dimensional turbulence (Ohkitani Reference Ohkitani1990; Laval et al. Reference Laval, Dubrulle and Nazarenko1999), where the large scales interact with the small ones in a non-local fashion. In addition, interactions are not symmetric with respect to a swap of indices $j, k$. This can be quantified by computing the coefficient

(3.5)\begin{equation} \left. \chi_i(n) = \sum_{j=1}^{n} \sum_{k=1}^{N} {\mathsf{N}}_{ijk} \right/ \sum_{j=1}^{N} \sum_{k=1}^{n} {\mathsf{N}}_{ijk}, \end{equation}

representing the relative dynamical importance of the subset of interactions $\mathrm {LL}+\mathrm {HL}$ and $\mathrm {LL}+\mathrm {LH}$. Figure 4(d) shows $\chi _i$ for $i = 1,10$ and $75$ as a function of the normalised cutoff $n$. The interaction subset $\mathrm {HL}$ is up to four times more important than the subset $\mathrm {LH}$. This is a consequence of the asymmetry of the projection coefficients ${\mathsf{Q}}_{ijk}$, which arise from the fact that the convective transport of structure $\boldsymbol {\boldsymbol {\phi }}_k(\boldsymbol {x})$ operated by the structure $\boldsymbol {\boldsymbol {\phi }}_j(\boldsymbol {x})$ is more intense when the modal structure $\boldsymbol {\boldsymbol {\phi }}_j(\boldsymbol {x})$ describes large-scale flow features.

Figure 4.

Magnitude of the average interaction tensor coefficients ${\mathsf{N}}_{ijk}$ for three POD modes across the spectrum, $i = 1, 10$ and 75 in panels (a), (b) and (c) respectively, for a model resolving 99 % of the fluctuation kinetic energy. Panel (d) shows the coefficient $\chi _i(n)$ as a function of the normalised cutoff $n$ for the same three modes.

We now consider energy analysis of a large, full-resolution $\mathrm {DFT}_f$ model constructed from five partitions of thirty time units as discussed in § 3.1. The model is composed of all $N=300$ modes, corresponding to $150$ distinct frequencies. We perform modal decomposition and energy analysis on each partition separately, and then average the mean energy transfer rate tensor ${\boldsymbol{\mathsf{N}}}$ over the five partitions. Figure 5(a) shows the mean transfer rate distribution for mode $i=100$. Energy interactions in the DFT model are very sparsely distributed on a thin horseshoe-shaped structure composed of three branches (denoted in the figure as L, C and U) of $2\times 2$ blocks, and all other mean energy transfer rates interactions are identically zero. This pattern results from the joint effect of the oscillatory nature of the temporal coefficients and the quadratic nonlinearity of system (2.3), which can only be satisfied by triads of modes having matching temporal wavenumbers. A less pronounced horseshoe-shaped distribution of the energy interactions has been previously observed in energy analysis of POD-based models of three-dimensional transitional boundary layers Rempfer & Fasel (Reference Rempfer and Fasel1994a,Reference Rempfer and Faselb). These authors noticed that low-energy modal structures resemble Fourier modes in the spanwise direction (justified by the spanwise periodic domain) and thus coefficients ${\mathsf{Q}}_{ijk}$ and energy interactions are non-zero only for specific triads of modes. In the present case, this pattern is determined exclusively by the temporal coefficients as the tensor ${\boldsymbol{\mathsf{Q}}}$ constructed from projection modes does not possess any sparsity structure and its coefficients have a similar statistical distribution to that obtained using the POD modes. This is illustrated in figure 6 showing maps of the first slice of the tensor ${\boldsymbol{\mathsf{Q}}}$ of the largest Galerkin models considered here, constructed from the POD and the DFT decompositions, in panels (a) and (b) respectively. We observe that no underlying structure is present except for the asymmetry already observed in the energy analysis in figure 4. This property is confirmed in the probability distribution of the coefficients, shown in panel (c).

Figure 5.

(a) Magnitude of the average interaction tensor ${\mathsf{N}}_{ijk}$ for $i=100$, with the three characteristics branches, showing that active interactions come in $2\times 2$ blocks corresponding to matching triads of modes. The small inset focuses on the interactions of branches $C$ and $U$. (b) Magnitude of the average interaction tensor (3.8) where the three branches of panel (a) have been unfolded on a larger plane spanned by the coordinates $l$ and $\eta$. The inset shows details of the interactions of the branch $U$ in the plane $\eta -l$.

Figure 6.

Maps of ${\mathsf{Q}}_{1jk}$ for Galerkin models constructed from the POD and the $\mathrm {DFT}_f$ decompositions, in panels (a) and (b), respectively. Panel (c) shows the probability distribution of all quadratic coefficients for these two models.

To facilitate the interpretation of the energy interaction pattern, we follow Rempfer & Fasel (Reference Rempfer and Fasel1994b) and Arbabi & Mezić (Reference Arbabi and Mezić2017) and define oscillatory modal structures

(3.6)\begin{equation} \boldsymbol{u}_l(t,\boldsymbol{x}) = a_{2l-1}(t){\boldsymbol \phi}_{2l-1}(\boldsymbol{x}) + a_{2l}(t){\boldsymbol \phi}_{2l}(\boldsymbol{x}), \end{equation}

numbered by the index $l$ and tracing fluid motion at a single frequency on a two-dimensional subspace. Their modal energy is

(3.7)\begin{equation} e_l(t) = \tfrac{1}{2}(a^{2}_{2l-1}(t) + a^{2}_{2l}(t)) + a_{2l-1}(t) a_{2l}(t) ({\boldsymbol \phi}_{2l-1}(\boldsymbol{x}),{\boldsymbol \phi}_{2l}(\boldsymbol{x})). \end{equation}

Numerical experiments show that, for large number of snapshots, pairs of modes ${\boldsymbol \phi }_{2l-1}(\boldsymbol {x})$ and ${\boldsymbol \phi }_{2l}(\boldsymbol {x})$ tend to be orthogonal. Hence, considering the evolution equation for the modal energy $e_l(t) \sim \frac {1}{2}(a^{2}_{2l-1}(t) + a^{2}_{2l}(t))$ leads to the condensed triadic interaction tensor $\hat {{\boldsymbol{\mathsf{N}}}}$ of size $(N/2, N/2, N/2)$ with entries

(3.8)\begin{equation} \hat{{\mathsf{N}}}_{lmn} = \sum_{i=l}^{l+1} \sum_{j=m}^{m+1} \sum_{k=n}^{n+1}{\mathsf{N}}_{ijk}, \end{equation}

lumping together the $2\times 2$ blocks of interactions at matching triads of figure 5(a). In addition, the three branches L, C and U can be unfolded and conveniently visualised on a two-dimensional plane spanned by the coordinate $l$, the modal structure index, and $\eta = m - n$, representing the distance in modal space between pairs of temporal wavenumbers. This unfolding process is shown in panel (b) of figure 5, and when repeated for all modal structures leads to the distribution shown in figure 7(a). In figure 7(b), we report the average transfer rate $\hat {{\mathsf{N}}}_{lmn}$ normalised with the total average transfer rate for each structure, the quantity

(3.9)\begin{equation} \hat{T}_l = \sum_{m=1}^{N/2} \sum_{n=1}^{N/2} \hat{{\mathsf{N}}}_{lmn}, \quad l = 1,\ldots, N/2, \end{equation}

to illustrate more clearly the relative strength of the interactions. In figure 7(c), the normalised mean transfer rate for $l=1$ and 150 is reported. Interactions between triads of pairs of DFT modes are organised in agreement with the physics of scale interactions previously discussed for POD models. In absolute terms, the most relevant interactions are clearly those located near the origin of the coordinates. These correspond to low-index modes where nonlinear interactions with other low-index modes dominate, while interactions with the high-index modes, for larger $\eta$, are less important. This suggests that a sparsification approach based on pruning the interactions involving the high-index modes, i.e. the small scales, would be effective. By contrast, for high-index modes, relevant energy interactions are organised in bands along the axes $m$ and $n$ and involve energy exchange between low-index modes and high-energy, high-index modes. This suggests that the dynamics of the small scales is driven primarily by non-local interactions with the largest structures of the flow and not by small-scale/small-scale quadratic interactions. The slight asymmetry visible in panel (c) arises from the structure of the coefficients tensor ${\mathsf{Q}}_{ijk}$ and has the same physical origin as that observed in figure 4 for the POD model.

Figure 7.

(a,b) Absolute and relative strength of the energy interactions between pairs of DFT modes for a model with $N=300$ visualised on the plane $m,n$, with the additional coordinates $l$ and $\eta$. (c) Relative energy interactions for the first and last mode pairs.

3.3. Sparsification of POD-based models

We now apply the methodology presented in § 2 to three POD-based models resolving $90\,\%$, $95\,\%$ and $99\,\%$ of the kinetic energy, respectively (see table 1 for details). Because the size of the database matrix ${\boldsymbol \varTheta }({\boldsymbol{\mathsf{A}}})$ grows quadratically with the number of modes, the number of possible interactions $q$ can easily become larger than the number of available snapshots $N_t$, resulting in an underdetermined regression problem and overfitting. This is a well-understood issue in data analysis and requires cross-validation techniques to ensure the statistical reliability of the result (Friedman et al. Reference Friedman, Hastie and Tibshirani2008). In this work, we employed $K$-fold cross-validation, using typically $K =10$. Briefly, the database is first divided into $K$ folds. The model is trained using $K -1$ blocks and the reconstruction error $\epsilon$ of (2.13) is obtained from the fold that was left out. This procedure is iterated over all folds, obtaining the mean and the standard deviation of $\epsilon$.

Table 1.

Normalised cumulative energy distribution $e(n)$ for POD and DFT modes, where the latter are sorted by energy content ($\mathrm {DFT}_e$) or by frequency ($\mathrm {DFT}_f$).

Figures 8(a)–8(c) show the sparsification curves on the $\rho$–$\epsilon$ plane for the three POD models considered. The mean of $\epsilon$ across the folds is displayed as a thick black line, while the grey dashed line indicates plus or minus one standard deviation. These curves have been obtained by solving problem (2.12) using strategy S1 and progressively increasing the regularisation weight. When low weights are used, dense systems with good prediction accuracy are obtained. The opposite is true for large weights, identifying points in the left part of the graph characterised by low density and poor prediction accuracy. As postulated in § 2, the curves show a sweet spot at around $\rho \approx 0.2$, displaying a plateau for $\rho \gtrsim 0.2$, while the error $\epsilon$ grows quickly when $\rho \lesssim 0.2$. These results indicate that it is possible to prune $\sim$80 % of the quadratic interactions in model (2.3) without influencing the average prediction accuracy. The red dashed line represents the reconstruction error obtained with a naive sparsification approach. The approach consists in pruning coefficients of ${\mathsf{Q}}_{ijk}$ in the area denoted as $\mathrm {HH}$ in figure 4(a). The approach exploits the structure of ${\mathsf{N}}_{ijk}$ and is therefore referred to as ‘greedy’. By varying $n \in (1, N)$, models with different sparsity are obtained, with $n=N$ corresponding to the original projection model projection (indicated as a red square in figure 8). The shapes of the sparsification curves for the greedy approach are similar to those obtained with the $l_1$ regression. This is a direct consequence of the existence of a subset of most relevant energy interactions. However, the reconstruction error obtained from the greedy method is generally higher than that obtained by the $l_1$ regression, since the optimisation procedure involved in the $l_1$ approach modulates the strength of the remaining interactions by tuning the active quadratic coefficients, minimising the prediction error. As we show later in § 3.5 dedicated to analysing the model performance in time integration, this difference will have a marked effect on the long-term temporal stability of the models.

Figure 8.

The $\rho$–$\epsilon$ curves for three POD models resolving 90 %, 95 % and 99 % of the kinetic energy in panels (a), (b) and (c), respectively. The black line represents the cross-validated error averaged over $K =10$ folds. The dashed grey lines represent plus/minus one standard deviation of the cross-validated error calculated over the folds. The red dashed line shows the reconstruction error obtained with the greedy approach. The squares indicate the global reconstruction error of the Galerkin model obtained directly from projection.

Regardless of the approach, the mean reconstruction error decreases as the resolved energy increases, moving from panel (a) to panel (c), as more modes participate in capturing the dynamics of velocity fluctuations. In addition, larger models can be more effectively sparsified, as the sparsification curve drops more rapidly. This results from the non-local structure of energy interactions shown in figure 4. When one additional low-energy mode is included, the number of relevant interactions to be retained in the model only grows as ${O}(N)$ and not as ${O}(N^{2})$, i.e. all non-local interactions with the rest of the hierarchy denoted as $\mathrm {LL}, \mathrm {LH}$ and $\mathrm {HL}$ in figure 4(a). Since the total number of possible interactions grows as ${O}(N^{3})$, larger models can be more effectively sparsified. This is conceptually in agreement with the observations of Taira et al. (Reference Taira, Nair and Brunton2016) on the sparsification properties of discrete vortex models. We also observe that the mean prediction error does not necessarily decrease monotonically when the density increases. This phenomenon is particularly visible for the model in panel (b) but all models reproduce the same behaviour. This is a symptom that the number of available snapshots (1500) is potentially not large enough for the number of coefficients ($q = N \times (N+1) / 2 + N+1 = 2926$ for the model in panel (c)) and overfitting would have occurred if no cross-validation had been performed.

3.4. Energy interactions identified by the regression and conservation properties

To visualise the sparsity pattern identified by the regression as the regularisation weight in (2.12) is increased (constant for all modes in strategy S1), we introduce the tensor $\boldsymbol {\gamma }$ with entries $\gamma _{ijk}$ defined as the value of the regularisation weight at which the corresponding coefficient ${\mathsf{Q}}_{ijk}$ is shrunk to zero by the LASSO. Figures 9(a)–9(c) show three slices of $\boldsymbol {\gamma }$ for modes $i=1, 10$ and $75$, respectively, for the largest POD model considered, capturing $99\,\%$ of the total fluctuation kinetic energy. The first interactions to disappear are the small-scale/small-scale interactions. Increasing the penalisation, interactions that are local in modal space are progressively pruned, leaving only non-local interactions involving triadic exchanges with the low-index modes for large penalisations. Interestingly, this pattern does not change qualitatively nor quantitatively as the modal index $i$ increases. In fact, a comparable number of interactions is retained across the hierarchy and the governing equations of all modes are sparsified by an equal amount. Hence, sparsification has not produced mode truncation, which would have occurred if all coefficients of some low-energy modes had been shrunk to zero by the LASSO. This behaviour can be justified by noting that the mean square acceleration $\| \dot {{\boldsymbol{\mathsf{A}}}}_i\|_2^{2}$ of the POD amplitude coefficients varies weakly with $i$. In fact, the sparsification pattern does not change significantly when strategy S2 is used. In figure 10(a–c) the base ten logarithm of the mean energy interaction tensor ${\mathsf{N}}_{ijk}^{s}$ computed as in (2.9) with the sparse coefficient tensor ${\boldsymbol{\mathsf{Q}}}^{s}$ is shown for the same three modal indices as in figure 9. The data refer to sparse models with $\rho =0.3$, located nearby the sweet spot of the curves in figure 8(c). It can be observed that the sparsified model has a pattern of interactions resembling that of the dense model in figure 4. However, weak interactions and the associated flow physics have been pruned. It is also clear that the asymmetry of the interaction pattern observed in figure 4 and the physical mechanism that originates it are invisible to the regression and the interaction pattern in figure 10 is now symmetric with respect to a swap of the indices $j, k$.

Figure 9.

Distribution of the base ten logarithm of $\gamma _{ijk}$ for $i=1,10$ and $75$, in panels (a), (b) and (c), respectively, for the POD model resolving 99 % of the fluctuation kinetic energy.

Figure 10.

Base ten logarithm of the sparsified interaction tensor ${\mathsf{N}}^{s}_{ijk}$ for $i=1,10$ and $75$ in panels (a), (b) and (c), respectively, for a POD model resolving $99\,\%$ of the fluctuation kinetic energy and $\rho = 0.3$.

Despite the aggressive pruning, the sparse models reproduce fairly accurately the overall structure of the intermodal energy budgets. In the present flow configuration, the convective nonlinearity is energy conserving and Galerkin models should obey the relation $\sum _{i=1}^{N} T_i = 0$ as $N\rightarrow \infty$, with $T_i = \sum _{j=1}^{N}\sum _{k=1}^{N} {\mathsf{N}}_{ijk}$ the time-averaged energy transfer rate to/from mode $i$. For finite-dimensional approximations, this property is not satisfied exactly and the residual of the summation can be taken as a measure of the overall energy balance. Figure 11(a) shows such residual for the $l_1$ sparsified models (empty circles) and for the models obtained with the greedy approach (empty squares), as a function of the density. The residual is normalised by the root mean square value of the rate of change of the integral fluctuation energy. Note that the greedy model at $\rho =1$ is the model obtained directly from projection. It can be observed that the energy conservation error is relatively small, of the order of $10^{-3} \div 10^{-4}$. For large densities, it is larger than that of the projection model, because the regression tunes model coefficients to minimise the mean square error on the modal accelerations and does not enforce this physical constraint directly. The energy conservation residual decreases for sparser models and is ten times smaller than the projection model, owing to the lower number of active coefficients that participate in the regression. Figure 11(b) shows the distribution of the time-averaged energy transfer rate associated with mode $i$ for the $l_1$ sparsified model at $\rho =0.3$ (red crosses) and the dense model obtained from projection (empty circles). Data are reported every second mode. For the projection model, the net energy transfer is negative for the first few modes and changes sign at $i \sim 10$. Physically, this trend suggests that the first few modes extract energy from the mean flow and feed the dissipative high-index modes via triadic interactions. The $l_1$ sparsified model correctly reproduces this global trend,even though no constraints have been introduced (Loiseau & Brunton Reference Loiseau and Brunton2018). As discussed in § 2, it is argued that this is a general properties of data-driven techniques relying on optimisation ideas, such as the LASSO, which naturally reproduce invariants and conservation properties embedded in the data to a level defined by noise levels. For instance, Taira et al. (Reference Taira, Nair and Brunton2016) used network-theoretic ideas to sparsify connections in a discrete vortex model and observed that sparsification conserves the invariants of discrete vortex dynamics.

Figure 11.

(a) The normalised energy conservation error as a function of the density $\rho$, for the greedy sparsification approach (empty squares) and the $l_1$ based approach (empty circles). (b) The total net energy transfer rate $T_i$ as a function of the modal index $i$ for two POD models resolving $99\,\%$ of the kinetic energy, with coefficients identified from projection and for a $l_1$ sparse model. One every two data points are reported.

3.5. Long-term temporal behaviour of the $l_1$ sparse systems

We now turn our attention to the long-time behaviour of the sparsified models under temporal integration. We consider results for models resolving $95\,\%$ of the turbulent kinetic energy as an illustrative example, and use the projections of the POD modes onto one of the DNS snapshots to obtain initial conditions. Results are shown in figure 12. Panel (a) shows the temporal evolution of the integral fluctuation kinetic energy $E(t)$, as defined in (2.7), for the dense projection model and the sparse models obtained with the $l_1$ and greedy approach, with same density $\rho = 0.3$. The temporal evolution obtained with these models is compared against the fluctuation kinetic energy from DNS. For the projection model, the integral kinetic energy grows substantially in the first 40 time units and then saturates on a value that is approximately two orders of magnitude larger than what observed in DNS. The over-prediction occurs because the truncation of the small scales in the ansatz (2.2) leads to a significant imbalance of the production–dissipation budget within the model (Noack, Papas & Monkewitz Reference Noack, Papas and Monkewitz2005; Noack et al. Reference Noack, Schlegel, Ahlborn, Mutschke, Morzyński, Comte and Tadmor2008, Reference Noack, Morzynski and Tadmor2011; Balajewicz et al. Reference Balajewicz, Dowell and Noack2013). A qualitatively similar behaviour, if not worse, is then necessarily observed for the model sparsified with the greedy approach, since neglecting weak interactions alone does not cure the original dissipation problems. Conversely, the $l_1$ model is able to predict the correct average fluctuation kinetic energy and has excellent long-term stability properties, despite this being not enforced in the regression procedure (see Fick et al. Reference Fick, Maday, Patera and Taddei2018). We argue that this is due to the fact that the $l_1$ procedure performs a ‘prune-then-calibrate’ approach, where weak interactions are first pruned and the remaining active coefficients are then tuned in the optimisation involved in (2.12) to match the reference dynamics. It is evident from these results that this second step is key to obtain accurate long-term behaviour. Figures 12(b)–12(d) show a shorter segment of state space trajectory from these models projected onto three different pairs of modes. We omit the data for greedy approach since orbits quickly drift out of the visible range and this case is thus qualitatively similar to that from the dense model. It is clear that the trajectory of the $l_1$ model remains in the same volume of state space occupied by the DNS projections, while the projection model drifts away to a different region of state space, over-predicting the integral kinetic energy. This is more effectively visualised in panel (e), displaying the average modal energy $\lambda _i = \overline {a_ia_i}$ as a function of the modal index. Data are obtained from averaging long trajectories after the initial transient in panel (a) has completed. The projection and greedy models predict much larger energies across the entire spectrum, while the $l_1$ correctly predicts the correct decay of modal energies. It is, of course, not possible to guarantee that $l_1$ sparsified models of generic turbulent flows will have good long-term stability (Schlegel & Noack Reference Schlegel and Noack2015), but the present results constitute evidence that this is realistically possible on a non-trivial problem. Finally, we have observed in animations of the reconstructed flow fields using the spatial modes and the temporal coefficients from numerical integration (see supplementary movie) that characteristic flow features, such as the erratic burst of corner vortices and the evolution of coherent structures in the shear layer bounding the main vortex, are also well reproduced by the $l_1$ model, providing a realistic flow reconstruction over long time horizons.

Figure 12.

(a) Temporal evolution of the turbulent kinetic energy $E(t)$ from DNS compared to that obtained from temporal integration of the $l_1$, greedy and Galerkin projection models. Sparse models have $\rho = 0.3$. (b)–(d) State space projections onto three different mode pairs. Data for the greedy model are omitted as the trajectory quickly leaves the visible range. (e) Average modal energy ($\lambda _i = \overline {a_ia_i}$) from DNS and long-time integration of the models considered in panel (a).

3.6. Sparsification set-up for DFT-based models

Before moving to the sparsification of DFT-based models, we briefly discuss three technicalities arising from the oscillatory nature of the temporal coefficients.As an illustrative example, we consider a small-sized model constructed with $N=26$ modes and perform sparsification as discussed in § 2, with a relatively small regularisation weight ($\gamma = 10^{-14}$, strategy S1). The first key result is that all the entries of the constant and quadratic coefficient tensors ${\boldsymbol{\mathsf{C}}}$ and ${\boldsymbol{\mathsf{Q}}}$ are set to zero, while the linear tensor ${\boldsymbol{\mathsf{L}}}$ has a characteristic bidiagonal structure, shown in figure 13(a). The system identified by the regression is equivalent to a set of $N/2$ decoupled linear oscillators in the form

(3.10)\begin{equation} \begin{bmatrix} \dot{a}_{2l-1} \\ \dot{a}_{2l} \end{bmatrix} = \omega_l \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} a_{2l-1} \\ a_{2l} \end{bmatrix} \quad l = 1,\ldots, N/2, \end{equation}

coupling pairs of temporal coefficients oscillating at the same angular frequency $\omega _l = 2{\rm \pi} l / T$, with $T$ being the observation time. The eigendecomposition of the tensor ${\boldsymbol{\mathsf{L}}}$ is trivial. Eigenvalues are all imaginary and come in pairs that are integer multiples of the fundamental frequency $\omega _1 = 2 {\rm \pi}/T$. While this result is consistent with recent ideas on Koopman operator theory (Mezić Reference Mezić2013), where nonlinear dynamics is modelled with a linear system of larger dimension, all information on nonlinear energetic interactions has been lost in the process since the nonlinear part of the system has been completely eliminated by the regression. This result is due to the fact that, when temporal coefficients are sine/cosine pairs, there is a column of ${\boldsymbol \varTheta }({\boldsymbol{\mathsf{A}}})$ that is exactly parallel to the target $\dot {{\boldsymbol{\mathsf{A}}}}_i$, since time differentiation is equivalent to a permutation of sine/cosine pairs. As pointed out in Brunton et al. (Reference Brunton, Noack and Koumoutsakos2019) incorporating and enforcing known flow physics is a challenge and opportunity for machine learning algorithms. In order to address this first aspect, we introduce a physically motivated approach based on considerations of the time averaged energy budget of system (2.7). Since the temporal coefficients have zero mean and are uncorrelated in time, we should obtain

(3.11)\begin{equation} {\mathsf{L}}_{ii}\overline{a_i a_i} + \sum_{j=1}^{N} \sum_{k=1}^{N} {\mathsf{Q}}_{ijk}\overline{a_ia_ja_k} = 0, \quad i = 1, \ldots, N \end{equation}

i.e. only the diagonal element of the linear term participates in the mean power budget. Hence, for the sparsification of DFT-based models we use a modified database matrix that only contains the column associated with the diagonal part of the linear term.

Figure 13.

(a) Entries of the linear tensor ${\boldsymbol{\mathsf{L}}}$ identified by the unconstrained regression. (b) Singular values $\sigma _i$ of the full database matrix $\boldsymbol {\varTheta }({\boldsymbol{\mathsf{A}}})$ of (2.10) (red circles), and of the reduced matrix (black crosses) obtained by keeping only the subset of columns corresponding to active interactions on the three branches of figure 5. One every five singular values are shown for clarity.

The second aspect is that for DFT models the database matrix $\boldsymbol {\varTheta }({\boldsymbol{\mathsf{A}}})$ is not full rank and some of the columns of this matrix are linearly dependent. In this case, the LASSO is known to select one column at random (according to the particular ordering of the columns) and sets to zero regression coefficients of the other linearly dependent columns (Tibshirani Reference Tibshirani2013; Hastie, Tibshirani & Wainwright Reference Hastie, Tibshirani and Wainwright2015). Machine learning techniques often come without guarantees for robustness (Brunton et al. Reference Brunton, Noack and Koumoutsakos2019), implying that physical insight obtained with these tools might be questionable. To avoid this problem, we constructed a reduced database matrix ${\boldsymbol \varTheta }({\boldsymbol{\mathsf{A}}})$ containing only columns corresponding to the interactions on the three branches of figure 5. The reduced database matrix is full rank, as can be seen in figure 13(b), showing the singular values of the full database matrix defined by (2.10) and of the reduced matrix. The important consequence is that the solution of the LASSO problem (2.12) is unique (Tibshirani Reference Tibshirani1996), and can be thus compared with the available physical knowledge of scale interactions in turbulent flows. This reduction is not strictly necessary, since the $l_1$ regression identifies this pattern anyway for fairly small regularisation weights. However, this has the advantage that the computational complexity of sparsifying the entire Galerkin model only grows as ${O}(N^{2})$ instead of ${O}(N^{3})$, as for POD models, because the reduced database matrix contains a number of interactions equal to $q = 2(N+1)$ at most. More importantly, because of the greatly reduced number of free coefficients cross-validation techniques to avoid over-fitting become unnecessary.

The third aspect of DFT-based models is that, as anticipated, the number of modes is not uniquely defined by the energy resolution but depends on the overall observation time. Long observation times would be beneficial to reach statistical significance but would result in low-energy/low-frequency modes that do not contribute significantly to the overall dynamics. In practice, we have divided the original dataset into $M$ partitions and performed DFT for each of them separately. Then, we stacked vertically the modal acceleration matrices and the reduced database matrices from the partitions and solved (2.12) for a common coefficients vector.

3.7. Sparsification of DFT-based models

We now move to the sparsification of DFT-based models. We focus primarily on the structure of energy interactions identified by the regression and leave long-term temporal stability considerations aside. In fact, the DFT produces modal structures by assuming a priori their temporal behaviour, i.e. harmonic motion, and the meaning of a time-domain analysis is thus conceptually unclear.

We introduce the modified density $\rho _{DFT}$ spanning the range $[0, 1]$ and representing the number of active coefficients with respect to the total number of active interactions on the three branches of figure 5. For large models, the approximation $\rho _{DFT} \approx 2/3\rho N$ can be used. In figure 14(a), sparsification curves for three models obtained with observation times $T=10$, $30$ and $50$ (with $M=15$, $5$ and $3$ partitions of the full dataset, respectively), at full energy resolution, are reported. Strategy S1, where the regularisation weight is maintained constant for all modes is used. We observe that the error decreases monotonically with the observation time. This is a consequence of the larger number of frequencies that interact quadratically to reconstruct the original DNS acceleration data. For the larger model obtained at $T=50$, $70\,\%$ of the triadic interactions can be pruned with no major effects on the overall prediction error. If the full coefficient tensor ${\boldsymbol{\mathsf{Q}}}$ is considered, this correspond to a remarkably low density of 0.0015. Figure 14(b) shows the sparsification curves for models obtained with observation time $T = 30$, for three different energy resolutions, $e(n) = 0.9$, $0.95$ and $0.99$. Interestingly, we notice that the curves do not present a plateau for high densities as opposed to the full resolution mode shown in panel (a) and the POD sparsification curves of figure 8. This is the combined effect of the dramatic decrease in the number of modes at lower energy resolutions (see table 1) and the inherent efficient description of energy interactions in DFT-based models compared to POD.

Figure 14.

(a) Sparsification curves for models obtained by three different observation times and resolving $100\,\%$ of the kinetic energy. (b) Sparsification performed with $T=30$ with three different energy resolutions $e(n)$.

We now compare strategies S1 and S2 on the full resolution model obtained with observation time $T=30$. Results of this analysis are reported in figure 15. Panels (a–c) and (d–f) are obtained with the strategy S1 and S2, respectively. Panel (a) shows the tensor $\hat {\boldsymbol {\gamma }}$, obtained by processing and visualising the full tensor $\boldsymbol {\gamma }$ using the same technique utilised for the interaction tensor $\hat {{\boldsymbol{\mathsf{N}}}}$ in figure 5. Panel (b) shows the density of individual ordinary differential equations for a selected number of modal structures as a function of the overall model density $\rho _{DFT}$, while the sparsified interaction tensor $\hat {{\boldsymbol{\mathsf{N}}}}^{s}$ for $\rho _{DFT} = 0.7$ is shown in panel (c). When the regularisation weight is maintained constant, the sparsification pattern emerging from the tensor $\hat {\boldsymbol {\gamma }}$ follows the distribution of the mean energy transfer rate of figure 7(a). In particular, despite the signature of non-locality is still visible in the pattern, the sparsification is highly skewed across the spectrum because the equations for high-frequency modes are excessively sparsified for moderate penalisations as opposed to those of low-frequency, high-energy modes.

Figure 15.

(a)–(c) Strategy S1; (d)–(f) strategy S2. Panels (a,d) show the distribution of $\hat {\boldsymbol {\gamma }}$. Panels (b,e) show the trend of the modal density $\rho _l$ against the global density $\rho _{DFT}$ for four different modes in different parts of the spectrum. Panels (c,f) show the energy interaction tensor ${\boldsymbol{\mathsf{N}}}^{s}$ of the sparsified system.

This behaviour is better seen in the individual density curves in panel (b). Specifically, the density $\rho _l$ of the last mode pair ($l=150$) drops quite pronouncedly to much lower density than average at $\rho _{DFT} \approx 0.5$. Panel (d) shows the sparsification pattern obtained with the second strategy. We observe that, in this case, the interactions are retained according to their relative strength producing a sparsification pattern that follows the relative energy transfer rate reported in figure 7(b). This results in a more balanced sparsification across the spectrum, where the modal density $\rho _l$ decreases more uniformly for all modes as the global density is decreased, as shown in panel (e). The mean energy transfer rate of the models sparsified using the two strategies, with $\rho _{DFT} = 0.7$, is reported in panels (c,f). Globally, the structure and intensity of energy interactions is preserved by the LASSO, although strategy S1 has more aggressively sparsified the high-index modes and truncated the equations of the last five pairs of modes. Although not shown here, all DFT models, regardless the strategy, have similar energy conservation properties as for POD modes, as illustrated in figure 11.

As a final remark, we have observed in sparsification of larger DFT models that, although the LASSO is able to successfully identify the dominant subset of energy interactions, the complexity of the optimisation problem makes an accurate reconstruction of the numerical values of the system coefficients challenging. This is due to the spectral properties of the database matrix $\boldsymbol {\varTheta }$ which deteriorate as the number of modes considered grows (Cordier et al. Reference Cordier, El Majd and Favier2010). A potential solution to this issue would be to use elastic-net regression (Friedman et al. Reference Friedman, Hastie and Tibshirani2008) which combines an $l_1$ term with an $l_2$ (Tikhonov) penalisation. This would provide a better trade-off between sparsification and stability of the reconstructed coefficients.