2.1. POD

2.1.1. General formulation

Proper orthogonal decomposition is arguably the most popular tool for representation and analysis of turbulent flow fields (Lumley Reference Lumley1967).

The POD method allows us to derive an orthogonal basis for the (sub)space of the fluctuations of a multidimensional quantity ${f}$ of finite variance. One can show that a basis for the space of fluctuations, defined as ${f}'(t) := {f}(t) - \langle\,{f}\rangle$, with $\langle \cdot \rangle$ the statistical mean, is given by the set of elements $\left \{\boldsymbol {\phi }_{{n}}\right \}_{{n}}$, eigenvectors of the following eigenvalue problem (Holmes et al. Reference Holmes, Lumley and Berkooz1996):

(2.1)\begin{equation} C \boldsymbol{\phi}_{{n}} = \lambda_{{n}} \boldsymbol{\phi}_{{n}}, \end{equation}

with $\lambda _{{n}}$ the eigenvalue and $C \in \mathbb {R}^{{N_x} \times {N_x}}$ the empirical two-point covariance matrix,

(2.2)\begin{equation} C = \frac{1}{{{N_s}}} \sum_{{i}=1}^{{{N_s}}}{{f}'\left(t_{i}\right) {f}'\left(t_{i}\right)}, \end{equation}

with $\left \{t_{i}\right \}_{i}$ the time instants for which the field ${f}$ is available. Some conditions on the temporal sampling scheme apply for the empirical covariance $\widehat {C}$ to be an accurate approximation of $C$ (Holmes et al. Reference Holmes, Lumley and Berkooz1996). Proper orthogonal decomposition modes are identified as the eigenvectors $\boldsymbol {\phi }_{{n}}$.

2.1.2. Method of snapshots

The above method is a quite natural implementation of the underlying Hilbert–Schmidt decomposition theory. However, the algorithmic complexity associated with the eigenvalue problem (2.1) scales as $O({{N_s}} \, {N_x}^2)$, where the number of field instances ${{N_s}}$ was assumed to be lower than the size ${N_x}$ of the discrete field, ${{N_s}} \le {N_x}$. For large field vectors (large ${N_x}$), the computational and memory cost is hence high. For this widely encountered situation, a possible workaround was suggested in Sirovich (Reference Sirovich1987) and consists of solving the following eigenvalue problem:

(2.3)\begin{equation} \tilde{C} \, \boldsymbol{a}_{{n}} = \lambda_{{n}} \boldsymbol{a}_{{n}}, \quad \boldsymbol{a}_{{n}} \in \mathbb{R}^{{N_s}}, \end{equation}

with

(2.4)\begin{equation} \tilde{C} _{{i}, {{i}'}} \propto \left\langle\,{f}'\left(t_{i}\right), {f}'\left(t_{{i}'}\right)\right\rangle_\varOmega, \quad \forall \ {i}, {{i}'} \in \left[1, {{N_s}}\right] \subset \mathbb{N} \end{equation}

and $\langle \cdot , \cdot \rangle _\varOmega$ the Euclidean inner product. Since the correlation matrix $\tilde {C}$ is Hermitian, its eigenvalues are real and non-negative, $\lambda _{{n}} \ge 0$, $\forall \, {n}$, and its eigenvectors $\left \{\boldsymbol {a}_{{n}}\right \}_{{n}}$ are orthogonal and can be made orthonormal in an Euclidean sense, $\boldsymbol {a}_{{n}}^{\mathrm {T}} \, \boldsymbol {a}_{{n}'} \propto \delta _{{n}, {n}'}$, with $\delta$ the Kronecker delta. The spatial POD modes are finally retrieved via projection as follows:

(2.5)\begin{equation} \boldsymbol{\phi}_{{n}} = \lambda_{{n}}^{{-}1/2} \, {F}' \, \boldsymbol{a}_{{n}}, \quad \forall \, {n}, \end{equation}

where the ${i}$th column of the matrix ${F}'$ is the snapshot ${f}'_{i}$.

The algorithmic complexity is now $O({{N_s}}^3)$ and scales much better than the standard POD approach ($O({{N_s}} \, {N_x}^2)$) in the usual situation where ${{N_s}} \ll {N_x}$. In this work, we rely on this so-called method of snapshots to implement POD.

Formally the decomposition of the snapshot matrix ${F}'$ is equivalent to a singular value decomposition (SVD)

(2.6)\begin{equation} {F}' = \varPhi \varSigma {A}^{\mathrm{T}}, \end{equation}

where $\varPhi$ is the matrix constituted by the ${n}$ columns $\boldsymbol {\phi }_{{n}}$.

Here, ${A}$ is the matrix containing the ${n}$ columns $\boldsymbol {a}_{{n}}$ and $\varSigma$ is a diagonal matrix whose entries are $\lambda _{{n}}^{-1/2}$. The snapshot matrix can thus be decomposed into a snapshot-mode matrix $A$ and into a cell-mode matrix $\varPhi$. The spatial modes or structures can be seen as latent variables allowing optimal reconstruction of the data in the $L_2$ norm or an equivalent. The decomposition can be truncated to retain only the ${N_T}$ largest values corresponding to the ${N_T}$ first columns of each matrix.

2.2. PLSA

In all that follows we will consider a collection of ${{N_s}}$ scalar fields $\{\,{f}_{i} \}_{{i}=1, \ldots , {{N_s}}}$. Each field is of dimension ${N_x}$ and consists of either positive or zero integer values on each grid cell. The key point is to interpret the value of the field ${f}_{i}$ measured on a grid cell ${l}$ of the snapshot ${i}$ as the activity level of the grid cell ${i}$, or in other words the number of times (occurrences) the grid cell ${i}$ can be considered as active. The representation of a snapshot is therefore equivalent to a list of cell indices ${i}$, each of which appears ${f}_{i}$ times.

This representation assumes that each snapshot ${f}_{i}$ consists of a mixture of structures $\boldsymbol {z}_{{n}}$.

Probabilistic latent semantic allocation adds a probabilistic flavour as follows:

1. (i) given a snapshot ${f}_{i}$, the structure $\boldsymbol {z}_{{n}}$ is present in that snapshot with probability $p(\boldsymbol {z}_{{n}}|\,{f}_{i})$;

2. (ii) given a structure $\boldsymbol {z}_{{n}}$, the grid cell $\boldsymbol {x}_{l}$ is reported as active (i.e. appears in the list) with probability $p(\boldsymbol {x}_{l}|\boldsymbol {z})$.

Formally, the joint probability of seeing a given snapshot ${f}_{i}$ and finding a grid cell $\boldsymbol {x}_{l}$ in the list of cells representing the snapshot is

(2.7)\begin{equation} p(\,{f}_{i},\boldsymbol{x}_{l})= p(\,{f}_{i})\sum_{{n}} p(\boldsymbol{z}_{{n}}|\,{f}_{i}) p(\boldsymbol{x}_{l}|\boldsymbol{z}_{{n}}), \end{equation}

where $p(\,{f}_{i})$, $p(\boldsymbol {z}_{{n}}|\,{f}_{i})$ and $p(\boldsymbol {x}_{l}|\,{f}_{i})$ are the parameters of the model – $p(\,{f}_{i})$ is the probability of obtaining such a snapshot ${f}_{i}$ and is constant in our case, $p(\,{f}_{i})=1/{{N_s}}$; $p(\boldsymbol {z}_{{n}}|\,{f}_{i})$ and $p(\boldsymbol {x}_{l}|\boldsymbol {z}_{{n}})$ can be inferred using the expectation–maximization (EM) algorithm of Dempster, Laird & Rubin (Reference Dempster, Laird and Rubin1977).

Using Bayes’ theorem, $p(\,{f}_{i}, \boldsymbol {x}_{l})$ can be equivalently written as

(2.8)\begin{equation} p(\,{f}_{i},\boldsymbol{x}_{l})= \sum_{{n}} p(\boldsymbol{z}_{{n}}) p(\boldsymbol{x}_{l}|\boldsymbol{z}_{{n}}) p(\,{f}_{i}|\boldsymbol{z}_{{n}}). \end{equation}

This alternative formulation shows a direct link between PLSA model and POD model (as mentioned above, POD is called latent semantic allocation, or LSA in text mining). The formulations in (2.6) and (2.8) are similar:

1. (i) the value of $f$ at cell ${n}$ for the snapshot ${i}$ represents the number of times cell ${n}$ will appear in the list of cells corresponding to snapshot ${i}$, and is therefore proportional to the probability of detecting cell ${n}$ in snapshot ${i}$;

2. (ii) the structure probability $p(\boldsymbol {z}_{{n}})$ is the equivalent of the diagonal matrix entry $\varSigma _{{n}}$;

3. (iii) the probability of the snapshot ${f}_{i}$ given the structure $\boldsymbol {z}_{{n}}$, $p(\,{f}_{i}| \boldsymbol {z}_{{n}})$ corresponds to the snapshot-mode matrix entry $A_{{i}, {n}}$;

4. (iv) the probability to detect the cell $\boldsymbol {x}_{l}$ as active given the structure $\boldsymbol {z}_{{n}}$, $p(\boldsymbol {x}_{l} | \boldsymbol {z}_{{n}})$, corresponds to the matrix entry $\varPhi _{{l}, {n}}$.

3.1. The probabilistic framework of LDA

Latent Dirichlet allocation extends PLSA to address its limitations. Its specificity is as follows.

1. (i) The introduction of a probabilistic model for the collection of snapshots: each snapshot is now characterized by a distribution over the structures which will be now called motifs.

2. (ii) The use of Dirichlet distributions to model both motif–cell and snapshot–motif distributions.

   The Dirichlet distribution is a multivariate probability distribution over the space of multinomial distributions. It is parameterized by a vector of positive-valued parameters $\boldsymbol {\alpha }=(\alpha _1,\ldots ,\alpha _N)$ as follows:

   (3.1)\begin{equation} p\left(x_1, \ldots, x_N ; \alpha_1, \ldots, \alpha_N\right) =\frac{1}{B(\boldsymbol{\alpha})} \prod_{{n}=1}^{N} x_{{n}}^{\alpha_{{n}}-1}, \end{equation}
   where $B$ is a normalizing constant, which can be expressed in terms of the Gamma function $\varGamma$,
   (3.2)\begin{equation} B(\boldsymbol{\alpha})= \frac{\prod_{{n}=1}^{N} \varGamma(\alpha_{{n}})}{\varGamma(\sum_{{n}=1}^{N} \alpha_{{n}})}. \end{equation}

The support of the Dirichlet distribution is the set of $N$-dimensional discrete distributions, which constitutes the $N-1$ simplex. Introduction of the Dirichlet distribution allows us to specify the prior belief about the snapshots. The Bayesian learning problem is now to estimate $p(\boldsymbol {z}_{{n}}, {f}_{i})$ and $p(\boldsymbol {x}_{l}, \boldsymbol {z}_{{n}})$ from ${F}$ given our prior belief $\boldsymbol {\alpha }$, and it can be shown that Dirichlet distributions offer a tractable, well-posed solution to this problem (Blei et al. Reference Blei, Ng and Jordan2003).

3.2. Principles of LDA

Latent Dirichlet allocation is therefore based on the following representation.

1. (i) Each motif $\boldsymbol {z}_{{n}}$ is associated with a multinomial distribution $\boldsymbol {\varphi }_{{n}}$ over the grid cells ($p(\boldsymbol {x}_{l}|\boldsymbol {z}_{{n}})= \varphi _{{l}, {n}}$). This distribution is modelled with a Dirichlet prior parameterized with an ${N_x}$-dimensional vector ${\boldsymbol {\beta }}$. The components $\beta _{l}$ of ${\boldsymbol {\beta }}$ control the sparsity of the distribution: values of $\beta _{l}$ larger than unity correspond to evenly dense distributions, while values lower than unity correspond to sparse distributions. In all that follows, we will assume a non-informative prior, meaning that ${\boldsymbol {\beta }} = \beta \boldsymbol {1}_{N_x}$.

2. (ii) Each snapshot ${f}_{i}$, is associated with a distribution of motifs $\boldsymbol {\theta }_{i}$ such that $\theta _{{n},{i}} = p(\boldsymbol {z}_{{n}}|\,{f}_i)$. The probabilities of each motif sum to unity in each snapshot. This distribution is modelled with an ${N_T}$-dimensional Dirichlet distribution of parameter $\boldsymbol {\alpha }$. The magnitude of $\boldsymbol {\alpha }$ characterizes the sparsity of the distribution (low values of $\alpha _{{n}}$ correspond to snapshots with relatively few motifs). The same assumption of a non-informative prior leads us to assume $\boldsymbol {\alpha }=\alpha \boldsymbol {1}_{N_T}$.

3.3. Snapshot generation using LDA

The generative process performed by LDA with ${N_T}$ motifs is the following.

1. (i) For each motif $\boldsymbol {z}_{{n}}$, a cell–motif distribution $\boldsymbol {\varphi }_{{n}}$ is drawn from the Dirichlet distribution of parameter $\beta$.

2. (ii) For each snapshot ${f}_{i}$:

   1. (a) a snapshot–motif distribution $\boldsymbol {\theta }_{{i}}$ is drawn;

   2. (b) each intensity unit $1 \le {j} \le N_{i}$ where $N_{i}$ is the total intensity with $N_{i}= \sum _{{l}} f_{{l}, {i}}$ is then distributed among the different cells as follows. Firstly, a motif $\boldsymbol {z}_{{n}}$ is first selected from $\boldsymbol {\theta }_{i}$ (motif $\boldsymbol {z}_{{n}}$ occurs with probability $\theta _{{n}, {i}}$ in the snapshot). Secondly, for this motif, a cell ${l}$ is chosen among the cells using $\varphi _{{l}, {n}}$ and the intensity associated with cell ${l}$ is incremented by one.

The generative process can be summarized as follows.

Algorithm 1

LDA Generative Model.

The snapshot–motif distribution $\boldsymbol {\theta }_{{i}}$ and the cell–motif distribution $\boldsymbol {\varphi }_{{n}}$ are determined from the observed ${f}_{i}$. They are, respectively, ${N_T}$- and ${N_x}$-dimensional categorical distributions. Finding the distributions $\boldsymbol {\theta }_{{i}}$ and $\boldsymbol {\varphi }_{{n}}$ that are most compatible with the observations is an inference problem that can be solved by either a variational formulation (Blei et al. Reference Blei, Ng and Jordan2003) or a Markov chain Monte Carlo (known as MCMC) algorithm such as a Gibbs sampler (Griffiths & Steyvers Reference Griffiths and Steyvers2002). In the variational approach, the objective function to minimize is the Kullback–Leibler divergence. The solution a priori depends on the number of motifs and on the values of the Dirichlet parameters $\alpha$ and $\beta$.

From a computational perspective, both Markov chain Monte Carlo and variational methods require the performance of several passes over each snapshot and cell to estimate the motif distribution, yielding a computational complexity of the order of (${N_x} {{N_s}} {N_T}$). In practice, stochastic approximations are used to speed up estimation, as is done in the implementation used for our experiments. It allows us to process snapshots one at a time or in small batches, to perform stochastic gradient steps, and achieves convergence much faster. Latent Dirichlet allocation has been shown to scale to millions of snapshots (${{N_s}} > 10^6$), with tens of thousands of cells, and hundreds to thousands of motifs (Hofmann et al. Reference Hofmann, Blei, Wang and Paisley2013).

3.4. A simple example

We now give a short illustration of LDA for a one-dimensional example. We consider five basis functions represented in figure 1(a). These basis functions can be interpreted as motifs $\boldsymbol {\varphi }_{{n}}$. Using the generating algorithm described in the previous section, we created 5000 snapshots as random mixtures of the motifs. Each snapshot was obtained from 3000 individual draws (the number of draws was kept constant for the sake of simplicity). For each draw, we first sampled a topic ${n}$ from a Dirichlet distribution of parameter 0.2, and we then determine a cell by sampling the corresponding distribution $\boldsymbol {\varphi }_{{n}}$. Examples of snapshots can be seen in figure 2. Application of LDA to the collection of snapshots for a number of five motifs yields the motifs represented in figure 1(b), which can be seen to be a relatively good approximation of the five original ones. It can be seen that the motifs are very different from POD eigenmodes, represented in figure 3.

Figure 1.

(a) Original motifs; (b) motifs recovered through application of LDA.

Figure 2.

Selected one-dimensional snapshots.

Figure 3.

Proper orthogonal decomposition modes extracted from the collection.

3.5. Remarks

We conclude this section with two remarks.

1. (i) Latent Dirichlet allocation can generalize to new snapshots more easily than PLSA, due to the snapshot–motif distribution formalism. In PLSA, the snapshot probability is a fixed point in the dataset, which cannot be estimated directly if it is missing.

   In LDA, the dataset serves as training data for the Dirichlet distribution of snapshot–motif distributions. If a snapshot is missing, it can easily be sampled from the Dirichlet distribution instead.

2. (ii) An alternative viewpoint can also be adopted in interpreting the LDA in the form of a regularized matrix factorization (MF) method. This is further discussed in Appendix A.

4.1. Numerical configuration

The idea of this paper is to apply this methodology to snapshots of turbulent flows in order to determine latent motifs from observations of $Q_-$ events. We will consider the configuration of turbulent channel flow at a moderate Reynolds number of $R_{\tau }= u_{\tau } h/ \nu = 590$ (Moser, Kim & Mansour Reference Moser, Kim and Mansour1999; Muralidhar et al. Reference Muralidhar, Podvin, Mathelin and Fraigneau2019), where $R_{\tau }$ is the Reynolds number based on the fluid viscosity $\nu$, channel half-height $h$ and friction velocity $u_{\tau }$. Wall units based on the friction velocity and fluid viscosity will be denoted with a subscript $_+$. The streamwise, wall-normal and spanwise directions will be referred to as $x,y$ and $z$, respectively. The horizontal dimensions of the numerical domain are $(2 {\rm \pi}, {\rm \pi})h$. Periodic boundary conditions are used in the horizontal directions. The resolution of $(256)^3$ points is based on a uniform spacing in the horizontal directions, which corresponds to streamwise and spanwise spacings of $\varDelta _{x+}=7.24$ and $\varDelta _{z+}=3.62$, and a hyperbolic tangent stretching function for the vertical direction. More details can be found in Podvin & Fraigneau (Reference Podvin and Fraigneau2017). The configuration is shown in figure 4. More details about the numerical simulation can be found in Muralidhar et al. (Reference Muralidhar, Podvin, Mathelin and Fraigneau2019).

Figure 4.

Numerical domain $D$. The shaded surfaces correspond to the two types of planes used in the analysis. The volume considered for 3-D analysis is indicated in bold lines.

4.2. LDA inputs

In this section, we introduce the different parameters of the study. The python library scikit-learn (Pedregosa et al. Reference Pedregosa2011) was used to implement LDA. The sensitivity of the results to these parameters will be examined in a subsequent section.

We first focus on two-dimensional (2-D) vertical subsections of the domain, then present 3-D results. The vertical extent of the domain of investigation was the half-channel height. Since this is an exploration into a new technique, a limited range of scales was considered in the horizontal dimensions: the spanwise dimension of the domain was limited to 450 wall units. The streamwise extent of the domain was in the range of 450–900 wall units. The number of realizations considered for 2-D analysis was ${{N_s}}=800$, with a time separation of 60 wall time units. The number of snapshots was increased to 2400 for 3-D analysis.

The scalar field ${f}$ of interest corresponds to $Q_-$ events. It is defined as the positive part of the product $-u'v'$ , where fluctuations are defined with respect to an average taken over all snapshots and horizontal planes. The LDA procedure requires that the input field consists of integer values: it was therefore rescaled and digitized and the scalar field $f$ was defined as

(4.1)\begin{equation} f = [A \tau_{-}], \end{equation}

where $\tau _-= \max ( -u'v',0 )$ and $[\cdot ]$ represents the integer part. The rescaling factor $A$ was chosen in order to yield a sufficiently large, yet still tractable, total intensity. In practice we used $A=40$, which led to a total intensity $\sum _{i} \sum _{{l}} f_{{l},{i}}$ of approximately $10^8$ for plane sections and corresponds to 1516 discrete intensity levels. The effect of the rescaling factor will be examined in a subsequent section.

Latent Dirichlet allocation is characterized by a user-defined number of motifs ${N_T}$, a parameter $\alpha$ which characterizes the sparsity of prior Dirichlet snapshot–motif distribution, and a parameter $\beta$ which characterizes the sparsity of the prior Dirichlet motif–cell distribution. Results were obtained assuming uniform priors for $\alpha$ and $\beta$ with a default value of $1/{N_T}$. The sensitivity of the results to the priors will be evaluated in § 5.2.

4.3. LDA outputs

For a collection of ${{N_s}}$ snapshots and a user-defined number of motifs ${N_T}$, LDA returns ${N_T}$ motif–cell distributions $\boldsymbol {\varphi }_{{n}}$ and ${{N_s}}$ snapshot–motif distributions $\boldsymbol {\theta }_{i}$. Each motif is defined by a probability distribution $\boldsymbol {\varphi }_{{n}}$ associated with each grid cell. It is therefore analogous to a structure or a portion of structure since it contains spatial information – note, however, that its definition is different from standard approaches. The motif–snapshot distribution $\boldsymbol {\theta }_{i}$ characterizes the prevalence of a given motif in the snapshot.

As will be made clear below, the motifs most often consist of single connected regions, although occasionally a couple of distinct regions were identified. In most cases, the motifs can thus be characterized by a characteristic location $\boldsymbol {x}^c$ and a characteristic dimension in each direction $L_j$, $j \in \{x,y,z\}$.

To determine these characteristics, we first define for each motif a mask associated with a domain $D$. The origin of the domain was defined as the position, $\boldsymbol {x}_{m}$ corresponding to its maximum probability $p_{m} = \boldsymbol {\varphi }_{{n}} (\boldsymbol {x}_{m})$. The dimensions of the domain in each direction (for instance $L_x$) were defined as the segment extending from the domain origin over which the probability remained larger than $1\,\%$ of its maximum value $p_{m}$. The position and characteristic dimension of a motif for instance in the $x$-direction are then defined as

(4.2)\begin{gather} x^c = \frac{\int_{D} x \boldsymbol{\varphi}_{{n}} {\mathrm{d}} D }{\int_{D} \boldsymbol{\varphi}_{{n}} {\mathrm{d}} D }, \end{gather}

(4.3)\begin{gather}L_x^2 = 4 \frac{\int_{D} (x - x^c)^2 \boldsymbol{\varphi}_{{n}} {\mathrm{d}} D }{\int_{D} \boldsymbol{\varphi}_{{n}} {\mathrm{d}} D }. \end{gather}

Analogous definitions can be given for $y^c$ and $z^c$.

6.1. Reconstruction

We now examine how the flow can be reconstructed using LDA. In all that follows, without loss of generality, we will focus on one of the cross-flow planes examined in § 5, specifically the cross-section at $x=0$ of dimensions $d_{y+}=590$ and $d_{z+}=450$. As described in the algorithm presented in § 3, both the motif–snapshot and the cell–motif distributions can be sampled for the total intensity $N_{i} = \sum _{l} f_{{l},{i}}$ in the ${i}$th snapshot. This total intensity is defined as the rescaled integral value of the Reynolds stress (digitized and restricted to $Q_-$ events) over the plane. Since results were found to be essentially independent of the rescaling, we can make the simplifying assumption that $N_{i}$ is large enough so that the distribution $\varphi _{{n}}$ is well approximated by the samples. For a given total intensity $N_{i}$, a reconstruction of the ${i}$th snapshot can then be obtained at each grid cell $\boldsymbol {x}_{l}$ from

(6.1)\begin{equation} \tau^{R\text{-}LDA}(\boldsymbol{x}, t_{i}) = \frac{1}{A} f_i(\boldsymbol{x}) = \frac{N_{i}}{A} \sum_{{n}=1}^{{N_T}} \theta_{{n}, {i}} \varphi_{{n}}(\boldsymbol{x}), \end{equation}

where

1. (i) $\varphi _{{n}}(\boldsymbol {x})$ is the motif–cell distribution;

2. (ii) the snapshot–motif distribution $\theta _{{n}, {i}}$ represents the likelihood of motif $\boldsymbol {z}_{{n}}$ in the ${i}$th snapshot.

It seems natural to compare this reconstruction with the POD representation of the field which has a similar expression,

(6.2)\begin{equation} \tau^{R\text{-}POD}(\boldsymbol{x},t_{i})= \sum_{{n}=0}^{N_{POD}-1} a_{{n},{i}} \phi_{{n}}(\boldsymbol{x}), \end{equation}

where

1. (i) $\phi _{{n}}(\boldsymbol {x})$ are the POD eigenfunctions extracted from the autocorrelation tensor $C_{{i}, {i}'}$ obtained from the ${{N_s}}$ snapshots;

2. (ii) $a_{{n},{i}}$ corresponds to the amplitude of the ${n}$th POD mode in the ${i}$th snapshot.

The first six fluctuating POD modes are represented in figure 17. We note that the zeroth POD mode represents the temporal average of the field. As expected, the fluctuating POD modes consist of Fourier modes in that spanwise direction (due to homogeneity of the statistics), and their intensity reaches a maximum at around $y_+ \simeq 25$.

Figure 17.

Contour plot of the first six fluctuating normalized POD spatial modes; contour values go from $-0.03$ to $0.03$. Negative values are indicated by dashed lines. $(a)\ n = 1, \ (b)\ n = 2, \ (c)\ n = 3, \ (d)\ n = 4, (e)\ n = 5, \ (\,f)\ n = 6. $

If the number of POD modes is equal to the number of motifs ${N_T}$, by construction, POD will provide a better representation of the statistics at least up to second order (Holmes et al. Reference Holmes, Lumley and Berkooz1996). We note that, in terms of computational requirement, POD is less expensive than LDA, as it requires solving an SVD problem versus implementing an iterative EM algorithm (Dempster et al. Reference Dempster, Laird and Rubin1977). However, the performance of the EM algorithm can be improved, in particular with online updates (Hofmann Reference Hofmann1999).

In terms of storage, a reconstructed snapshot requires $N_{POD}$ modes for POD and ${N_T}$ topics for LDA. However, storage reduction could be obtained in the case of LDA by filtering out the motifs with a low probability $\theta _{{n}, {i}}$, i.e. lower than a threshold $\kappa$. We note that, in this case, it is necessary to store the indices ${n}$ of the motifs as well as the value of $\theta _{{n}, {i}}$, so that if $n$ modes (respectively, topics) are kept, storage will consist of $2 n$ variables per snapshot. We see that storage reduction can be achieved if the fraction of retained modes $\eta =n/{N_T}$ is sufficiently small. The LDA storage data length per snapshot $2 \eta {N_T}$ should then be compared with the POD data length $N_{POD}$.

For ${N_T}=96$, choosing a threshold of $\kappa =0.015$ resulted in less than 8 % difference between the filtered and unfiltered LDA reconstructions (the $L_2$ norm was used). The average value for $\eta$ was $0.2$, which means that the number of POD modes that would represent a storage equivalent to that of LDA with ${N_T}=96$ is $N_{POD} \simeq 2 \eta {N_T} \simeq 40$. We note that the total storage cost should further take into account the size of the LDA basis $\{\boldsymbol {z}_{{n}}\}_{{n}}$, which will be larger than the POD basis $\{\boldsymbol {\phi }_{{n}}\}_{{n}}$ since they are, respectively, equivalent to ${N_T}$ and $N_{POD}$ fields. However, efficient storage of the LDA basis can be achieved by making use of the limited spatial support of $\boldsymbol {z}_{{n}}$, in particular for motifs located close to the wall.

In the remainder of this section we will compare a filtered LDA reconstruction of 96 motifs (where values of $\theta _{{n}, {i}}$ lower than $\kappa =0.015$ are excluded from the reconstruction), with a POD representation of $N_{POD}=48$ modes, which captures approximately 75 % of the total energy. Figure 18 compares an instantaneous field with its LDA reconstruction and its POD reconstruction. A more general assessment is provided by figure 19, which shows the correlation coefficient between each snapshot and its reconstruction based on POD as well as that based on LDA. Although POD appears to be slightly superior, as can be expected from its optimality property, the correlation coefficients are very close with respective average values of 0.75 for LDA and 0.77 for POD.

Figure 18.

Instantaneous Reynolds stress field (limited to $Q_-$ events) (a) True field; (b) LDA-reconstructed field using 96 modes; (c) POD-reconstructed field using 48 modes.

Figure 19.

Distribution of the correlation coefficient between each original snapshot and its reconstruction based on (a) LDA or (b) POD.

6.2. Generation

Latent Dirichlet allocation is a generative model, so it is straightforward to generate synthetic snapshots by sampling from distributions $\theta$ and $\varphi$ for a total intensity $N_{{i}}=\sum _{{l}} f_{{l}, {i}}$, which is modelled as a Poisson process with the same mean and standard deviation as the original database.

In contrast, POD is not a generative model per se. We will use a simplified version of the probabilistic extension of POD (PPCA) derived by Tipping & Bishop (Reference Tipping and Bishop1999), which is presented in Appendix B, where we will make the additional assumption that no noise is present in the model, POD-based synthetic fields will be reconstructed from deterministic spatial POD modes $\boldsymbol {\phi }_{{n}}$ and random POD amplitudes $\boldsymbol {a}_{{n}}$ which are assumed to be Gaussian variables. Examination of figure 20, which represents the distribution of the first fluctuating POD coefficients ${n} \ge 1$, suggests that it is quite acceptable as a first approximation to assume Gaussian distributions for the amplitudes $\boldsymbol {a}_{{n}}$ – alternatively, the amplitudes could be sampled from the empirical distributions. The amplitude of the zeroth mode, which corresponds to the average of the field over the snapshots, will be assumed to be constant for all snapshots.

Figure 20.

Histograms of the normalized amplitudes of the first six fluctuating POD modes and comparison with a sampled Gaussian distribution. $(a)\ n = 1, \ (b)\ n = 2, \ (c)\ n = 3, \ (d)\ n = 4, \ (e)\ n = 5, \ (\,f)\ n = 6.$

We can therefore compare the databases reconstructed from and generated with LDA with those obtained from POD. The generated databases consist of ${{N_s}}$ snapshots corresponding to arbitrary instants $\tilde {t}_{i}$. Overall, the statistics of five different databases can be compared:

1. (i) the true database $\tau _-(y,z,t_{i})$ corresponding to the actual values of the $Q_-$ events;

2. (ii) the POD-reconstructed (R-POD) or POD-projected database,

   (6.3)\begin{equation} \tau_-^{R\text{-}POD}(y,z,t_{i}) = \sum_{{n}=0}^{N_{POD}-1} a_{{i},{n}} \phi_{{n}}(y,z), \end{equation}
   where $\phi _{{n}}$ are the POD eigenfunctions and $a_{{i},{n}}$ are the amplitudes of the ${n}$th POD mode in the ${i}$th snapshot;

3. (iii) the POD-generated (G-POD) database,

   (6.4)\begin{equation} \tau_-^{G\text{-}POD}(y,z,\tilde{t}_{i}) = \sum_{{n}=0}^{N_{POD}-1} \tilde{a}_{{i},{n}} \phi_{{n}}(y,z), \end{equation}
   where $\tilde {a}_{{i},0}=\langle a _{{i},0}\rangle$, with $\langle \cdot \rangle$ the average over all snapshots and $\tilde {a}_{{i},{n}}$, ${n} \geq 1$, centred Gaussian random variables with variance $\langle \tilde {a}_{{i},{n}}^2\rangle$;

4. (iv) the LDA-reconstructed database (R-LDA),

   (6.5)\begin{equation} \tau_-^{R\text{-}LDA}(y,z,t_{i}) = \frac{N_{i}}{A} \sum_{{n}=1}^{{N_T}} \theta_{{n}, {i}} \varphi_{{n}}(y,z), \end{equation}
   where $N_{i}$ is the total intensity measured in the ${i}$th snapshot, $\theta _{{n}, {i}}$ is the distribution of motif ${n}$ on the ${i}$th snapshot and $\varphi _{{n}}(y,z)$ is the identified distribution of the cell at $(y,z)$ on motif ${n}$.

5. (v) the LDA-generated database (G-LDA),

   (6.6)\begin{equation} \tau_-^{G\text{-}LDA}(y,z,\tilde{t}_{i}) = \frac{\tilde{N}_{i}}{A} \sum_{{n}=1}^{{N_T}} \tilde{\theta}_{{n},{i}} \varphi_{{n}}(y,z), \end{equation}
   where $\tilde {N}_{i}$ is the total intensity, which is sampled from a Poisson process, $\varphi _{{n}}(y,z)$ is the identified distribution of the cell at $(y,z)$ on motif ${n}$ and $\tilde {\theta }_{{n},{i}}$ is sampled for each ${n}$ from the empirical distribution $\theta _{{n}, {i}}$ over the snapshots.

Figure 21 shows the statistics of the different databases as a function of the wall distance. Averages are taken over all snapshots and in the streamwise direction. The mean value of the Reynolds stresses is correctly recovered by all methods. The second-order statistics are slightly better recovered by the POD-reconstructed and POD-generated snapshot sets, but both LDA approaches also capture a significant portion of the variance. The POD databases capture 75 % of the total variance, while the reconstructed and generated LDA databases capture 68 % and 60 % of the variance, respectively. Again, the superiority of POD over LDA is expected since it is optimal in the L2 sense, however, it should be noted that the gain is relatively small. Figure 22 shows the vertical spatial autocorrelation of $\tau _-$ defined as $R(y, y')=\langle \tau _-(x,y,z,t)\tau _(x,y',z,t)\rangle$ (where $\langle \cdot \rangle$ represents an average taken in time and in the spanwise position). We can see that the generated LDA autocorrelation is very similar to its reconstructed POD counterpart, which shows that the LDA synthetic fields capture as much as the spatial structure as the POD reconstructed ones. We note that the autocorrelation at large separations is well reproduced by all datasets.

Figure 21.

Statistics of the different databases averaged over the spanwise direction and the number of snapshots: (a) mean value; (b) standard deviation.

Figure 22.

Spatial autocorrelation of the Reynolds stress (limited to $Q_-$ events) in the vertical direction at different heights. The average is taken over snapshots and in the spanwise direction. $(a)\ y_{0+} = 33, (b)\ y_{0+} = 163.$

Figure 23 shows histograms of the fields at different heights. We note that unlike the LDA approach, which is a non-negative decomposition (since it is based on probabilities), some negative values are observed for the POD approach, even though the original field values considered are always positive. We can see that at different wall distances the POD-reconstructed database reproduces well the distribution of the original database, but the POD-generated database does not. This failure is due to the fact that although POD amplitudes are uncorrelated by construction, they are not independent. We note that the same failure was observed when sampling the POD coefficients from their data-observed distributions $a_{{i}, {n}}$ instead of Gaussian processes. In contrast, both reconstructed and generated LDA methods yield very similar distributions, which reproduce the main features of the original Reynolds stress values, such as the intermittency (sharp peak at zero) and the asymptotic decay for positive values. The LDA is therefore able to reproduce strongly non-Gaussian statistics, which is not possible with its POD equivalent.

Figure 23.

Histograms of the Reynolds stress (limited to $Q_-$ events) corresponding to the different databases at different heights. $(a,\!b)\ y_{+} = 19, \ (c,\!d)\ y_{+} = 61, \ (e,\!f)\ y_{+} = 157, \ (g,\!h)\ y_{+} = 343.$