Generative adversarial networks to infer velocity components in rotating turbulent flows

Abstract Inference problems for two-dimensional snapshots of rotating turbulent flows are studied. We perform a systematic quantitative benchmark of point-wise and statistical reconstruction capabilities of the linear Extended Proper Orthogonal Decomposition (EPOD) method, a nonlinear Convolutional Neural Network (CNN) and a Generative Adversarial Network (GAN). We attack the important task of inferring one velocity component out of the measurement of a second one, and two cases are studied: (I) both components lay in the plane orthogonal to the rotation axis and (II) one of the two is parallel to the rotation axis. We show that EPOD method works well only for the former case where both components are strongly correlated, while CNN and GAN always outperform EPOD both concerning point-wise and statistical reconstructions. For case (II), when the input and output data are weakly correlated, all methods fail to reconstruct faithfully the point-wise information. In this case, only GAN is able to reconstruct the field in a statistical sense. The analysis is performed using both standard validation tools based on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}L2 spatial distance between the prediction and the ground truth and more sophisticated multi-scale analysis using wavelet decomposition. Statistical validation is based on standard Jensen–Shannon divergence between the probability density functions, spectral properties and multi-scale flatness. Graphical abstract


Introduction
Understanding and predicting turbulent processes in geophysical and laboratory set-ups is key for many important questions [1][2][3][4].Although observation technologies are progressing, the quality and quantity of available data so far is still inadequate in many respects.The most important limitation are due to measurements sensitivity, spatial and temporal sparsity [5][6][7].As a result, different tools have been developed to reconstruct the whole 'dense' field from sparse/noisy/limited measurements.Gappy Proper Orthogonal Decomposition (GPOD) is based on a linear reconstruction in eigenmodes [8], and it was later improved and applied in [9] for the simple idealized and prototypical case of a flow past a circular cylinder with random spatio-temporal gappiness.More recently, inspired by the success of Convolutional Neural Networks (CNNs) in computer vision tasks, a series of proof-of-concepts studies have used a Generative Adversarial Network (GAN) [10] to reconstruct two-dimensional (2D) snapshots of three-dimensional (3D) rotating turbulence with a large gap [11,12].In our previous work [12], POD-and GAN-based reconstructions have been systematically compared in terms of the point-wise error and turbulent statistical properties, for gaps with different sizes and geometries.For super-resolution tasks, CNN and GAN have been applied to recover high-resolution laminar or turbulent flows from low-resolution coarse data in space and time [13][14][15][16].Moreover, CNN or GAN have been recently proposed to reconstruct the 3D velocity fields from several 2D sections [17] or a cross-plane of unpaired 2D velocity observations [18].
In this paper, we deal with different inference problems, connected to gappiness in the diversity of physical fields that can be measured.In particular, we ask how much one can make use of the measurement of one component of the 3D turbulent velocity field to predict another one.The problem is of course important for many field and laboratory applications where sensors/observations can collect fluctuations only in some preferred directions, e.g.Particle Image Velocimetry (PIV) [19,20] and wind observation from satellite infrared images [21][22][23][24][25]. Inference problems can be attacked by Extended Proper Orthogonal Decomposition (EPOD) using linear decomposition on a predefined basis of independent functions [26].Recently, for a turbulent open-channel flow, it has been shown that the 2D velocity-fluctuation fields at different wall-normal locations can be predicted with a good statistical accuracy from the wall-shear-stress components and the wall pressure with a fully non-linear CNN [27] and GAN [28].Here we ask a more complex multi-objective question, checking how much one can infer the unkown field in terms of (i) the point-to-point L 2 error and (ii) the statistical properties, e.g. on the basis of Jensen-Shannon (JS) divergence between the probability density function (PDF) of the ground truth and the one of the predicted velocity component.
In order to be specific with some important applications, we concentrate on the case of 3D turbulent flows under rotation, a system with plenty of physical interest and of high relevance Fig. 1 Energy spectra of different velocity components and the total energy spectrum from the TURB-Rot database [39].The Kolmogorov dissipative wave number, kη = 32, is determined as the scale where the total energy spectrum starts to decay exponentially.
to geophysical and engineering problems [29][30][31].Under strong rotation, the flow tends to become quasi-2D, with the formation of large-scale coherent vortical structures parallel to the rotation axis [32][33][34].Moreover, strong non-Gaussian, 3D and intermittent fluctuations exist at small scales because of the forward energy cascade [35][36][37][38].The tendency towards quasi-2D configurations implies that the out-plane velocity component, u z , behaves close to a passive scalar advected by the in-plane components, u x and u y .This can be further shown in Fig. 1 with the energy spectra of different velocity components from the database used in this study.It is obvious that the energy of u x or u y is dominant compared with that of u z over the whole range of scales, which indicates that u z is almost decoupled from the leading dynamics.Therefore, u x and u y are well correlated with each other while u z is expected to show less correlations (Fig. 2).In this work, we assess the potential of EPOD, CNN and GAN methods to infer the velocity components on 2D slices perpendicular to the rotation axis.We study two scenarios with different difficulties as shown in Fig. 2: (I) Using one in-plane component, u x , to predict the other, u y , and (II) using an in-plane component, u x , to predict the out-plane one, u z .
The organization of this paper is as follows.In Sect.2, we describe the databset used to evaluate the inference tools and introduce the EPOD and GAN-based methods.The CNN used in this study is taken from the generator part of the GAN.Results of above tools on the two inference tasks are given in Sect.3. Finally, we present the conclusions in Sect. 4.

Methodology 2.1 Dataset
To evaluate different tools on the inference of velocity components, we use the TURB-Rot database [39], obtained from Direct Numerical Simulations (DNS) of the Navier-Stokes equations for incompressible rotating fluid in a 3D periodic domain, which can be written as where u is the incompressible velocity, Ω = Ω ẑ is the rotation vector, p is the modified pressure in the rotating frame,the forcing f acts at scales around k f = 4 and it is generated from a second-order Ornstein-Uhlenbeck process [40,41].To enlarge the inertial range, the dissipation ν∆u is modeled by a hyperviscous term ν h ∇ 4 u.Moreover, a linear friction term β∆ −1 u is used at large scales to reduce the formation of a large-scale condensate [37].The Rossby number in the stationary regime is Ro = E 1/2 /(Ω /k f ) ≈ 0.1, where E is the flow kinetic energy.The Kolmogorov dissipative wave number, k η = 32, is chosen as the scale where the total energy spectrum E(k) starts to decay exponentially (Fig. 1).Given the domain length L 0 , the integral length scale is L = E/ kE(k) dk ≈ 0.15L 0 and the integral time scale is T = L/E 1/2 ≈ 0.185.Refer to [39] for more details of the simulation.
Although all inference tools in this study can be applied to 3D data, here we restrict to 2D horizontal slices in order to make contact with the geophysical observation, PIV experiments and to restrict the amount of data to be used for training (see [14,15,17,18] for a few applications to 3D datasets).To extract data from the simulation, we first sampled snapshots of the full 3D velocity field during the stationary regime.Snapshots are chosen with a large temporal interval ∆t s = 5.41T to decrease their correlations in time.There are 600 snapshots sampled over 3243T at early times for training and validation, while 160 snapshots sampled over 865T at later times.The time separation between the two samplings for training/validation and testing is more than 3459T .Second, the resolution of the sampled fields is downsized from 256 3 to 64 3 by a Galerkin truncation in Fourier space: where the cut-off wave number is chosen to be the Kolmogorov dissipative wave number, such as to eliminate only fully dissipative degrees of freedom where the flow can be approximated to be linear with a good approximation.This is needed to balance the request to reduce the amount of data to be analyzed, without reducing the complexity of the task.For each downsized configuration, we selected 16 x-y planes at different z levels and each plane is augmented to 11 (for training/validation) or 8 (for testing) different ones using the periodic boundary condition [39].After independent random shuffles of the two sets of planes, the dataset is in Train/Validation/Test split as follows, 84480/10560/20480, corresponding to 73.1%, 9.1% and 17.7% respectively.

EPOD inference
Let H V = C(Ω, V ) denote the space of continuous functions from the spatial 2D domain, Ω, to a vector space of physical quantities, V = S ⊕ G, where S and G are the disjoint subspaces composed by the known quantities (S) and those to be inferred (G).Specifically, define u V ∈ H V , u S ∈ H S and u G ∈ H G , then we have (see Fig. 2) for the inference task (I), and for the inference task (II), The first step of EPOD method is to compute the correlation matrix where with • we denote an average over the configurations selected for the training set.Then we solve the eigenvalue problem where n = 1, . . ., N S and N S is the number of physical points in the domain Ω, to obtain the eigenvalues σ n and the POD eigenmodes φ (n) S (x).It is easy to realize that any realization of the observed fields can be decomposed as where the POD coefficients are given by the scalar product: Furthermore, exploiting the orthogonality of the POD eigenmodes one can also write the following identity [26]: The idea of the Extended POD modes is to generalize the exact relation ( 9) also to the case where on the r.h.s.we use the quantities to be inferred, G when there are more than one components/quantities measured or to be inferred.
i.e. by replacing u S with u G but keeping using the same coefficients b n defined in (8): Once completed the training protocol and defined the set of EPOD modes (10) one can start the inferring procedure for any new configuration outside the training set, by taking the measured components u S (x), calculating the coefficients (8) and defining the predicted/inferred field as

GAN-based inference with context encoders
Context encoders embedding in GAN [42] have been used to reconstruct the missing data in square gaps and random gappiness of different sizes for the same turbulent database in our previous works [11,12].For the inference tasks of velocity components here studied, we use a GAN architecture consisting of a generator and a discriminator as shown in Fig. 3.The generator is a functional GEN (•) first encoding the supplied measurements, u S , to a latent feature representation, and second it generalizes by a decoder the latent representation to produce a candidate for the missing measurements, GEN (u S ) =u G is the missing true turbulent configuration.The discriminator plays as a 'referee' functional D(•) which takes either u G and outputs the probability that the supplied configuration is sampled from the true turbulent dataset.The loss function of the generator is where the L 2 loss is defined as the mean squared error (MSE), given the predicted and original data, u G and u G , averaged over the spatial domain of area A(I).The hyper-parameter λ adv is called the adversarial ratio and the adversarial loss is where p(u S ) is the probability distribution of the known input measurements in the training set and p p (u G ) is the probability distribution of the predicted fields from the generator.The discriminator is trained simultaneously with the generator to maximize the cross entropy between the ground truth and generated samples Here p t (u G ) is the probability distribution of the real -ground truth-fields, u G .It is possibile to show that the generative-adversarial training with λ adv = 1 in (12) minimizes the JS divergence between the generated and the true probability distributions, JSD(p t p p ) [10,43].Therefore, the adversarial loss helps the generator to produce predictions with correct turbulent statistics, while the generator, if trained alone (λ adv = 0), would minimize only the L 2 loss (12) which is mainly sensitive to the large energy-containing scales.Compared with the linear EPOD, GAN not only takes advantage of the non-linear expression of the generator, but also can optimize a loss considering both the MSE error and the generated probability distribution.More details of the GAN are provided in Appendix A, including the architecture, hyperparameters and the training schedule.

Results
In this section, we systematically compare EPOD, CNN and GAN-based methods for the two tasks of velocity component inference.The CNN has the same architecture as the generator of the GAN, thus it corresponds to our GAN with λ adv = 0 (without the discriminator); In this way we are also able to judge the importance of adding/removing the objective to have a small JS divergence with the ground truth for our inferred field.For the decision of the adversarial ratio for the GAN, readers can refer to Appendix A.
To quantify the inference error, we define the normalized mean squared error (MSE) as where is the spatially averaged L 2 error for one flow configuration and • represents now the average over the test data.The normalization factor is defined as where G is defined similarly.The form of E u G makes sure that a prediction with too small or too large energy gives a large MSE.
We use JS divergence to evaluate the PDF of the inferred velocity components.For two probability distributions, P (x) and Q(x), JS divergence measures their similarity and is defined as where M = 1 2 (P + Q) and is the Kullback-Leibler (KL) divergence.If the two probability distributions are close, it gives a small JS divergence, and vice versa.

Prediction at large scales and small scales
For the inference task (I), in Fig. 4 we present an inference experiment using a test data never showed during training.The aim is to use the velocity component u x (1st row, 1st column) to predict the velocity component u y , of which the ground truth and predictions from the different tools are shown in the other columns of the 1st row.In the 2nd row we show their gradients in x-direction, ∂u y /∂x.EPOD only keeps the correlated part with the given information and the prediction looks meaningful because of the correlation between u x and u y .With the capability of expressing non-linear correlations, CNN predicts better results which are close to the ground truth.However, CNN predictions are blurry without the high-frequency information, while GAN can generate realistic turbulent configurations with the benefit of adversarial training.
To analyze the predictions quantitatively, Fig. 5(top) shows the MSE(u y ) and the JS divergence between PDFs of the original and the predicted u y , which is denoted as JSD(u y ) = JSD(PDF(u ).We divide the test data into batches of size 128 (or 2048), calculate the MSE (or JS divergence) over these batches and we indicate with the error bound its range of fluctuation.The EPOD approach has a MSE(u y ) around 0.6, while CNN and GAN have smaller values around 0.1.Fig. 6 shows the PDF of the spatially averaged L 2 error, ∆ uy , over different configurations, the peaks of which give consistent results with MSEs.The JS divergence of the inferred component can be explained together with Fig. 6(b), which shows PDFs of the predicted velocity component compared with the original data.EPOD has the largest value of JSD(u y ) and the predicted PDF(u y ) has the correct shape but deviates from the original one at around |u y | = σ(u y ), where σ(u y ) represents the standard deviation of the original data.CNN and GAN generate nearly perfect PDFs with small values of JSD(u y ).Moreover, GAN is better with the help of the discriminator, which also slightly increases MSE(u y ).Fig. 5(bottom) shows the MSE and JS divergence between PDFs for the gradient of u y , ∂u y /∂x.The MSE(∂u y /∂x) behaves similarly as MSE(u y ) that EPOD gives a large value while CNN and GAN have smaller close values.However, for the gradient both EPOD and CNN give large values of JS divergence and only GAN predicts a small value, indicating the importance of adversarial training on the high-order quantities.

Multi-scale prediction error
Here, we evaluate the prediction from a multi-scale perspective with the help of wavelet analysis [44][45][46].For a field defined on a uniform grid of size 2 N × 2 N , 2D wavelet decomposition gives that where ūy is the mean value and is the wavelet contribution at wave number k j = 2 j , corresponding to the length scale 1/k j .Given that σ ∈ {x, y, d}, c j,ix,iy is the wavelet coefficient and ψ (x) j,kx,ky (x, y) = ψ j,kx (x)φ j,ky (y), ψ j,kx,ky (x, y) = φ j,kx (x)ψ j,ky (y), ψ (d) j,kx,ky (x, y) = ψ j,kx (x)ψ j,ky (y), (24) where φ(•) and ψ(•) are the Haar scaling function and associated wavelet, respectively.To measure the inference error at different scales, we define the normalized wavelet mean squared error (W-MSE) as which is the MSE between wavelet contributions at k j of the predicted and original data, following the definition (16).We stress that the normalization factor given by ( 18) and ( 19) is different for different k j .Fig. 7 displays the W-MSE obtained from different methods.It shows that at all scales EPOD has the largest prediction error.CNN and GAN produce close values of W-MSE, which is large at small scales where k j = 2 6 , while varies between 0 to 0.4 at 1 ≤ k j ≤ 2 5 .The wave number k j with the minimum W-MSE, i.e. the scale where the maximum percentage prediction is reached, corresponds to the scale in correspondence of the maximum of the spectrum, as shown by comparing with Fig. 8(a).

Spectra and Flatness
To further study the statistical properties of different predictions, in Fig. 8(a) we have measured the energy spectrum for the different inferred fields, namely where k = (k x , k y ) is the wave number, ûy (k) is the Fourier transform of u y (x) and û * y (k) is its complex conjugate.Since EPOD only extracts the correlated part with the supplied information, it predicts an energy spectrum with a similar shape but with smaller energy at all wave numbers compared with the original one.CNN benefits from the non-linear properties and predicts the correct energy spectrum at small wave numbers.With the help of the discriminator, GAN can predict close energy spectrum to the original one at all wave numbers.Fig. 8(b) plots the flatness where δ r u y = u y (x + r, y) − u y (x, y).The results are consistent with those of energy spectrum where EPOD fails at all scales, CNN predicts correct flatness at large scales and GAN has satisfying results at all scales.

Inference task (II)
Now we move to the inference task (II) which is more challenging as illustrated in the introduction.Fig. 9 presents an inference experiment similar as Fig. 4 but for the task (II).It is obvious that GAN predictions are realistic and well correlated with the ground truth.However, although the predicted structures are correct, GAN predictions can have wrong values which can be even with opposite signs in some vortices (5th column).The EPOD method is not able to infer meaningful results due to the complexity of the task.The

Conclusion
We studied a problem with practical applications in geophysical and engineering problems, which is using one velocity component to infer another one for 2D snapshots of rotating turbulent flows.Moreover, this problem is also with theoretical interest, connected to the feature ranking of turbulence: identifying which degrees of freedoms are dominant in turbulent flows.Linear (EPOD) and non-linear (CNN & GAN) methods are systematically compared by two inference tasks with different complexities, which depend on the correlation between the input component and the one to be inferred.The EPOD method conducts POD analysis with the supplied component and gives the correlated part as the prediction.CNN & GAN methods are based on fully non-linear approximators, without and with an adversarial component.
For the simpler inference task (I), EPOD generates meaningful inference which is not satisfying in terms of the MSE and the JS divergence with real data.An obvious improvement is reached after the usage of non-linear mapping, where CNN gives good predictions with small MSE and JS divergence with the ground truth.Compared with CNN, GAN makes the inference more realistic with fine structures, with the cost of slightly increasing the MSE.Besides, the training cost is also more expensive for GAN.For the inference task (II), the supplied component is not well correlated with the one to be inferred.The EPOD method cannot make meaningful inference because of this complexity.However, with the nonlinear capability CNN and GAN can recognize the positions of coherent structures, although they cannot predict correct values (or even signs) in the vortices.Moreover, the usage of discriminator is very important for this task.Indeed, CNN only predicts some smooth blobs while GAN generates turbulent configurations correlated with the ground truth.Finally, we also showed that GAN predictions always have good statistical properties.
Note that EPOD, CNN and GAN-based methods are straightforward to the inference task with multiple components/quantities supplied and/or to be inferred.The inference task where was also conducted and it shows very similar results with task (II), because of the high correlation between the two in-plane components.Finally, we stress that this study is for the case where the turbulent flow is fully resolved.Future study can focus on conducting inference with noisy and/or Fourier-filtered data, which is also important for practical applications such as PIV.
[23] J Sheng, E Malkiel, and J Katz.Using digital holographic microscopy for simultaneous measurements of 3d near wall velocity and wall shear stress in a turbulent boundary layer.Experiments in fluids, 45(6):1023-1035, 2008.

Fig. 2
Fig.2Examples of two inference tasks on 2D slices of 3D turbulent rotating flows: (I) Using ux to infer uy and (II) using ux to infer uz.The rotation axis is along zdirection.Note that the in-plane components, ux and uy, are strongly correlated among each other while the outplane component, uz, is less correlated with both ux and uy.

Fig. 3
Fig. 3 Architecture of generator and discriminator for the velocity inference tasks.Each convolution (up-convolution) layer has a kernel size of 4 followed by a Leaky Rectified Linear Unit (ReLU) activation function, except that the last layer of the generator uses a Tanh activation function.We can increase channels in the layer u S or the layers u (p) G and u (t)

Fig. 4
Fig. 4 Prediction visualization of the inference task (I) by the different tools for an instantaneous field.In the 1st row, the input velocity component, ux, is shown in the 1st column, while the 2nd to 5th columns show the ground truth and the inferred velocity component, uy, obtained from EPOD, CNN and GAN.The corresponding gradient in x-direction, ∂uy/∂x, is shown in the 2nd row.

Fig. 7 W
Fig. 7 W-MSE for the velocity component to be inferred, uy, at different wave numbers k j .Results are obtained from the different tools for the inference task (I).

Fig. 8
Fig.8(a) The energy spectrum and (b) the flatness of predictions and the ground truth of the velocity component, uy.The unit of r is the grid width, wg = 2π/64.Results are obtained from the different tools for the inference task (I).

Fig. 9
Fig.9Prediction visualization of the inference task (II) by the different tools for an instantaneous field.In the 1st row, the input velocity component, ux, is shown in the 1st column, while the 2nd to 5th columns show the ground truth and the inferred velocity component, uz, obtained from EPOD, CNN and GAN.The corresponding gradient in x-direction, ∂uz/∂x, is shown in the 2nd row.