An evolve-then-correct reduced order model for hidden fluid dynamics

In this paper, we put forth an evolve-then-correct reduced order modeling approach that combines intrusive and nonintrusive models to take hidden physical processes into account. Specifically, we split the underlying dynamics into known and unknown components. In the known part, we first utilize an intrusive Galerkin method projected on a set of basis functions obtained by proper orthogonal decomposition. We then formulate a recurrent neural network emulator based on the assumption that observed data is a manifestation of all relevant processes. We further enhance our approach by using an orthonormality conforming basis interpolation approach on a Grassmannian manifold to address off-design conditions. The proposed framework is illustrated here with the application of two-dimensional co-rotating vortex simulations under modeling uncertainty. The results demonstrate highly accurate predictions underlining the effectiveness of the evolve-then-correct approach toward realtime simulations, where the full process model is not known a priori.


I. INTRODUCTION
In fluid mechanics community, flow systems are often characterized by excessively large spatio-temporal scales. This puts a severe restriction on efficient deployment of practical applications which usually require near realtime, and many-query responses. Traditional full-order simulations cannot achieve this since their computational cost cannot be handled even with the largest supercomputers. Therefore, reduced order modeling emerges as a natural choice to tackle the computationally expensive problems with acceptable accuracy [1][2][3][4][5][6][7][8]. In reduced order models (ROMs), the evolution of the most important and relevant features is tracked rather than simulating each point in the flow field. These features represent the dynamics of the underlying flow patterns, also called basis functions or modes.
Proper orthogonal decomposition (POD) is a common technique to extract the dominant modes that contribute most to the total system's energy [9]. POD, coupled with Galerkin projection (GP), has been used for many years to formulate ROMs for dynamical systems [10][11][12][13]. In these ROMs, the full-order set of equations is projected on the reduced space resulting in a dynamical system (in modal coefficients) with much lower order than the full order model (FOM). Clearly, this type of ROMs is limited to systems for which we have access to the governing equations. That is why they are called intrusive ROMs. However, in many situations, there is a discrepancy between the governing equations and the observed system. This might result from the approximation of underlying phenomena, incorrect parameterization, or insufficient information about source terms. This is, in particular, evident in atmospheric and geophysical sys- * osan@okstate.edu tems where many complex phenomena and source terms interact with each other in different ways [14][15][16].
Conversely, pure data-driven techniques, also known as nonintrusive, solely depend on observed data to model the system of interest. Therefore, more complicated processes can be modeled without the need to formulate them with mathematical expressions [17]. Machine learning (ML) tools have shown substantial success in the fluid community identifying the underlying structures and mimicking their dynamics [18][19][20][21][22][23]. However, end-to-end modeling with ML, especially deep learning, has been facing stiff opposition, both in academia and industry, because of their black-box nature, lack of interpretability and generalizability which might produce nonphysical results [24][25][26]. A perspective on machine learning for advancing fluid mechanics is available in a recent review article [20].
Hybridization of both of the aforementioned approaches is therefore sought to maximize their pros and mitigate their cons [27]. Several research efforts have been devoted to achieving this hybridization, for example, in the form of closure modeling [28][29][30][31][32], accelerating simulations [33], and enforcing physical laws by tailoring loss functions [34][35][36]. In this paper, we address the problem of hidden physics, or unknown source terms by utilizing a neural network with the assumption that observed data is a manifestation of all interacting mechanisms and sources. In particular, a long-short term memory (LSTM) architecture, a variant of recurrent neural networks, is used to account for the unknown physics by learning a correction term representing the discrepancy between the physical model and actual observations. This is the part where intrusive ROM fails. Meanwhile, the generalizability of the model under different operating conditions is retained by employing intrusive (physics-based) ROM to represent the core dynamics. In other words, the underlying physics is divided into two parts, known part (core physics) modeled by intrusive ap-proaches (e.g., Galerkin ROM) and unknown part (hidden physics) modeled by nonintrusive approaches (e.g., LSTM). A similar framework was presented in another study where we proposed a modular hybrid analysis and modeling (HAM) approach to account for hidden physics in parameterized systems [37]. The main difference between the two studies is the mechanism of adding a correction to the Galerkin projection model. In the previous study, the GP model is corrected dynamically with the LSTM model at each timestep and the corrected modal coefficients are used to predict the future state recursively. In the present study, we propose an evolve-thencorrect (ETC) approach to account for the hidden physics where the GP and LSTM models are segregated. First, the GP model is used to evolve the initial state to the final time based on the known core dynamics. Then, a correction from the LSTM model is added to the time instant of interest. This means that the GP model evolves with uncorrected modal coefficients and the correction is enforced statically as a final post-processing step.
Relevant to our ETC approach, an evolve-then-filter approach has been introduced by Wells et al [38] as ROM regularization. In their study, GP ROM is evolved for one timestep, after which a spatial filter is applied to filter the intermediate solution obtained in the evolve step. This filtering reduces the numerical oscillation of the flow variables (i.e., adds numerical stabilization to ROM). More recently, Gunzburger et al [39] proposed an evolve-filterrelax approach for uncertainty quantification of the timedependent Navier-Stokes equations in convection dominated regimes. This is similar to the evolve-then-filter approach with the additional step of relaxation which averages the unfiltered and filtered flow variables to control the amount of numerical dissipation introduced by the filter. In the present paper, we put forth a nonintrusive memory embedding neural network model to take hidden physics into account for superparameterized systems. An illustration of the proposed framework is provided in solving fluid flow problems dynamically governed by the two-dimensional Navier-Stokes equations. We highlight that our ETC is modular in the sense that it does not require any change in the legacy codes and can be considered as an enabler for parametric model order reduction of partial differential equations.

II. EVOLVE-THEN-CORRECT APPROACH
In this study we consider the nonlinear dynamical system parameterized with a parameter µ which has the form u t (x, t; µ, κ) = F(x, t; u; µ) + Π(x, t; u; µ, κ), (1) where u is the state of the system, the subscript t denotes the temporal derivative, F is the dynamical core of the system governing the known processes parameterized by µ, and Π includes the unknown physics. The un-known physics encompasses deviation between the modeled and observed data resulting from several factors such as empirical parameterizations, and imperfect knowledge about the physical processes. The κ in our study refers to the control parameter modeling the interaction between dynamical core of the system and hidden physics.
We use the proper orthogonal decomposition (POD) to extract the dominant modes representing the above nonlinear dynamical system. We collect the data snapshots u 1 , u 2 , . . . , u N ∈ R M at different time instances. In POD, we construct a set of orthonormal basis functions that optimally describes the field variable of the system. We form the rectangular matrix A ∈ R M ×N consisting of the the data snapshot u n as its column. Then we use singular value decomposition (SVD) to compute the left and right singular vectors of the matrix A. In matrix form the SVD can be written as where W ∈ R M ×N , Σ ∈ R N ×N , and V ∈ R N ×N . The W and V contains the left and right singular vectors which are identical to the eigenvectors of AA T and A T A, respectively. Also, the square of singular vales are equal to the eigenvalues, i.e., λ k = σ 2 k . The vectors w k (also the eigenvectors of AA T ) are the POD basis functions and we denote them as φ k in this text. The POD basis functions are orthonormal (i.e., φ i , φ j = δ ij ) and are computed in the optimal manner in the L 2 sense [40,41]. The full order model (FOM) for the dynamical system can be approximated using these POD basis functions as follow, where R of the number of retained basis functions such that R << N and a k are the time dependent modal coefficients. The POD basis functions minimizes the meansquare error between the field variable and its truncated representation. Moreover, it also minimizes the number of basis functions required to describe the field variable for a given error [42]. The number of retained modes is usually decided based on their energy content. Using these retained modes we can form the POD basis set Φ = {φ k } R k=1 to build the ROM. POD is often complemented with Galerkin projection (GP) to reduce higher-dimensional partial differential equations (PDEs) into reduced-order ordinary differential equations (ODEs). To get GP equations, first we use a linear superposition of approximated field variable given by Equation 3 into the governing equation of the physical system. Then we apply inner product of the resulting equation with the orthonormal basis function φ k . Therefore, we need the complete information about the governing equation of the physical system to form the GP equations that can describe the system accurately. We do not know about the hidden physics given by the source term Π in Equation 1. Therefore, we cannot derive a fully intrusive GP model for such dynamical system. In this study, we use machine learning algorithm to model this source term due to hidden physics and GP equations are derived based on the dynamical core of the system F. The GP equations for the system with linear and nonlinear operator can be written aṡ or more explicitly, where L and N are the linear and nonlinear operator of the physical system. We limited our formulation in considering quadratic nonliterary without loss of generality by using R modes. We use third-order Adams-Bashforth (AB3) to numerically integrate Equation 5. In the discrete sence, the update formula can be written as where s and β q are the constants corresponding to AB3 scheme, which are s = 2, β 0 = 23/12, β 1 = −16/12, and β 2 = 5/12. We can obtain the true projection modal coefficients by simply projecting the field variable onto the basis functions and can be written as where the angle-parentheses refers to the Euclidean inner product defined as x, y = x T y = M i=1 x i y i . The true projection modal coefficients includes the hidden physics and its interaction with the dynamical core of the system. The GP modal coefficients given by Equation 5 does not model this effect. We can then define the correction as The above correction term can be learned with datadriven machine learning algorithm [28,43]. We employ long-short term memory architecture [44], a variant of recurrent neural network to learn this correction. LSTM are particularly suitable for time-series prediction, since they can use the information about the previous state of the system to predict the next state of the system. We train our LSTM neural network to learn the mapping from GP modal coefficients to the correction term, i.e., where C is the correction given by Equation 8. We use three lookbacks to be consistent with AB3 scheme during training (please see the details in [45]). Since, GP modal coefficients are used as input features to the LSTM network, the parameter µ governing the system's behavior is taken implicitly into account. Once the model is trained, we can correct the GP modal coefficients with LSTM based correction to approximate true projection modal coefficients.

ETC ROM = GP ROM + LSTM Correction
The ROM requires the basis functions which were obtained using POD. The dataset for the POD is generated using FOM. Running FOM for every ROM is opposite to the motivation behind using ROM. Hence, we should be able to generate the POD basis from existing POD basis sets using interpolation. The training data is generated for different values of parameters µ 1 , . . . , µ P . We compute the separate basis set Φ 1 , . . . , Φ P for each of these parameter. We employ Grassmann manifold interpolation [46,47] to compute the POD basis set for test parameter which is not included in the training set. The graphical illustration of the Grassmann manifold interpolation is provided in Figure 1 and the procedure is described in Algorithm 1.
where the trigonometric operators apply only to the diagonal elements. We demonstrate the performance of evolve-thencorrect (ETC) framework discussed in Section II for the two-dimensional Navier-Stokes equations. We use vorticity-streamfunction formulation of Navier-Strokes equations and it can be written as

III. NUMERICAL RESULTS
where ω is the vorticity defined as ω = ∇ × u, u = [u, v] T is the velocity vector, and ψ is the streamfunction. We use vortex merger as the test example. In this test example, a pair of co-rotating vortices are separated from each other by some distance. These vortices induce fluid motion and strongly interact in the merging process. Vortex merger is extensively studied in the two-dimensional context as it explains the average inverse energy cascades and the direct enstrophy cascade observed in twodimensional turbulence [48]. The initial condition for the vortex-merger test case is given as where their centers are initially located at (x 1 , y 1 ) = (3π/4, π) and (x 2 , y 2 ) = (5π/4, π). We utilize an arbitrary array of Taylor-Green vortices as the source term that represents a perturbation field (i.e., referring to hidden physics). The source term in Equation 14 is given below where F (t) = γ e −t/Re and η = 3. We use computational domain (x, y) ∈ [0, 2π] with periodic boundary conditions. We generate data snapshots for Re = [200, 400, 600, 800] with 256 2 spatial grid and a time-step of 0.01 from time t = 0 to t = 20. We test the ETC framework for out-of-sample condition at Re = 1000. The linear and nonlinear operators in GP equations for two-dimensional Navier-Stokes equation are where φ ω k and φ ψ k refer to POD basis functions of the vorticity and streamfunction fields, respectively [49]. We retain 8 basis functions (i.e., R = 8) as it captures more than 99.95% energy for all Reynolds number included in training. We illustrate ETC framework for two different amplitudes of the source term, γ = 0.01 and γ = 0.1.
We illustrate the ETC approach for two different magnitudes of the source term, γ = 0.01 and γ = 0.1. We use Re = 800 as the reference Reynolds number for both test cases. Figure 2 shows the true basis functions and Grassmann interpolated basis functions at Re = 1000. The Grassmann manifold interpolation procedure can compute correct basis functions with very little deviation (especially for basis φ 7 and φ 8 ) from true basis functions. We train the LSTM network with two hidden layers and 80 cells. Our experiments with different hyperparameters show that the results are not highly sensitive to the choice of hyperparameters. Figure 3 shows the evolution of vorticity modal coefficients for γ = 0.01. The GP modal coefficients are different from the true projection modal coefficients even though the 8 modes are capturing more than 99.95% of the energy. This difference is due to the source term not modeled by GP equations. In the ETC approach, we correct the GP modal coefficients and we get an excellent agreement with the true modal coefficients after correction. The vorticity modal coefficients obtained with the ETC approach is almost the same as the true projection modal coefficients. To show the difference between vorticity modal coefficients at all train Reynolds number and the test Reynolds number, we plot time series of modal coefficients for all Reynolds number in the same plot as reported in Figure 4. We can observe that the evolution of vorticity modal coefficients for the test Reynolds number is different from train Reynolds number and the ETC approach is able to produce the correct trajectory of modal coefficients at test Reynolds number. This also shows that the LSTM has not simply memorized the correction from train Reynolds number, but it has learned the dependence of modal coefficients on the parameter (through input features) governing the physical system.  Figure 5 displays the evolution of modal coefficients for γ = 0.1. We see that there is a large deviation between GP modal coefficients and true projection modal coefficients due to a large magnitude of the source term. The ETC approach can correct GP modal coefficients and produce the modal coefficients trajectories close to the true projection. Since the source term is very large and has the magnitude the same as the main field, we do not get the same level of accuracy as the test case with γ = 0.01, especially near the end time (i.e., t = 15 − 20). Figure 6 presents the vorticity modal coefficients at all train Reynolds number and the test Reynolds number.
Finally, Figure 7 shows the three-dimensional plot of the reconstructed vorticity field computed with FOM, true projection, GP, and ETC modal coefficients. We do not see a significant difference between the threedimensional vorticity field computed with GP, and ETC modal coefficients for γ = 0.01. However, the orientation of two vortices is not correctly captured by GP modal coefficients in comparison to the true projection. The ETC approach correctly predicts the orientation of two vortices close to the true projection. There is a large difference between the reconstructed vorticity field for GP and true projection at γ = 0.1. The magnitude of the reconstructed vorticity field with GP modal coefficients is very small compared to the true projection modal coefficients and it is not able to produce the vorticity field due to the Taylor-Green vortices. The ETC approach correctly predicts the vorticity field with sufficient accuracy and also produces the vorticity field due to Taylor-Green vortices. The correct prediction of the vorticity field for large source term illustrates the advantage of the ETC approach for physical systems where there is a discrepancy between the modeled and observed data due to modeling assumptions, imperfect parameterizations, and insufficient knowledge about the physical system.

IV. CONCLUDING REMARKS
We present an evolve-then-correct approach for reduced order modeling of parameterized systems with hidden information. This hidden information can be attributed to unknown source term, imperfect parameterizations, or incorrect modeling assumptions. We correct the physics-based model (GP model in this study) with the hidden information recovered by ML algorithms (LSTM neural network in our framework).
We demonstrate the performance of the ETC approach for two-dimensional Navier-Stokes equations used to simulate the merging of co-rotating vortices. Our numerical experiments with two different magnitudes of source term yield highly accurate solutions close to true projection results. We also show that the LSTM algorithm has effectively learned the relationship between the unknown source term and the parameter governing the system's behavior through GP modal coefficients as input features.
We generate the synthetic source term using an array of Taylor-Green vortices in our numerical experiments and the results of the ETC approach are promising. However, there are several research directions that we plan to pursue in the future to better understand the capability and limitations of the ETC approach. One of the future research directions is to test the ETC approach with more realistic noisy data obtained from sensors. Disclaimer. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.