# Full trajectory optimizing operator inference for reduced-order modeling
using differentiable programming

Surya Chakrabarti 111Ph.D. Student<EMAIL_ADDRESS>The Ohio State
University, Mechanical and Aerospace Engineering Department, Columbus, Ohio,
43210  Datta Gaitonde222John Glenn Professor<EMAIL_ADDRESS>The Ohio
State University, Mechanical and Aerospace Engineering Department, Columbus,
Ohio, 43210  Arvind T. Mohan333Scientist, Computational Physics and Methods
Group, Los Alamos National Laboratory, Los Alamos, New Mexico 87545  Daniel
Livescu444Scientist, Computational Physics and Methods Group, Los Alamos
National Laboratory, Los Alamos, New Mexico 87545

###### Abstract

Accurate and computationally inexpensive Reduced Order Models (ROMs) capable
of forecasting turbulent flows can facilitate rapid design iterations and thus
prove critical for predictive control strategies in engineering problems.
Galerkin projection based Reduced Order Models (GP-ROMs), derived by
projecting the Navier-Stokes equations on a truncated Proper Orthogonal
Decomposition (POD) basis, are popular because of their low computational
costs and theoretical foundations in the governing fluid flow equations.
However, the accuracy of traditional GP-ROMs degrades over long time
prediction horizons. While calibration can somewhat mitigate this issue, the
procedure is cumbersome and requires insight into the full order solver. To
address these issues, we extend the recently proposed Neural Galerkin
Projection (NeuralGP) data driven framework to compressibility-dominated
transonic flow, considering a prototypical problem of a buffeting NACA0012
airfoil governed by the full Navier-Stokes equations. The algorithm maintains
the form of the ROM-ODE obtained from the Galerkin projection; however
coefficients are learned directly from the data using gradient descent
facilitated by differentiable programming. This effectively blends the
strengths of the physics driven GP-ROM and purely data driven neural network-
based techniques, resulting in a computationally cheaper model that is easier
to interpret. While both NeuralGP and traditional GP-ROMs yield accurate
forecasts over short time prediction horizons, the former accurately tracks
the reference POD coefficients over longer prediction horizons, in time ranges
where the calibrated GP-ROM displays the tendency to diverge in amplitude and
frequency. We show that the distinct long-term behavior of the blended
NeuralGP method minimizes a more rigorous full trajectory error norm compared
to a linearized error definition optimized by the calibration procedure. We
also find that while both procedures stabilize the ROM by displacing the
eigenvalues of the linear dynamics matrix of the ROM-ODE to the complex left
half-plane, the NeuralGP algorithm adds more dissipation to the trailing POD
modes resulting in its better long-term performance. The results presented
highlight the superior accuracy of the NeuralGP technique compared to the
traditional calibrated GP-ROM method.

##  1 Introduction

Incorporating a robust, active feedback control for turbulent flows is a
promising research problem with the potential to greatly improve the design of
flight vehicles and more efficient energy systems [1, 2, 3, 4]. The goal of
controlling chaotic and nonlinear systems like turbulent flow is greatly
enabled by the capability to making predictions for several time instants into
the horizon in a fast, inexpensive and accurate manner. High fidelity
simulations of fluid flows, while accurate, can be prohibitively expensive for
efficient forecasting and prediction-based control techniques. This is
particularly the case for 3D high Reynolds number turbulent flows, where the
number of degrees of freedom can exceed billions, and requires substantial
resources on supercomputers.

The increased push towards developing control strategies for turbulent flows
has led to some focus on the development of Reduced Order Models (ROM) that
retain the accuracy of high fidelity simulations while reducing the
computational resources needed for predictions made over time horizons
spanning several non-dimensional flow times [5, 6]. One of the more popular
strategies uses Galerkin projection of the Navier-Stokes equations onto a set
of truncated orthogonal basis vectors, which are usually obtained from a
Proper Orthogonal Decomposition (POD) of the flow fluctuations. The resulting
projected equations may be solved at relatively little computational cost,
while capturing the dominant physics of the flow as discussed in Refs. [7, 5,
8, 9, 10]. Several of these exploit assumptions such as isentropic flow which
can constrain parameters where they may be used, such as moderate Mach
numbers. ROMs for compressible flow problems governed by the full Navier-
Stokes equations, while possible [11], have not been as widely explored.

Althought their theoretical formulation is directly predicated on the Navier
Stokes equations, traditional GP-ROMs suffer from structural instabilities [5,
12] that severely limit their prediction capabilities and result in model
divergence over longer time horizons. This has prompted the development of
several optimization techniques aimed at reducing prediction errors and
stabilizing the ROM solutions. Examples include non-linear constrained
optimization [13], least squares minimization [14] and Tikhonov regularization
[15]. Cordier et al. [15] compared several calibration methodologies and
concluded that the technique based on Tikhonov regularization outperformed all
other methods. Despite the application of state-of-the-art calibration
techniques however, long-term predictions from GP-based ROMs remain unreliable
[16], necessitating newer methods for improved prediction performance.
Furthermore, the construction of a GP-ROM is an intrusive process, requiring
knowledge of the governing equations and the discretization schemes employed
in the full order model [17]. Additionally, the introduction of new problems
with additional physics, such as reacting flows, requires a re-derivation of
the ROM Ordinary Differential Equation (ODE), further limiting the
generalizability of such methods.

Recent studies have leveraged deep learning-based approaches, such as neural
networks (NNs), to create data-driven alternatives to GP-ROMs with the goal of
improving long-term predictions from ROMs, while also avoiding the
intrusiveness of traditional ROM methodologies. Several variants of NNs exist,
with the Long Short-Term Memory (LSTM) architecture [18] being popular in
turbulence [19, 20, 21, 22] owing to its strengths in time-series prediction.
Although the literature reports that NNs have considerable promise, they face
important limitations to practical use because of their black-box nature, lack
of interpretability and the requirement of relatively large amounts of
training data. To avoid these problems associated with NNs for ROMs,
Peherstorfer and Willcox [17] proposed a non-intrusive operator inference
method which retains the form of the ROM-ODE obtained from the Galerkin
Projection, but seeks to estimate its coefficients directly using training
data. While a prior knowledge of the form of the ROM-ODE is required for
implementing this method, it can be easily obtained owing to the properties of
the Galerkin projection. For a full order model given by a Partial
Differential Equation (PDE) with a polynomial non-linearity, the Galerkin
projection onto a reduced basis preserves the structure of the equations and
results in a ROM-ODE with a similar form as the original full order PDE [23].

The operator inference method, however, relies on a linear optimization
technique that operates on training data comprising discretely sampled points
on the trajectories of the projected reduced state and as such, does not
account for the errors incurred over the full predicted trajectory obtained by
integrating the ROM-ODE. Recently, Mohan et al. [16] proposed a hybrid
learning approach called Neural Galerkin Projection (NeuralGP) that uses deep
learning to derive the coefficients for a ROM-ODE of a known form. The method
uses Differentiable Programming (henceforth referred to as DiffProg) - a
programming paradigm that enables software to perform full differentiation via
automatic differentiation [24, 25, 26]. This facilitates calculation of
gradients of the solutions from generalized ODE/PDE solvers with respect to
the ODE coefficients, allowing for their optimization using gradient-descent
based methods. We tested this algorithm for an isentropic flow, and while it
showed considerable promise yielding superior results compared to an
isentropic GP-ROM, the simplifying assumptions of constant entropy
significantly limits practical applications.

The present paper thus seeks to extend the NeuralGP method to a high-speed
compressible flow problem governed by the full Navier-Stokes equations and
compare the resulting ROM-ODEs and their predictive capability against the
more widely studied GP-ROMs. As a test case, we examine transonic buffet flow
over a NACA0012 airfoil, representing a well-examined aerospace application
often used as a test bed for feedback control [27, 28] and ROM [11] studies.
We summarize the numerical method used for the simulation in Sec. 2 while a
detailed description of the flowfield features are provided in Sec. 3. Two
different ROMs for the flowfield are then constructed using the GP-ROM and
NeuralGP, respectively. We outline a description of the procedure used to
construct the ROM along with key details of the mathematical derivation of
each method in Sec. 4. We evaluate the short and long-term predictions of the
NeuralGP in Sec. 5 and identify benefits relative to prior approaches In Sec.
6, we analyze the ROM-ODEs from a dynamical systems perspective, with the goal
of explaining factors that contribute to improvement in prediction
performance. Finally, we summarize our contributions in Sec. 7.

##  2 Numerical Method

The non-dimensionalized compressible Navier-Stokes equations cast in the
strong conservation form are solved in a curvilinear coordinate $(\xi,\eta)$
frame with $J$ being the Jacobian of the curvilinear transformation given by
$\bigg{(}J=\frac{\partial(\xi,\eta,\tau)}{\partial(x,y,t)}\bigg{)}$ [29].

$\frac{\partial}{\partial\tau}\bigg{(}\frac{\overrightarrow{Q}}{J}\bigg{)}+\frac{\partial\hat{F}}{\partial\xi}+\frac{\partial\hat{G}}{\partial\eta}=\frac{1}{Re}\bigg{[}\frac{\partial\hat{F_{\nu}}}{\partial\xi}+\frac{\partial\hat{G_{\nu}}}{\partial\eta}\bigg{]}$
(1)

where $\overrightarrow{Q}$ refers to a vector comprised of the conserved
variables. Thus, $\overrightarrow{Q}=[\rho,\rho u,\rho v,\rho E]^{T}$. The
vectors $\hat{F}$, and $\hat{G}$ correspond to the inviscid fluxes and the
vectors $\hat{F_{\nu}},\hat{G_{\nu}}$ constitute the viscous fluxes. Thus,

$F=\left[\begin{array}[]{c}\rho\hat{U}\\\ \rho u\hat{U}+\xi_{x}p\\\ \rho
v\hat{U}+\xi_{y}p\\\ (\rho E+p)\hat{U}\\\ \end{array}\right]\quad and\quad
F_{\nu}=\left[\begin{array}[]{c}0\\\ \xi_{x_{i}}\sigma_{i1}\\\
\xi_{x_{i}}\sigma_{i1}\\\ \xi_{x_{i}}(u_{j}\sigma_{ij}-\Theta_{i})\\\
\end{array}\right]$ (2)

The contravariant velocity component is given by $\hat{U}=\xi_{x}u+\xi_{y}v$
and the specific energy density is

$E=\frac{T}{\gamma(\gamma-1)M^{2}_{\infty}}+\frac{1}{2}(u^{2}+v^{2})$ (3)

The deviatoric stress ($\sigma_{ij}$) and heat flux vector ($\Theta_{i}$) are
given by:

$\sigma_{ij}=\mu\bigg{(}\frac{\partial\xi_{k}}{\partial x_{j}}\frac{\partial
u_{i}}{\partial\xi_{k}}+\frac{\partial\xi_{k}}{\partial x_{i}}\frac{\partial
u_{j}}{\partial\xi_{k}}-\frac{2}{3}\frac{\partial\xi_{l}}{\partial
x_{k}}\frac{\partial u_{k}}{\partial\xi_{l}}\delta_{ij}\bigg{)}$ (4)
$\Theta_{i}=-\frac{1}{(\gamma-1)M^{2}_{\infty}}\bigg{(}\frac{\mu}{Pr}\bigg{)}\frac{\partial\xi_{j}}{\partial
x_{i}}\frac{\partial T}{\partial\xi_{j}}$ (5)

The $\hat{G}$ and $\hat{G_{\nu}}$ vectors have the same structure, with $\xi$
variations replaced by $\eta$ variations. Stokes’ hypothesis is assumed for
the bulk viscosity coefficient i.e $\lambda=-2/3\mu$. Further details of the
derivations may be found in Vinokur [30]. Flow velocities and density are
normalized by the free stream velocity ($u_{\infty}$) and density
($\rho_{\infty}$), respectively, time is non-dimensionalized by $c/u_{\infty}$
(where $c$ is the chord length of the airfoil), while the pressure is non-
dimensionalized by $\rho_{j}U_{\infty}^{2}$; the perfect gas relationship then
reads $p=\rho T/\gamma M_{\infty}^{2}$, where $M_{\infty}$ refers to the
freestream Mach number. All plots shown in the paper are in non-dimensional
units. A constant Prandtl number of $0.72$ is assumed along with a constant
ratio for the specific heats $(\gamma=1.4)$. Sutherland’s law is used to model
the temperature dependence of viscosity.

The Roe scheme is used for the inviscid fluxes [31] with a third-order
reconstruction using the MUSCL approach [32]. Viscous terms are discretized
with a second-order central differencing technique. The temporal integration
is performed explicitly using a $4^{th}$ order Runge-Kutta scheme with a time
step size of $\Delta t=0.0001c/u_{\infty}$. The use of the above methodology
along with validation studies for flow past airfoils and cylinders can be
found in Ranjan et al. [33].

##  3 Overall Flowfield Features

2

(a) Computational Domain

(b) Instantaneous snapshot

Figure 1: Details of the computational domain and flow physics shown through
an instantaneous snapshot of the flowfield

We choose a transonic flow past a NACA0012 airfoil at a zero angle of attack
for the development of the ROMs. This flowfield has been widely explored using
both numerical [34, 35] and experimental [36] techniques to analyze the
unsteadiness in the wake including compressibility effects. At low Mach
numbers ($M_{\infty}<0.7$), an undulation develops in the wake close to the
trailing edge, which transitions to a pair of symmetrical recirculation
bubbles with an increase in the freestream Mach number. In the transonic
regime ($0.7<M_{\infty}<0.9$), a shock wave forms close to the trailing edge
because of the acceleration of the flow up to supersonic speeds over the
curved airfoil surface. The supersonic region is terminated by an unsteady
shock wave, which induces a separation of the boundary layer and leads to the
growth of a von Kármán instability. These instabilities result in an
alternating vortex shedding resembling a von Kármán vortex street downstream
of a cylinder.

Following Bourguet et al. [11], we choose a freestream Mach number of $0.85$
and a Reynolds number of $Re=5{,}000$ for the present simulation. The
transonic nature of the flow produces shock-induced unsteadiness in the
airfoil wake. A stretched O-grid comprising $150{,}000$ points is used to
discretize the domain, shown in fig. 1(a). The far-field boundaries are placed
at $150$ chord lengths in all directions and are modeled using characteristic
boundary conditions [37]. Figure 1b shows an instantaneous snapshot of the
flowfield. The Mach number contours highlight the acceleration of the flow
over the airfoil surface, giving rise to a region of supersonic flow around
the airfoil. A shock forms towards the trailing edge, as seen in the $\nabla
p$ contours, terminates this supersonic region. Vorticity contours display the
growth of the von Kármán instability resulting in a vortex street in the wake.
Shock oscillation and the coupled flapping shear layer induce strong pressure
fluctuations resulting in significant entropic gradients in the flow. These
gradients violate the assumptions for isentropic flow, which are only valid
for low-to-moderate Mach numbers and have been instrumental in the derivation
of several GP-based ROMs [5]. This flowfield is thus a challenging testbed for
evaluating the performance of the new NeuralGP algorithm, which has to date
only been demonstrated for an isentropic flow over an open cavity [16].

2

(a) Normal turbulent stresses

(b) Cross-streamwise velocity fluctuation ($v^{\prime}$) spectra at $x/c=2.0$

Figure 2: (a) Second order statistics of the flowfield shown using the normal
turbulent stresses; (b) Power Spectral Density (PSD) of velocity fluctuation
spectra on the airfoil axis

The strength, spatial extent and period of the vortex shedding downstream of
the airfoil are displayed using second order statistics. Figure 2(a). shows
normal turbulent stresses downstream of the airfoil. The streamwise
fluctuation energy ($\overline{u^{\prime}u^{\prime}}$) exhibits two symmetric
lobes and is more concentrated close to the trailing edge of the airfoil while
the cross-streamwise fluctuations ($\overline{v^{\prime}v^{\prime}}$) are more
extended in the streamwise direction. Significant anisotropy is evident in the
turbulent stresses with the cross-streamwise fluctuation energy
($\overline{v^{\prime}v^{\prime}}$) being an order of magnitude larger than
the streamwise fluctuations. The peak values of $0.0912$ and $0.167$ for the
streamwise and cross-streamwise normal stresses, respectively, are within
$0.01$ of the values obtained in the benchmark simulations of Bouhadji and
Braza [35] for the same flow conditions. The frequency of the vortex shedding
is assessed using the power spectrum of the cross-streamwise velocity
fluctuations. Figure 2b. shows the $v^{\prime}$ spectrum on the axis of the
airfoil one chord length downstream of the trailing edge (at $x/c=2.0$). The
fundamental vortex shedding frequency (annotated as $f=1.412u_{\infty}/c$) and
two super-harmonics dominate the spectrum. The fundamental shedding frequency
lies within one frequency resolution bin of the values reported by Bouhadji
and Braza [35] for the same operating conditions, lending confidence in the
accuracy of these simulations in capturing the dominant flow dynamics.

##  4 Reduced Order Model

The ROM construction approach relies on a set of reduced basis functions
optimized to represent the flowfield. Typically, POD modes are selected as the
optimal basis [11, 5, 17] owing to their ability to capture the maximum
variance of the flow fluctuations and their mutual orthogonality, which
simplifies the Galerkin projection of the governing equations. The process
used to compute the POD basis is outlined, followed by a description of the
ROM construction.

###  4.1 Proper Orthogonal Decomposition Procedure

The following transformation of the primitive variables is applied prior to
POD mode calculation:

$\\{u,v,w,p,\rho\\}\longrightarrow q=\\{1/\rho,u,v,w,p\\}$ (6)

This facilitates the Galerkin projection of the full Navier-Stokes equations
onto a reduced POD basis to construct the ROM [11]. The corresponding spatial
POD basis functions
($\Phi_{i}=\\{\Phi_{i}^{1/\rho},\Phi_{i}^{u},\Phi_{i}^{v},\Phi_{i}^{w},\Phi_{i}^{p}\\}^{T}$)
are then obtained using $N_{t}$ snapshots of the flowfield sampled at $\Delta
t=5\times 10^{-3}$. The snapshot POD method of Sirovich [38] is used to
compute the POD spatial modes, the details of which, including the spatial and
temporal inner products employed, are provided for completeness in Appendix A.

Figure 3: Effect of the number of snapshots used for the POD procedure on the
modal energy for the leading six POD modes

The convergence of the POD modes is first established because of its
substantial influence on the stability of the resulting ROMs. This is achieved
by examining the effect of the number of snapshots $N_{t}$ used for the POD
procedure on the modal energies for the six leading POD modes, as shown in fig
3. For smaller values of $N_{t}$, the modal energies of the two leading POD
modes are a strong function of the number of snapshots considered. This
follows from the fact that the leading modes correspond to the lowest
fundamental vortex shedding frequency (fig 2) and thus, fully converging their
statistics requires snapshots distributed over a sufficiently long time
horizon. With an increase in $N_{t}$, however, the modal energies of all six
leading POD modes become invariant with respect to a change in $N_{t}$. Thus,
to ensure an appropriately converged POD basis, we use $N_{t}=2{,}000$
snapshots to compute the POD modes.

Figure 4: Fractional as well as the cumulative energy
($\sum_{i}^{N_{POD}}\lambda_{i}/\sum_{i}^{N_{t}}\lambda_{i}$) content of the
individual POD modes given by eigenvalues of the cross correlation matrix

2

(a) Streamwise velocity

(b) Pressure

Figure 5: Leading spatial POD mode shapes for the 5 leading modes

Figure 4 shows the eigenvalues of the $20$ leading POD modes as a fraction of
the energy in each mode, along with their respective cumulative energy
content. A rapid drop in the energy content of the higher POD modes is readily
observable; as $99.3\%$ of the fluctuation energy is retained with a
truncation beyond the 6 leading POD modes. Additionally, the leading modes
occur in pairs, as is common in flows with dominant convective instabilities.
Figure 5(a) and (b) display the leading spatial POD modes for streamwise
velocity and pressure fluctuations; for brevity, only odd numbered modes are
shown among the repeated mode pairs. The mode shapes represent convecting
vortices in the wake, with the primary mode reflecting the spacing of the
vortices in the von Kármán vortex street formed in the wake. Their
antisymmetric pattern has been observed in previous studies on flows
exhibiting a dominant von Kármán instability [39]. The higher modes represent
increasingly finer spatial scales and correspond to superharmonics of the
primary vortex shedding frequency.

###  4.2 Galerkin projection and calibration

A Galerkin projection is performed of the full 3D compressible Navier-Stokes
equations on the truncated POD basis comprising the leading $N_{POD}=6$ modes.
The spatial inner product used for the projection of the equations is
identical to that used for the classical POD algorithm (detailed discussion in
the context of equation 25 in Appendix 4.1):

$\left<\frac{\partial\tilde{q}}{\partial
t},\Phi_{i}\right>_{\Omega}+\left<\frac{\partial\tilde{F}_{j}}{\partial\xi_{j}},\Phi_{i}\right>_{\Omega}=\left<\frac{\partial\tilde{F}^{vis}_{j}}{\partial\xi_{j}},\Phi_{i}\right>_{\Omega}$
(7)

where the Einstein summation convention has been used. Galerkin projection
yields the following set of ordinary differential equations for each time
coefficient of the leading POD modes:

$\displaystyle\begin{split}&\dot{a}_{i}=C_{i}+\sum_{j=1}^{N_{pod}}L_{ij}a_{j}+\sum_{j=1}^{N_{pod}}\sum_{k=1}^{N_{pod}}Q_{ijk}a_{j}a_{k}\\\
&a_{i}(0)=\left<\tilde{q}_{0}-\overline{\tilde{q}},\Phi_{i}\right>_{\Omega}\end{split}$
(8)

As discussed previously, the time integration of the above model suffers from
instabilities [40, 41, 42, 43] resulting in a significant prediction error,
alteration of oscillation periods, and, finally, a divergence of the temporal
integration. To obtain a stable prediction that minimizes prediction errors,
in the present study we implement a calibration method similar to that used by
Bourguet et al. [11]. The calibration procedure modifies the behavior of the
ROM equations using calibration matrices ($L_{ij}^{cal}$ and $C_{i}^{cal}$)
added to the original ROM coefficient matrices ($L_{ij}$ and $C_{i}$). This
results in the following equation:

$\dot{a}_{i}=\left(C_{i}+C_{i}^{cal}\right)+\sum_{j=1}^{N_{pod}}(L_{ij}+L_{ij}^{cal})a_{j}+\sum_{j=1}^{N_{pod}}\sum_{k=1}^{N_{pod}}Q_{ijk}a_{j}a_{k}=f_{i}(L^{cal},C^{cal},\mathbf{a})$
(9)

We obtain the calibration matrices by solving an optimization problem that
seeks to minimize the ROM prediction error subject to the constraint of
calibration cost. Thus, the following cost function is minimized:

$\mathcal{J}(C^{cal},L^{cal},\theta)=\theta\mathcal{E}(C^{cal},L^{cal})+(1-\theta)\mathcal{C}(C^{cal},L^{cal})$
(10)

where $\mathcal{E}(C^{cal},L^{cal})$ and $\mathcal{C}(C^{cal},L^{cal})$
represent normalized ROM error and calibration cost, respectively. The two
terms are blended using the coefficient $\theta\in\left(0,1\right)$.

The choice of error metric and calibration cost has a significant bearing on
the calibrated ROM equations and merit a discussion. The ROM prediction error
is defined as:

$E(C^{cal},L^{cal})=\sum_{i=1}^{N_{pod}}\int_{0}^{T}\left(a_{i}^{pod}(t)-a_{i}^{pod}(0)-\int_{0}^{t}f_{i}(C^{cal},L^{cal},a^{pod})dt^{\prime}\right)^{2}dt.$
(11)

Here, $a_{i}^{pod}$ indicates the time coefficient of the $i^{th}$ POD mode
observed directly from the flow field, referred to hereafter as the reference
POD time coefficient. Thus, the ROM prediction error is calculated as the
difference between the change in the $i^{th}$ reference POD time coefficients
($a^{POD}_{i}$) and that predicted by the RHS of the ROM equation when using
the reference POD time coefficients 9 over a time horizon $t=0$ to $t$.
Alternative definitions for the ROM prediction error, such as the difference
between the instantaneous time rate of changes of $a_{i}$s and that predicted
by the RHS of the governing equation, may also be used [15]:

$E^{alt}(C^{cal},L^{cal})=\sum_{i=1}^{N_{pod}}\int_{0}^{T}\left(\left(\frac{da_{i}^{pod}}{dt}-f(C^{cal},L^{cal},a^{pod})\right)^{2}dt\right)$
(12)

Crucially, these expressions represent a linearized definition of the ROM
error which only consider the differences in the rate of change of POD
coefficients ($\dot{a}$) versus that predicted by the ROM at each time step.
While this choice results in a linearization of the optimization problem, it
ignores the cumulative effects of the small differences in the predicted
$\dot{a}$ over the entire calculated trajectory. The consequences of this
choice are analyzed in greater detail in section 6.

The error term $E$ is normalized using the ROM prediction error with no
calibration to obtain $\mathcal{E}$:

$\mathcal{E}(C^{cal},L^{cal})=\frac{E(C^{cal},L^{cal})}{E(0,0)}$ (13)

The calibration cost is defined as a measure of the calibration matrices
relative to the system dynamic matrices. Here, a sum of an $L_{2}$ norm of the
$C_{i}$ vector and a Frobenius norm of the $L_{ij}$ matrix is considered:

$\mathcal{C}(C^{cal},L^{cal})=\frac{||C^{cal}||^{2}+||L^{cal}||^{2}}{||C||^{2}+||L||^{2}}$
(14)

The calibration procedure then solves the optimization problem to minimize
$\mathcal{J}$ for a given value of $\theta$ resulting in specific calibration
matrices $L^{cal}$ and $C^{cal}$ obtained as a function of the blending
parameter $\theta$.

Minimizing $\mathcal{J}$ is equivalent to solving the following linear system
of size $N_{POD}+1$:

$A^{cal}(K^{cal}_{i})^{T}=b^{i}$ (15)

where $K^{cal}=[C^{cal},L^{cal}]$ and the matrices $A^{cal}$ and $b^{i}$ are
given as:

$\displaystyle A_{ij}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\left(\int_{0}^{t}a_{i}^{\star}dt^{\prime}\right)\left(\int_{0}^{t}a_{j}^{\star}dt^{\prime}\right)dt+\theta^{\star}\delta_{ij}$
(16) $\displaystyle b_{i}^{k}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\left(a_{k}^{pod}(t)-a_{k}^{pod}(0)-\int_{0}^{t}f_{k}\left(0,0,a^{pod}\right)dt^{\prime}\right)\left(\int_{0}^{t}a_{i}^{\star}dt^{\prime}\right)dt,$
(17)

with $a^{\star}=[1,a_{1},a_{2},\cdots,a_{NPOD}]$ and

$\theta^{\star}=\frac{1-\theta}{\theta}\frac{E(0,0)}{||C||^{2}+||L||^{2}}$
(18)

###  4.3 Neural Galerkin ROM using Differentiable Programming

The GP formulation of equation 8 gives rise to $C$ (Constant), $L$ (Linear)
and $Q$ (Quadratic) terms, that are then modified for greater stability via an
optimization and calibration process from the data. We remark that the data-
driven nature of eq. 9, while being structured as an ODE, makes it a good
candidate to be cast instead as a DiffProg problem. Specifically, the $C$ and
$L$ coefficients have a significant impact on the accuracy and long-time
stability of the GP-ROM [44, 11]. Since this entails considerable manual
intervention as calibration, we seek to automate the discovery of these
coefficients with DiffProg, and directly use the trajectories of the reduced
state ($a_{i}(t)$) from the high resolution simulation as the training data.

In the DiffProg problem setup, the “unknown" operators in equation 8 are
approximated with trainable parameters $p$ in a Neural Network (NN), such that
$C_{i}\approx C_{i}^{p}$ and $L_{ij}\approx L_{ij}^{p}$, so that the equation
is now,

$\displaystyle\begin{split}&\dot{a}_{i}=C_{i}^{p}+\sum_{j=1}^{N_{pod}}L_{ij}^{p}a_{j}+\sum_{j=1}^{N_{pod}}\sum_{k=1}^{N_{pod}}Q_{ijk}a_{j}a_{k}\\\
&a_{i}(0)=\langle\tilde{q}_{0}-\overline{\tilde{q}},\Phi_{i}\rangle_{\Omega}\end{split}$
(19)

Since the trainable parameters $p$ have to be learned, they are typically
initialized with random numbers, as is standard in deep learning literature
[45]. In the forward pass of our setup, we solve Eqn. 19 using standard ODE
solvers to estimate $a_{i}(t)$ trajectories over a stipulated time period
($T$).

The first estimate of the predicted trajectories will necessarily be
inaccurate since the operators are essentially random, and the rest of the ODE
propagates it in time. This error in trajectory is captured in the loss
function, which is defined as the $L_{2}$ norm of the error between the
calculated trajectories and those obtained from the high-resolution
simulation:

$\displaystyle L(\mathbf{a}(t_{0}))$
$\displaystyle=\sum_{i=1}^{N_{POD}}\int_{0}^{T}\bigg{\\{}{a_{i}^{ref}}(t_{1})-\int_{0}^{t_{1}}f_{i}(\mathbf{a},t)dt\bigg{\\}}^{2}dt_{1}$
(20)
$\displaystyle=MSE(ODESolve(\mathbf{a}(t_{0}),t_{0},f,\theta)-\mathbf{a}_{ref})$
(21)

where $f_{i}(\mathbf{a},t)$ represents the RHS of Eqn. 19 and $a_{i}^{ref}$
refers to the reference POD time coefficients. Crucially, the solution of the
ODE problem with given parameters $C_{i}^{p}$ and $L_{ij}^{p}$ is an integral
part of the loss function and the prediction error is calculated for the
entire trajectory over a time period (from $t=0$ to $T$). In the backward
pass, the NN parameters $p$ are updated by backpropagating the error through
the ODE. The updated $p$ is slightly more accurate than the previous iteration
and used in the forward pass for the next iteration, where the loss is
computed again. This forward-backward pass cycle is repeated till the error
has reduced to a satisfactory value, which is user-defined and based on the
accuracy requirements of the application. We outline this process in Fig. 6,
as it has some distinct advantages compared to other methods in literature:

Figure 6: Schematic of NeuralGP training by treating the linear and constant
ROM-ODE coefficients as trainable parameters

  * •

Both $C_{i}$ and $L_{ij}$ are treated as trainable parameters in the NN and
they interact in a physical manner governed by the form of the projected
Navier-Stokes equations.

  * •

DiffProg allows for NN models to be back propagated seamlessly thru the ODE
solver via adjoint sensitivity analysis and a combination of forward/reverse
mode Automatic Differentiation (AD). The resulting abstract, powerful
framework allows for arbitrary neural network models to be tightly integrated
within the context of the well-known GP equations.

Thus, contrary to the previously adopted deep learning approaches for for
fluid flows [19, 46, 21, 47] that use black-box neural networks employing
standardized non-linearities, the NeuralGP approach leverages the power of the
DiffProg paradigm to derive coefficients to a known non-linearity (the
Galerkin ROM-ODE).

Our method may thus be viewed as an extension to the non-intrusive ROM
techniques proposed by Peherstorfer and Willcox [17], Ghattas and Willcox [23]
which seek to obtain the physically consistent operator matrices for a known
form of the ROM. These methods use training data composed of ROM trajectories
calculated directly from the full-order model sampled at discrete points and
solve a linear optimization problem to obtain the ROM operator matrices. In
contrast, in the current approach, the operators are learned by minimizing a
more rigorous prediction error defined over the entire trajectory obtained by
integrating the ROM. This rigorous definition of the error results in robust
learned data-driven models while retaining the key advantages of physically
derived ROM ODEs [23]. Two features are noteworthy. First, the computational
cost of training is very low since retention of physics as its Galerkin ODE
structure allows leveraging of fast, cheap ODE solvers instead of expensive,
complex deep NN architectures. These deep networks, like LSTMs and CNNs, are
black boxes that need to represent the entire physics, which requires
significant computing resources (several GPUs) to learn even modest-sized
datasets [19, 48, 49]. Second, the algorithm attempts to learn a physically
consistent system governing equation rather than fitting the data to a generic
model for prediction; this alleviates concerns regarding over-training the
neural network for a dataset at the cost of generalizability. As noted
earlier, the use of differentiable programming approach has demonstrated
substantial advantages for isentropic flows over Galerkin
projection/calibration methods [16]. Here we demonstrate that the method shows
similar improvement for more general flows where isentropic behavior is not
observed.

##  5 Evaluation of NeuralGP-based ROM

The performance of the proposed approach in obtaining stable ROM coefficients
is evaluated using short term and long-term predictions. The method uses
reference POD coefficients over a time horizon of $t=1.5c/u_{\infty}$,
corresponding to two full periods of the leading POD mode. For testing, the
calibrated ROM was also obtained using the same time horizon and a blending
parameter of $\theta=0.9$, corresponding to a $2.25\%$ calibration cost (eq.
14). The GP-ROM equations are shown in appendix B. The NeuralGP method
underwent $2{,}000$ training iterations to reduce the loss function by five
orders of magnitude, which required $\approx 228$ seconds on an Intel i7-8550U
octa-core CPU running at $4$ GHz. The inference cost for both the calibrated
GP-ROM as well as the NeuralGP ROM are similar, since they both require an
ODE-solve.

###  5.1 Short time prediction

Figure 7: Reference time coefficients $a_{i}(t)$ of the 6 leading POD modes
(original) along with predictions from a) calibrated ROM (Cal), b)
Uncalibrated ROM (noCal) and c) NeuralGP ROM, obtained using Galerkin
Projection. Uncalibrated ROM quickly diverges compared to the other two
approaches.

Figure 7 shows the time coefficients $a_{i}(t)$ for the six leading POD modes,
along with predictions from the uncalibrated GP-ROM, calibrated GP-ROM, and
the NeuralGP-ROM. Time coefficients of all leading modes exhibit a sinusoidal
variation with the repeated pairs of eigenvalues discussed in fig. 4 having
the same period but with a $\pi/2$ phase lag to maintain their orthogonality
in time. The frequencies corresponding to the time coefficients also increase
with the mode number, along with a progressive reduction in amplitude.

The uncalibrated ROM diverges relatively quickly from the reference POD time
coefficients, especially in the higher POD modes, with rapid growth in
amplitudes and a change of the oscillation periods. Even for the two leading
POD modes, reasonable results are only obtained over the first oscillation
period, after which we can observe a significant deviation from the reference
POD coefficients. Even with the relatively small additions to the $L_{ij}$ and
$C_{i}$ matrices (as indicated by the small calibration cost), the calibrated
ROM tracks the time coefficients obtained from the simulation accurately over
the short time horizon. Finally, the results of the NeuralGP ROM closely track
the original POD time coefficients, resulting in virtually indistinguishable
differences between the calibrated GP and NeuralGP ROM results from the
original POD time coefficients in the short time horizon. While the close
correspondence of the ROM predictions with the reference POD time coefficients
in the short time horizon is encouraging, it only serves as a sanity check for
the ROM methodologies. This is because the short time horizon over which the
predictions are being tested here is identical to that used for calibrating
and training the GP and NeuralGP ROMs, respectively. We now undertake a more
rigorous test of the ROMs by evaluating their performance over a longer term
prediction horizon.

###  5.2 Long time predictions

(a) Long term predictions

(b) Zoomed in predictions in the terminal prediction window

(c) Spectral distribution at the beginning and end of the prediction window
with the fundamental shedding frequency($St=1.42$) and its harmonics being
identified in each pair of POD modes

Figure 8: Long time prediction of the 4 leading POD modes using the calibrated
GP-ROM and NeuralGP ROMs along with a detailed view of the terminal prediction
window

The long term predictive capabilities of both ROMs are evaluated by
integrating the respective ROM-ODEs over a longer time horizon
($t\in(0,30c/u_{\infty})$). Figure 8(a) shows the energy-normalized long time
predictions of the three leading odd number modes obtained using both the
calibrated GP and NeuralGP ROMs compared to the original POD time
coefficients. The even modes exhibit similar behaviors as their corresponding
odd mode pairs, and are therefore not displayed here for brevity. The
predictions from both ROMs follow the original POD time coefficients closely
until $t\approx 10c/u_{\infty}$. Beyond this time, the NeuralGP predictions
track the original POD time coefficients closely, with no indication of
divergence, though the GP-ROM shows significant errors in prediction. The
terminal prediction window at the end of the long time prediction horizon is
shown in figure 8b. There is relatively little degradation in the performance
of the NeuralGP while the calibrated GP-ROM chiefly changes in the oscillation
periods, along with a slow growth in the amplitude. These gradual changes in
the frequencies of each time coefficient lead to substantial errors in the
predicted phase of the modes, as observed for modes $3$ and $5$ in the
terminal prediction window.

The performance of the NeuralGP ROM is further investigated using Power
Spectral Densities (PSD) of the time coefficients. To understand the
differences in the spectral content brought about by the longer term
integration of the ROMs, the analysis is performed on smaller time windows of
the long-time prediction horizon ($tU_{\infty}/c\in(0,150)$). Figure 8(c)
shows the differences in the spectral content of the three leading mode pairs
within time windows of $10c/U_{\infty}$ at the start
($tU_{\infty}/c\in(0,10)$) and end ($tU_{\infty}/c\in(135,145)$) of the long
time prediction horizon. The leading POD modes captures the fundamental
shedding frequency marked in fig. 2(b) while the superharmonics are captured
by the progressively higher POD mode pairs. The extents of the $x$-axes of the
PSD plots are adjusted to represent the correspondingly higher temporal
frequencies of the higher POD time coefficients. The GP-ROM prediction shows a
substantial difference in the spectral content over the course of the entire
prediction horizon. The calibrated GP-ROM yields much better results, but
there is nonetheless a shift in the frequency content and amplitude growth in
the time window close to the end of the long time prediction horizon
$t\in(135,145)$. The NeuralGP predictions maintain a near identical spectral
distribution matching closely with the spectral signature of the vortex
shedding (fig. 2(b)) irrespective of the time window selected. This appears to
be a systematic characteristic of NeuralGP, since it minimizes the full
trajectory error in a nonlinear fashion. We study this in detail in the next
section.

##  6 System dynamics analysis of the ROM-ODE

The attractive properties of the long term predictions of the NeuralGP-ROM are
now traced to aspects of its underlying mathematical formulation. As mentioned
in the previous section, the calibration procedure alters the ROM coefficients
obtained from the Galerkin projection of the Navier-Stokes equations on the
POD basis to stabilize the ROM predictions. This is achieved by solving an
optimization problem seeking to minimize a cost function ($\mathcal{J}$
defined in eq. 10), which is composed of a linearized definition of the
prediction error ($\mathcal{E}$ defined in eq. 13) and a calibration cost
penalty ($\mathcal{C}$ defined in eq. 14), balanced via a blending parameter
($\theta$). These parameters dictate the properties of the final calibrated
ROM dynamic matrices, therefore studying their variations can provide
important insights into the calibration process itself.

2

(a) Full Range of $\theta$($\theta\in(0,1)$)

(b) $\theta$ close to 1($\theta\in(0.9,1)$)

Figure 9: Variation of calibration cost ($\mathcal{C}$) and prediction error
($\mathcal{E}$) with blending parameter ($\theta$)

Figure 9 shows the variation of the calibration cost and linearized error as a
function of the blending coefficient ($\theta$). The plot also shows the full
trajectory error of the ROM calculated as the $L_{2}$ norm of the error
between the reference POD time coefficients and those obtained from the ROM
prediction over the entire calibration horizon similar to eq. 20. Increasing
values of $\theta$ represent a progressively diminishing relative importance
of the calibration cost in the optimization problem and thus, result in a
reduction of the linearized ROM prediction error and an increase in the
calibration cost, as expected. However, for all values of $\theta$, the full
trajectory error is several orders of magnitudes larger than the linearized
prediction error. This is because the latter, by definition (eq. 11), does not
account for the cumulative nature of errors in predicted trajectories when
integrating the ROM-ODE with inaccuracies in the coefficients. While both the
linearized and full trajectory error exhibit a similar variation (decreasing
rapidly with increasing values of $\theta$), this trend breaks down for
$\theta$ values very close to $\theta=1$ where the full trajectory error
increases even with a decrease in the linearized error metric, as shown in
fig. 9(b). This is likely the result of the near singular nature of the
problem for values of $\theta$ close to one (Appendix B). These variations of
the error metrics and calibration costs during the calibration procedure are
now compared to those obtained during the training process of the NeuralGP
algorithm.

To ensure a illustrative evaluation, we initialized the values for the
trainable ROM coefficients in the NeuralGP ROM to those obtained from the
Galerkin Projection procedure. While this is unnecessary for the convergence
of the DiffProg algorithm as shown by Mohan et al. [16] who used random
numbers to initialize the ROM-ODE coefficients, it affords key insights into
the learning process by enabling a meaningful tracking of parameters similar
to those used for the calibration method. To study this, the linearized and
full trajectory error (NeuralGP loss) as well as an equivalent calibration
cost are tracked throughout the learning process and compared to those
observed in the calibration procedure. Since the NeuralGP process does not
consider a calibration cost for the error minimization, we define an
equivalent calibration cost as:

$\mathcal{C}_{learn}=\frac{||C^{learn}-C^{GP}||^{2}+||L^{learn}-L^{GP}||^{2}}{||C^{GP}||^{2}+||L^{GP}||^{2}}$
(22)

where $C^{learn}$, $L^{learn}$ and $C^{GP}$, $L^{GP}$ refer to the learned and
Galerkin Projection based ROM-ODE coefficients respectively.

Figure 10 shows the variations in these performance metrics throughout the
training process of the Neural GP procedure.

Figure 10: Variation in prediction error ($\mathcal{E}$) and equivalent
calibration($\mathcal{C}$) costs with training iterations

The full trajectory error exhibits a monotonic decrease over the course of the
training, highlighting the ability of the DiffProg algorithm to converge
rapidly to the trainable parameters of the actual system dynamics governing
the ROM-ODE. At the end of the training, the full trajectory error obtained
using the NeuralGP algorithm is lower than that obtained in the calibrated-GP-
ROM with even the highest values of $\theta$. This illustrates the superior
performance of the NeuralGP algorithm in long-term predictions compared to
those obtained from the calibrated GP-ROM.

Some additional important differences in the variation of the error metrics
compared to their behavior for the calibration process (shown in fig. 9) can
also be observed. Except for the initial few learning iterations, at all
points during the learning process, the full trajectory error is lower than
the linearized definition. Furthermore, the calibration cost grows initially
to a value slightly greater than $0.025$ before marginally diminshing and
remaining constant throughout the learning process. Even with a constant value
of the calibration cost, however, both metrics of the prediction error
decrease by several orders of magnitudes. This is an important difference
compared to the calibration process, where the calibration cost increases
monotonically with a decrease in the prediction error. This indicates that the
Frobenius norm chosen for the calibration cost (eq. 14) is not fully
representative of the actual impact of the calibration matrix on the system
dynamics and may form a poor constraint on the calibration problem.
Interestingly, even without considering calibration cost i.e with a $\theta$
value close to one, the NeuralGP procedure remains more efficient in reducing
the achieved total trajectory error.

These trends highlight the two key advantages of NeuralGP technique: the
choice of the error metric being minimized and the consideration of a
calibration cost. Specifically, unlike the calibration procedure, the NeuralGP
procedure seeks to minimize the full trajectory error in prediction with no
constraint on the calibration cost. Besides the above general trends, an
interesting observation regarding the NeuralGP training procedure over the
first few training iterations is that the full trajectory error decreases
sharply with a rapid increase in the linearized error metric and the
calibration cost. This is distinct from the general trend of the full
trajectory and linearized error metrics having a similar qualitative variation
(decreasing or increasing together). This rapid change in the error metrics
over the first few iterations is followed by a more gradual change where both
error metrics exhibit a qualitatively similar decay while the calibration cost
increases over the subsequent learning iterations. The rate of decay in both
error metrics increases again at approximately thousand training iterations
and continues to decay more quickly for the rest of the training epoch. We
explore this behavior over the first few training iterations in a later
section.

The differences in the trends exhibited by the ROM prediction metrics for the
two methods are now explored further by considering their respective system
dynamics. To derive a theoretical basis for comparison between the two ROM
methodologies, their ROM-ODE coefficients are analyzed. First, the relative
importance of the terms on the RHS of eq. 8 were scrutinized over the short
term prediction horizon by individually calculating their contribution to the
temporal rate of change $\dot{a}$.

Figure 11: Relative magnitudes of the terms on the RHS of the ROM-ODE (eq. 8)

Figure 11 compares the magnitudes of the change in the ROM coefficients at
each time step ($da$) to the contributions from the constant ($dConst$),
linear ($dLin$) and quadratic ($dQuad$) terms for the first six POD modes. The
linear term has the largest impact on the system dynamics for all the leading
modes tested, contributing to most of the rate of change at each time step. In
comparison, the quadratic term has a contribution one order of magnitude
smaller than that of the linear term while the impact of the constant term is
smaller still. In view of the dominant influence of the linear term on the
ROM-ODE (eq. 8), we analyze the eigenvalues of the linear matrix ($L_{ij}$) to
explain the impact of calibration and the NeuralGP method on the system
dynamics.

2

(a) NeuralGP

(b) Calibration

Figure 12: Trajectories of the eigenvalues for both the calibration and
NeuralGP procedures.

The trajectories of the eigenvalues with a variation in $\theta$ and training
iterations for the calibrated GP-ROM and the NeuralGP method respectively are
shown in fig. 12. Throughout their trajectories, the eigenvalues occur as
conjugate pairs with the individual eigenvectors corresponding to the POD mode
pairs as labeled. The perfectly conjugate eigenvalues reflect the fact that
the leading POD mode pairs have the same temporal frequencies separated by a
$90^{o}$ phase difference. The $x$-axis $Re(Eig)$ in Fig 12 denotes the real
part of the eigenvalues, and $Re(Eig)>0$ in the positive half-plane indicates
that the associated modes are unstable, i.e., growing with time. Similarly,
$Re(Eig)<0$ in the negative half-plane indicates modes that are decaying with
time, while $Re(Eig)=0$ are neutral modes that are neither growing nor
decaying. The ROM comprises a linear combination of these modes, and their
individual stability characteristics dictate the overall ROM stability.
Therefore, the ROM would become inaccurate if all the modes are growing or
decaying. A ROM typically comprises modes that have both $Re(Eig)<0$ and
$Re(Eig)>0$, and the challenge is to construct the modes with stability values
that ensure the ROM does not diverge. The dominant modes also have a higher
impact on ROM stability than the lower modes, as can be inferred from Fig. 4.
As a result, both the calibration in GP-ROM and the differential programming
approach in NeuralGP aim to perform a balancing act between mode dominance and
its stability, to ensure ROM accuracy for as long as possible.

The physical mechanism associated with the instability of the uncalibrated GP-
ROM equations can be inferred from the corresponding eigenvalues of the
uncalibrated $L_{ij}$ matrix. We show this analysis in fig. 12 where cross-
markers denote the initial eigenvalues of $L_{ij}$ before learning. These
initial eigenvalues have positive real parts with progressively greater values
for higher POD modes. The greatest of these, corresponding to the mode $(5,6)$
pair, exhibit real parts as high as $\approx 2.0$. This is consistent with the
rapid destabilization of the $(5,6)$ mode pair within the first few
oscillation cycles in the uncalibrated GP-ROM predictions as shown in fig. 7.

Both calibrated GP and NeuralGP ROM stabilize the eigenvalues by “damping"
their real parts and making small corrections to their imaginary parts. We
plot these continuous changes to the eigenvalues with training as a trajectory
in Fig. 12(a), since it shows the intermediate eigenvalues obtained in the
path to the final optimized $L_{ij}$. However, the trajectories followed for
each ROM exhibit interesting differences. The NeuralGP method pushes the
eigenvalues corresponding to the leading modes towards the negative half of
the complex plane within a few initial training iterations (every tenth
training iteration shown with markers). Such rapid displacement of the
eigenvalues quickly decreases the full trajectory error, and increases the
equivalent calibration cost within the first few iterations of the training
process, as seen in Fig. 10. Following this swift change in the leading
eigenvalues during the initial training period, a reversal in their direction
of motion is evident. A similar reversal is also observed for the (3,4) mode
pair later during the training. This is in contrast to the typical GP-ROM
calibration process (Fig. 12(b)) which linearly moves the eigenvalues with
only relatively small adjustments for the leading mode pair.

The explanation for the distinct trajectories of the eigenvalues lies in the
differences in the optimization problem being solved in both cases. As
discussed briefly in the previous section, the calibration process studied
here (based on Tikhonov regularization) solves a linearized form of a non-
linear constrained optimization problem [15]. Thus, a linearized error
function (eq. 11) is minimized during the calibration process and the movement
of the eigenvalues also reflects the linearity of the optimization problem
being solved with linear unidirectional trajectories. In contrast, the
NeuralGP method solves a non-linear minimization problem addressing the
prediction error for the full trajectory across a prediction horizon. Since
this algorithm integrates the ROM-ODE at each training iteration over a short
time window, it captures the interactions between the modes. Thus, initially,
the NeuralGP minimizes the full trajectory error by excessively damping the
leading modes, while reducing the amplitudes of the highly unstable $(5,6)$
mode pair as well. This is a consequence of the coupled nature of the ROM-ODE
wherein damping the leading modes also has the effect of damping the higher
POD modes, due of their higher contribution to the RHS of the ROM equations.
For example, the ROM equation for the time coefficient corresponding to the
$a_{5}$ POD mode is given by:

$\frac{da_{5}}{dt}=C_{5}+L_{5,1}a_{1}+L_{5,2}a_{2}+L_{5,3}a_{3}+L_{5,4}a_{4}+L_{5,5}a_{5}+L_{5,6}a_{6}+\sum_{j=1}^{N_{POD}}\sum_{k=1}^{N_{POD}}Q_{5,j,k}a_{j}a_{k}$
(23)

Here, the overpowering contribution of the leading $a_{1}$ and $a_{2}$ modes
on the RHS is visible; reducing their magnitude has the effect of lowering the
contribution of the higher $a_{5}$ and $a_{6}$ POD modes as well. At the
beginning of the training interval, the fastest method to reduce the
prediction error is to damp the leading POD modes since correcting highly
unstable eigenvalues of the $(5,6)$ mode pair requires a much more significant
alteration of the dynamic matrix. However, this trend is quickly reversed
within a few training iterations when the inaccuracies generated by the
additional dissipation added to the $(1,2)$ modes become a dominant factor in
the total prediction error.

The position of the final stabilized eigenvalues in the complex plane also
exhibit distinct and important features of the NeuralGP method. Initially, the
most unstable eigenvalues corresponding to the $(5,6)$ mode pair are pushed
far into the negative half-plane. Small corrections apply to the eigenvalues
representing the leading POD modes ($1,2$ mode pair) whose stabilized
eigenvalues retain slightly positive real parts. This shows that the two
leading POD modes are, in fact, growing modes which are non-linearly saturated
to yield a stationary solution. Despite the overall similarities between the
stabilized eigenvalues obtained using the two methods, a crucial difference
emerges from fig. 12(b). The NeuralGP procedure stabilizes the initially most
unstable $(5,6)$ mode pair, shifting these values further into the negative
half-plane. Even with a high value of $(\theta=0.95)$ for the blending
parameter, the calibration procedure fails to stabilize this mode pair,
resulting in discrepancies in the long-term predictions observed in fig. 8.
The large errors in long-term predictions from relatively insignificant
changes in the eigenvalues serves to highlight the sensitivity of the ROM-ODE
integration to its coefficients because of an accumulation of small errors at
each time step.

##  7 Summary and Discussion

A method using differentiable programming is developed to construct a Machine
Learning-based Reduced Order Model (ROM) that alleviates concerns about the
black-box nature of the typical deep learning-based predictive tools while
improving interpretability and prediction confidence. The performance of the
ROM is demonstrated and evaluated on a compressible, transonic flow over a
NACA0012 airfoil, which contains shocks and vortices that can challenge the
accuracy of long-time forecasting. The results are compared with a traditional
approach based on a Galerkin projection (GP) of the full Navier stokes
equation on a truncated POD spatial basis, which is susceptible to
instabilities as reported in the literature.

The GP-ROM calibration procedure required for accurate predictions using the
uses optimization of linearized ROM prediction errors and incurs a calibration
cost of a Tikhonov regularization framework. The differentiable programming
method as used in the current method maintains the form of the ODE usually
obtained in the GP-ROMs but learns the coefficients directly from training
data using neural networks and gradient descent. Therefore, it effectively
blends the strengths of both the ML and GP-based techniques for predictive
modeling. The present approach displays several advantages over traditional
black-box neural networks typically used in machine learning-based flow
prediction tools (LSTMs or CNNs). These include lower requirements for
computational resources and training data while mitigating the risk of
overfitting for a given data set. Although the short-time predictions are
similar to GP-ROM methods, the new approach avoids the known divergence issues
of calibrated GP-ROMs even for the longest prediction horizons tested..

A dynamical systems perspective of the ROM-ODEs provides insights into the
distinct features of the proposed approach. In particular, an examination of
the eigenvalues identify the linear term in the ROM-ODE to be a dominant
contributor to the system dynamics. The linear dynamics matrix from the
uncalibrated GP-ROM-ODE contains eigenvalues on the right complex half-plane,
indicating an exponential growth in the predictions, consistent with the
observed instabilities of the ROM trajectories. Both the calibrated ROM and
the current NeuralGP procedures stabilized the unstable eigenvalues by
displacing them towards the left half of the complex plane with the final
eigenvalues occupying similar positions, except for the higher POD modes,
which are not adequately damped by the calibration procedure. Despite the
general similarity between the stabilized eigenvalues, substantial differences
become evident in the trajectories followed by the eigenvalues over the
training/calibration process. These differences are shown to result from the
NeuralGP algorithm effectively minimizing a more rigorous error norm. The data
driven nature, as well as lower error, bound achieved by the NeuralGP method,
shows considerable promise in blending deep learning with physically
interpretable methods and should facilitate the construction of ROMs for more
complex flows, including fully turbulent flowfields.

## Acknowledgements

A.T.M and D.L. are supported by the LDRD (Laboratory Directed Research and
Development) program at Los Alamos National Laboratory (LANL) under Contract
No. 89233218CNA000001 with the U.S. Department of Energy/National Nuclear
Security Administration under projects 20190058DR and 20220104DR.

##  Appendix A Details of the POD procedure

The spatial POD basis functions
($\Phi_{i}=\\{\Phi_{i}^{1/\rho},\Phi_{i}^{u},\Phi_{i}^{v},\Phi_{i}^{w},\Phi_{i}^{p}\\}^{T}$)
are sought that maximize the projection of the fluctuation fields
($q(\mathbf{x,t})$) and optimally represent the second order statistics of the
flow:

$\max_{\mathbf{\Phi(\mathbf{x})}}\frac{\overline{\langle
q(\mathbf{x},t),\Phi(\mathbf{x})\rangle}_{\Omega}}{\langle\Phi(\mathbf{x}),\Phi(\mathbf{x})\rangle_{\Omega}}$
(24)

Here, the $\overline{.}$ operator refers to a time averaging operation and the
$\langle.\rangle_{\Omega}$ refers to the spatial inner product defined on the
spatial domain $\Omega$, which, for two separate realizations of the state
vector ($q^{I}$ and $q^{II}$), is given by:

$\langle
q^{I},q^{II}\rangle_{\Omega}=\sum_{ivar=1}^{nvar}\frac{1}{\sigma_{ivar}^{2}}\int_{\Omega}q^{I}_{ivar}q^{II}_{ivar}d\mathbf{x}$
(25)

Thus, the contribution of the $i^{th}$ variable of the state vector
($q^{I}_{ivar}$) to the spatial inner product is normalized by
$\sigma_{ivar}^{2}$, the spatio-temporally averaged energy of the variable:

$\sigma_{ivar}^{2}=\int_{\Omega}\overline{q^{2}_{ivar}}d\mathbf{x}$ (26)

Using the above definition of the spatial inner product, the time cross
correlation matrix is defined as:

$K(t_{1},t_{2})=\langle
q(\mathbf{x},t_{1}),q(\mathbf{x},t_{2})\rangle_{\Omega}$ (27)

Following the snapshot POD method [38], the maximization of the flowfield
projections (equation 24) is equivalent to solving the following eigenvalue
problem to obtain the time variation of each spatial mode:

$\overline{K(t_{1},t)a_{i}(t)}=\lambda_{i}a_{i}(t_{1})$ (28)

The above time coefficients ($a_{i}(t)$) are used to obtain the spatial basis
functions by calculating the projection of the simulation snapshots on these
time coefficients:

$\Phi_{i}(\mathbf{x})=\frac{1}{\sqrt{\lambda_{i}}}\overline{\tilde{q}(\mathbf{x},t)a_{i}(t)}$
(29)

Finally, the simulation snapshots ($q$) can be expressed using the time
coefficients ($a_{i}(t)$) and the spatial basis functions
($\Phi_{i}(\mathbf{x})$) after a decomposition into mean ($\overline{q}$) and
fluctuating component ($\tilde{q}$). The optimality of the POD spatial basis
then allows for an error bound truncation using only the leading $N_{POD}$
modes:

$q=\overline{q}+\tilde{q}=\overline{q}+\sum_{i=0}^{\infty}a_{i}(t)\Phi_{i}(\mathbf{x})\approx\overline{q}+\sum_{i=0}^{N_{POD}}a_{i}(t)\Phi_{i}(\mathbf{x})$
(30)

##  Appendix B Details of the GP-ROM coefficients

The coefficients of the quadratic ROM-ODE (eq. 8) are given by:

$\displaystyle C_{i}$ $\displaystyle=\langle
F^{\alpha}_{11}-A^{\alpha}_{11},\Phi_{i}\rangle_{\Omega}$ (31) $\displaystyle
L_{ij}$ $\displaystyle=\langle
F^{\alpha}_{1(j+1)}+F^{\alpha}_{(j+1)1}-A^{\alpha}_{1(j+1)}-A^{\alpha}_{(j+1)1},\Phi_{i}\rangle_{\Omega}$
$\displaystyle Q_{ijk}$ $\displaystyle=\langle
F^{\alpha}_{(j+1)(k+1)}-A^{\alpha}_{(j+1)(k+1)},\Phi_{i}\rangle_{\Omega}$

where the Greek indices ($\alpha,\beta$) indicate a summation over the
coordinate directions. The $F_{ij}$ and $A_{ij}$ matrices are given by:

$A^{i}_{jk}=\begin{bmatrix}\Phi_{j}^{\star
u_{i}}\Phi_{k,i}^{\star(1/\rho)}-\Phi_{j}^{\star(1/\rho)}\Phi_{k,i}^{\star
u_{i}}\\\ \Phi_{j}^{\star u_{i}}\Phi_{k,i}^{\star
u_{1}}+\Phi_{j}^{\star(1/\rho)}\Phi_{k,i}^{\star p}\delta_{1i}\\\
\Phi_{j}^{\star u_{i}}\Phi_{k,i}^{\star
u_{2}}+\Phi_{j}^{\star(1/\rho)}\Phi_{k,i}^{\star p}\delta_{2i}\\\
\Phi_{j}^{\star u_{i}}\Phi_{k,i}^{\star
u_{3}}+\Phi_{j}^{\star(1/\rho)}\Phi_{k,i}^{\star p}\delta_{3i}\\\
\gamma\Phi_{j}^{\star p}\Phi_{k,i}^{\star u_{i}}+\Phi_{j}^{\star
u_{i}}\Phi_{k,i}^{\star p}\end{bmatrix}\text{ and
}F^{i}_{jk}=\begin{bmatrix}0\\\
\Phi_{j}^{\star(1/\rho)}\tau^{\star}_{1ik,i}\\\
\Phi_{j}^{\star(1/\rho)}\tau^{\star}_{2ik,i}\\\
\Phi_{j}^{\star(1/\rho)}\tau^{\star}_{3ik,i}\\\
\frac{\gamma\mu}{P_{r}}\left(\Phi_{j}^{\star
p}\Phi_{k}^{\star(1/\rho)}\right)_{,ii}+(\gamma-1)\Phi_{j,i}^{\star
u_{\alpha}}\tau^{\star}_{\alpha ik}\end{bmatrix}.$ (32)

Here, $\tau^{\star}_{ijk}=\mu(\Phi_{k,j}^{\star u_{i}}+\Phi_{k,i}^{\star
u_{j}}-2/3\Phi_{k,\alpha}^{\star u_{\alpha}}\delta_{ij})$ and $\Phi^{\star}$
refers to the augmented POD basis comprising of the mean flow and the spatial
POD modes: $\Phi^{\star}=[\overline{q},\Phi_{1},\cdots,\Phi_{N_{pod}}]$. The
matrices $C_{i}$, $L_{ij}$, and $Q_{ijk}$ of equation 8 can be pre-computed
using only the spatial POD modes. The integration of the ROM equation itself
requires $N_{POD}(1+N_{POD}+N_{POD}^{2})$ scalar products per time step and is
much cheaper than integrating the governing equations (equation 1).

## References

  * Gad-el Hak [2001] Gad-el Hak, M., “Flow control: The future,” _Journal of aircraft_ , Vol. 38, No. 3, 2001, pp. 402–418.
  * Bons et al. [2005] Bons, J. P., Hansen, L. C., Clark, J. P., Koch, P. J., and Sondergaard, R., “Designing low-pressure turbine blades with integrated flow control,” _Turbo Expo: Power for Land, Sea, and Air_ , Vol. 47268, 2005, pp. 1079–1091.
  * Samimy et al. [2007] Samimy, M., Kim, J. H., Kastner, J., Adamovich, I., and Utkin, Y., “Active control of high-speed and high-Reynolds-number jets using plasma actuators,” _Journal of Fluid Mechanics_ , Vol. 578, 2007, p. 305.
  * Wei and Freund [2006] Wei, M., and Freund, J. B., “A noise-controlled free shear flow,” _Journal of Fluid Mechanics_ , Vol. 546, 2006, p. 123.
  * Rowley et al. [2004] Rowley, C. W., Colonius, T., and Murray, R. M., “Model reduction for compressible flows using POD and Galerkin projection,” _Physica D: Nonlinear Phenomena_ , Vol. 189, No. 1-2, 2004, pp. 115–129.
  * Noack et al. [2011] Noack, B. R., Morzynski, M., and Tadmor, G., _Reduced-order modelling for flow control_ , Vol. 528, Springer Science & Business Media, 2011.
  * Rempfer [2000] Rempfer, D., “On low-dimensional Galerkin models for fluid flow,” _Theoretical and Computational Fluid Dynamics_ , Vol. 14, No. 2, 2000, pp. 75–88.
  * Lorenzi et al. [2016] Lorenzi, S., Cammi, A., Luzzi, L., and Rozza, G., “POD-Galerkin method for finite volume approximation of Navier–Stokes and RANS equations,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 311, 2016, pp. 151–179.
  * Benner et al. [2015] Benner, P., Gugercin, S., and Willcox, K., “A survey of projection-based model reduction methods for parametric dynamical systems,” _SIAM review_ , Vol. 57, No. 4, 2015, pp. 483–531.
  * Chinesta et al. [2016] Chinesta, F., Huerta, A., Rozza, G., and Willcox, K., “Model order reduction,” _Encyclopedia of Computational Mechanics, 2Eds., Wiley_ , 2016\.
  * Bourguet et al. [2011] Bourguet, R., Braza, M., and Dervieux, A., “Reduced-order modeling of transonic flows around an airfoil submitted to small deformations,” _Journal of Computational Physics_ , Vol. 230, No. 1, 2011, pp. 159–184.
  * Callaham et al. [2021] Callaham, J. L., Brunton, S. L., and Loiseau, J.-C., “On the role of nonlinear correlations in reduced-order modeling,” _arXiv preprint arXiv:2106.02409_ , 2021.
  * Bergmann [2004] Bergmann, M., “Optimisation aérodynamique par réduction de modèle POD et contrôle optimal: application au sillage laminaire d’un cylindre circulaire,” Ph.D. thesis, Institut National Polytechnique de Lorraine, 2004.
  * Couplet et al. [2005] Couplet, M., Basdevant, C., and Sagaut, P., “Calibrated reduced-order POD-Galerkin system for fluid flow modelling,” _Journal of Computational Physics_ , Vol. 207, No. 1, 2005, pp. 192–220.
  * Cordier et al. [2010] Cordier, L., El Majd, B. A., and Favier, J., “Calibration of POD reduced-order models using Tikhonov regularization,” _International journal for numerical methods in fluids_ , Vol. 63, No. 2, 2010, pp. 269–296.
  * Mohan et al. [2021] Mohan, A. T., Nagarajan, K., and Livescu, D., “Learning Stable Galerkin Models of Turbulence with Differentiable Programming,” _arXiv preprint arXiv:2107.07559_ , 2021.
  * Peherstorfer and Willcox [2016] Peherstorfer, B., and Willcox, K., “Data-driven operator inference for nonintrusive projection-based model reduction,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 306, 2016, pp. 196–215.
  * Hochreiter and Schmidhuber [1997] Hochreiter, S., and Schmidhuber, J., “Long short-term memory,” _Neural computation_ , Vol. 9, No. 8, 1997, pp. 1735–1780.
  * Mohan and Gaitonde [2018] Mohan, A. T., and Gaitonde, D. V., “A deep learning based approach to reduced order modeling for turbulent flow control using LSTM neural networks,” _arXiv preprint arXiv:1804.09269_ , 2018.
  * Ahmed et al. [2020] Ahmed, S. E., San, O., Rasheed, A., and Iliescu, T., “A long short-term memory embedding for hybrid uplifted reduced order models,” _Physica D: Nonlinear Phenomena_ , Vol. 409, 2020, p. 132471.
  * Deng et al. [2019] Deng, Z., Chen, Y., Liu, Y., and Kim, K. C., “Time-resolved turbulent velocity field reconstruction using a long short-term memory (LSTM)-based artificial intelligence framework,” _Physics of Fluids_ , Vol. 31, No. 7, 2019, p. 075108.
  * Maulik et al. [2020] Maulik, R., Lusch, B., and Balaprakash, P., “Non-autoregressive time-series methods for stable parametric reduced-order models,” _Physics of Fluids_ , Vol. 32, No. 8, 2020, p. 087115.
  * Ghattas and Willcox [2021] Ghattas, O., and Willcox, K., “Learning physics-based models from data: perspectives from inverse problems and model reduction,” _Acta Numerica_ , Vol. 30, 2021, pp. 445–554.
  * Griewank et al. [1989] Griewank, A., et al., “On automatic differentiation,” _Mathematical Programming: recent developments and applications_ , Vol. 6, No. 6, 1989, pp. 83–107.
  * Bartholomew-Biggs et al. [2000] Bartholomew-Biggs, M., Brown, S., Christianson, B., and Dixon, L., “Automatic differentiation of algorithms,” _Journal of Computational and Applied Mathematics_ , Vol. 124, No. 1-2, 2000, pp. 171–190.
  * Baydin et al. [2018] Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M., “Automatic differentiation in machine learning: a survey,” _Journal of machine learning research_ , Vol. 18, 2018.
  * Caruana et al. [2005] Caruana, D., Mignosi, A., Corrège, M., Le Pourhiet, A., and Rodde, A., “Buffet and buffeting control in transonic flow,” _Aerospace science and technology_ , Vol. 9, No. 7, 2005, pp. 605–616.
  * Ren et al. [2020] Ren, K., Chen, Y., Gao, C., and Zhang, W., “Adaptive control of transonic buffet flows over an airfoil,” _Physics of Fluids_ , Vol. 32, No. 9, 2020, p. 096106.
  * Anderson et al. [2016] Anderson, D., Tannehill, J. C., and Pletcher, R. H., _Computational fluid mechanics and heat transfer_ , CRC Press, 2016.
  * Vinokur [1989] Vinokur, M., “An analysis of finite-difference and finite-volume formulations of conservation laws,” _Journal of computational physics_ , Vol. 81, No. 1, 1989, pp. 1–52.
  * Roe [1981] Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” _Journal of Computational Physics_ , Vol. 43, No. 2, 1981, pp. 357–372.
  * [32] Morrison, J. H., “Flux difference split scheme for turbulent transport equations. 1990,” _AIAA Paper_ , ????, pp. 90–5251.
  * Ranjan et al. [2020] Ranjan, R., Unnikrishnan, S., and Gaitonde, D., “A robust approach for stability analysis of complex flows using high-order Navier-Stokes solvers,” _Journal of Computational Physics_ , Vol. 403, 2020, p. 109076.
  * Bouhadji and Braza [2003a] Bouhadji, A., and Braza, M., “Organised modes and shock–vortex interaction in unsteady viscous transonic flows around an aerofoil: Part I: Mach number effect,” _Computers & fluids_, Vol. 32, No. 9, 2003a, pp. 1233–1260.
  * Bouhadji and Braza [2003b] Bouhadji, A., and Braza, M., “Organised modes and shock–vortex interaction in unsteady viscous transonic flows around an aerofoil: Part II: Reynolds number effect,” _Computers & fluids_, Vol. 32, No. 9, 2003b, pp. 1261–1281.
  * Harris [1981] Harris, C. D., “Two-dimensional aerodynamic characteristics of the NACA 0012 airfoil in the Langley 8 foot transonic pressure tunnel,” 1981.
  * Whitfield [1984] Whitfield, D., “Three-dimensional unsteady Euler equations solution using flux vector splitting,” _17th Fluid Dynamics, Plasma Dynamics, and Lasers Conference_ , 1984, p. 1552.
  * Sirovich [1987] Sirovich, L., “Turbulence and the dynamics of coherent structures. I. Coherent structures,” _Quarterly of applied mathematics_ , Vol. 45, No. 3, 1987, pp. 561–571.
  * Bergmann and Brancher [2004] Bergmann, L. C., and Brancher, J., “Optimal rotary control of the cylinder wake using POD reduced order model. AIAA Paper 2004-2323,” _2nd AIAA Flow Control Conference_ , 2004.
  * Iollo et al. [2000] Iollo, A., Lanteri, S., and Désidéri, J.-A., “Stability properties of POD–Galerkin approximations for the compressible Navier–Stokes equations,” _Theoretical and Computational Fluid Dynamics_ , Vol. 13, No. 6, 2000, pp. 377–396.
  * Balajewicz and Dowell [2012] Balajewicz, M., and Dowell, E. H., “Stabilization of projection-based reduced order models of the Navier–Stokes,” _Nonlinear Dynamics_ , Vol. 70, No. 2, 2012, pp. 1619–1632.
  * Amsallem and Farhat [2012] Amsallem, D., and Farhat, C., “Stabilization of projection-based reduced-order models,” _International Journal for Numerical Methods in Engineering_ , Vol. 91, No. 4, 2012, pp. 358–377.
  * Akhtar et al. [2009] Akhtar, I., Nayfeh, A. H., and Ribbens, C. J., “On the stability and extension of reduced-order Galerkin models in incompressible flows,” _Theoretical and Computational Fluid Dynamics_ , Vol. 23, No. 3, 2009, pp. 213–237.
  * Nagarajan et al. [2013] Nagarajan, K. K., Cordier, L., and Airiau, C., “Development and application of a reduced order model for the control of self-sustained instabilities in cavity flows,” _Communications in Computational Physics_ , Vol. 14, 2013, pp. pp–186.
  * Narkhede et al. [2022] Narkhede, M. V., Bartakke, P. P., and Sutaone, M. S., “A review on weight initialization strategies for neural networks,” _Artificial intelligence review_ , Vol. 55, No. 1, 2022, pp. 291–322.
  * Shen et al. [2021] Shen, S.-L., Atangana Njock, P. G., Zhou, A., and Lyu, H.-M., “Dynamic prediction of jet grouted column diameter in soft soil using Bi-LSTM deep learning,” _Acta Geotechnica_ , Vol. 16, No. 1, 2021, pp. 303–315.
  * Srinivasan et al. [2019] Srinivasan, P. A., Guastoni, L., Azizpour, H., Schlatter, P., and Vinuesa, R., “Predictions of turbulent shear flows using deep neural networks,” _Physical Review Fluids_ , Vol. 4, No. 5, 2019, p. 054603.
  * Mohan et al. [2019] Mohan, A., Daniel, D., Chertkov, M., and Livescu, D., “Compressed convolutional LSTM: An efficient deep learning framework to model high fidelity 3D turbulence,” _arXiv preprint arXiv:1903.00033_ , 2019.
  * Mohan et al. [2020] Mohan, A. T., Tretiak, D., Chertkov, M., and Livescu, D., “Spatio-temporal deep learning models of 3D turbulence with physics informed diagnostics,” _Journal of Turbulence_ , Vol. 21, No. 9-10, 2020, pp. 484–524.