# Koopman-Based Approach to Non-intrusive Projection-Based Reduced-Order
Modeling with Black-Box High-Fidelity Models. Part II: Application
111Submitted to the editors

S. Ashwin Renganathan 222School of Aerospace Engineering Georgia Institute of
Technology
275 Ferst Drive NW, Atlanta GA 30332

###### Abstract

A methodology for non-intrusive, projection-based non-linear model reduction
originally presented by Renganathan et. al. (2018) [1] is further extended
towards parametric systems with focus on application to aerospace design.
Specifically, we extend the method for static systems with parametric geometry
(that deforms the mesh), in addition to parametric boundary conditions. The
main idea is to first perform a transformation on the governing equations such
that it is lifted to a higher dimensional but linear under-determined system.
This enables one to extract the system matrices easily compared to that of the
original non-linear system. The under-determined system is closed with a set
of model-dependent non-linear constraints upon which the model reduction is
finally performed. The methodology is validated on the subsonic and transonic
inviscid flow past the NACA0012 and the RAE2822 airfoils. We further
demonstrate the utility of the approach by applying it to two common problems
in aerospace design namely, derivative-free global optimization and parametric
uncertainty quantification with Monte Carlo sampling. Overall, the methodology
is shown to achieve accuracy upto 5% and computational speed-up of 2-3 orders
of magnitude as that of the full-order model. Comparison against another non-
intrusive model reduction method revealed that the proposed approach is more
robust, accurate and retains the consistency between the state variables.

## Nomenclature

$N$ = | degrees of freedom in the full model
---|---
$M$ = | number of flow snapshots of full model (varying parameters)
$\mathbf{N}$ = | non-linear operator
$\mathbf{L}$ = | linear operator
$R$ = | residual operator of full model
$\mathbf{u}$ = | state of the full model
$\tilde{\mathbf{u}}$ = | reduced state
$\mathbf{U}$ = | matrix with stacked $\mathbf{u}$ snapshots
$\mathbf{y}$ = | observable ($=g(\mathbf{u})$)
$\tilde{\mathbf{y}}$ = | reduced observable
$\mathbf{\theta}$ = | parameters of the FOM
$h$ = | equality constraint functions
$S$ = | number of PDEs in the full model
$O$ = | number of terms in full model operated by a differential term
$\mathbf{V},\mathbf{\Sigma},\mathbf{W}$ = | SVD matrices
$\mathbf{\Phi}$ & $\mathbf{\Phi}_{i}$ = | overall and observable-wise POD matrices
$k$ & $k_{i}$ = | reduced size of all and individual observables
$\mathcal{M}$ = | manifold
$P_{\infty},\rho_{\infty},a_{\infty},\mu_{\infty},\mathcal{M},\alpha$ = | free-stream parameters
$\rho,p,\gamma$ = | density, pressure and specific-heat ratio
$u,v$ = |  $x$ and $y$ velocity components
$H,E$ = | enthalpy and internal energy
$\nabla$ = | gradient operator
$C_{P},C_{d},C_{l}$ = | coefficient of pressure, drag and lift of an airfoil
$\mathcal{N}$ = | Gaussian distribuion
$\\{\ell,\mu,\sigma^{2}\\}$ = | hyperparameters of the Gaussian process model
$\mathbf{R}$ = | correlation matrix of Gaussian process

##  1 Introduction

In the design of complex aerospace engineering systems such as aircraft,
rotorcraft and spacecraft, the advent of Computer Aided Engineering (CAE) and
high performance computing has significantly contributed in reducing the time
and cost involved in the design process. Specifically, high-fidelity
mathematical models such as Computational Fluid Dynamics (CFD), have to a good
extent substituted expensive physical testing with computer simulations [2].
However when it comes to decision making in the design of aerospace systems,
such high-fidelity models have to be queried several thousands of times, which
is not practical since each query could take few days to weeks to compute even
on a supercomputer. Specifically, we focus on _real-time decision making_ via
solving _many-query_ problems involving _global optimization_ and _uncertainty
quantification_.

Partial Differential Equations (PDE) based models such as CFD, typically
require an iterative solution whose computational cost scales in polynomial
time with the number of degrees of freedom ($N$) in the spatial
discretization. Practical problems of aerospace interest, typically involve
$N\sim\mathcal{O}(10^{6})$ degrees of freedom, and therefore their application
in the many-query context is not feasible. Therefore a _surrogate model_ that
would trade a small amount of accuracy for a significantly large gain in
computational cost is needed. Such a model would enable reliable, real-time
decision making thereby offering a necessary paradigm shift in the existing
process in aerospace design.

We specifically focus on Reduced Order Modeling (ROM) or Model Order Reduction
(MOR) which retain the underlying physical characteristics of the high-
fidelity model (henceforth referred as Full Order Model (FOM)) by projecting
the actual governing equations onto a suitably chosen low-dimensional
subspace, for which Proper Orthogonal Decomposition (POD) [3] is a common
technique. ROMs are quite popular for this reason, however are limited to
situations only where there is access to the discrete-form of governing
equations of the FOM. In situations where the governing equations are
available as a black-box (such as in commercial codes), projection-based ROM
is not feasible and hence a _non-intrusive_ technique is necessary. In such
situations, the common approach taken is to circumvent the projection step all
together and use the POD to directly approximate the state whose generalized
coordinates are interpolated in the time/parameter space (see [4, 5, 6, 7]).
Such an approach is effective in the sense that it is purely data-driven and
is more general in its applicability. However they do not guarantee that the
resulting ROM still satisfies the actual governing equations and the
associated consistency between the state variables which is an important
characteristic of MOR. For instance, independent surrogate models for the
thermodynamic variables pressure ($p$), density ($\rho$) and temperature ($T$)
might not satisfy the equation of state for an ideal gas $p=\rho RT$ (where
$R$ here is the specific gas constant). There has also been work done in data-
driven discovery of the governing equations from black-box codes, which can
then be used for model reduction [8, 9]; which rely on trajectory data of the
state in addition to initial and boundary conditions. However, in this work we
are focused on static parametric systems where there is no trajectory data and
further the boundary conditions might not be explicitly available. In an
earlier work, Renganathan et al (2018) [1] showed that by lifting the system
to a higher dimension via the _Koopman_ theory [10], a linear but under-
determined system can be obtained which can then be closed with a set of non-
linear problem-specific constraints. They further showed that with such a
technique the discrete linear operator can be extracted by discretizing the
linear differential terms via a method such as finite volume method, at a cost
that scales linearly with the grid size, $N$. However, the method was
restricted only to parameters in the boundary conditions. Here, we further
extend that method to apply towards systems with parametric geometry, where
the grid varies at every parameter snapshot. We show the effectiveness of the
approach by generating a database of ROMs for a pre-determined set of
snapshots and later interpolating between the ROMs. Therefore, this work
demonstrates the methodology developed in [1] with specific applications in
aerospace design by extending them to more general form of parametrizations.
The baseline method for comparison is chosen to be a _POD+interpolation_ where
as mentioned before, the generalized coordinates of the POD basis set are
directly interpolated in the parameter space. The author believes that this is
the only feasible non-intrusive method there exists in the literature given
the same scope of the present study and hence is chosen as the baseline
method. Further details are provided in section 5.1.

The rest of the paper is organized as follows. The non-intrusive MOR method is
first outlined in a more generic form in section 2, followed by the
compressible Euler equations in section 3. The model validation is shown for
subsonic and transonic inviscid flow past airfoils in section 4, to
demonstrate its predictive capability. The application of the method to two
many-query problems in aerospace design is discussed in section 5. The
conclusion section summarizes the main findings and outlines some directions
for future work.

The following notation is followed through the rest of the paper. Scalar
quantities are denoted by regular-face fonts in both upper and lower case.
Vectors and matrices are represented by lower-case bold-face and upper-case
bold-face fonts. The same rule applies to vectors and matrices made of several
vectors and block matrices respectively. Any exception to these conventions
are clarified as and when they are created in the paper.

##  2 Koopman-based Non-Intrusive Reduced Order Modeling

The Koopman theory forms the basis of the present approach where a linear
representation of the non-linear system is obtained. The finite volume method
is used to discretize linear differential terms as they are well suited for
unstructured computational grids and are almost the standard in commercial CFD
codes. As mentioned before, this approach depends on generating a database of
ROMs corresponding to a pre-determined set of parameter snapshots which are
then interpolated for new realizations of the parameters outside of the
training set. The ROM interpolation draws from differential geometry in order
to address the manifold-embedding of the ROM system matrices. Such an approach
is essential in ensuring that the fundamental properties of the system
matrices are retained post-interpolation. An outline of the methodology is
provided here while the reader is referred to [11, 1] for more details.

Consider a static non-linear system representing the FOM and in its
discretized form

$\mathbf{N}(\mathbf{u})=0$ (1)

where $\mathbf{N}$ represents a non-linear operator acting on the state
variable $\mathbf{u}\in\mathbb{R}^{N}$; $N$ being the degrees of freedom of
the FOM. Let $g:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ be a function that
operates on the state (such as $g(\mathbf{u})=$
$\mathbf{u}^{2},\mathbf{u}^{3},e^{\mathbf{u}}$ etc.), then we state that

$\mathbf{N}(\mathbf{u})\equiv\mathbf{L}\left([g_{1}(\mathbf{u})^{\top},g_{2}(\mathbf{u})^{\top},\ldots,g_{O}(\mathbf{u})^{\top}]^{\top}\right)$
(2)

where, $\mathbf{L}$ is a linear operator acting on the the _lifted_ system
where the $g_{i}(\mathbf{u})$’s replace $\mathbf{u}$. We call each
$g_{i}(\mathbf{u})$ an _observable_ following the convention of other works on
the topic (particularly [12]) and the number of such observables $O$ in (2) is
dependent on the system under consideration as is illustrated in section 3. We
then decompose the linear operator in the above equation as

$\mathbf{L}[g_{1}(\mathbf{u})^{\top},g_{2}(\mathbf{u})^{\top},\ldots,g_{O}(\mathbf{u})^{\top}]^{\top}\approx\mathbf{A}[g_{1}(\mathbf{u})^{\top},g_{2}(\mathbf{u})^{\top},\ldots,g_{O}(\mathbf{u})^{\top}]^{\top}+\mathbf{b}_{a}=0$
(3)

which follows the discretization of linear PDEs where,
$\mathbf{b_{a}}\in\mathbb{R}^{ON\times 1}$ is the vector that arises due to
the discretization of boundary conditions in addition to lumping any source
terms and $\mathbf{A}\in\mathbb{R}^{N\times ON}$ is the matrix resulting from
the discretization of the differential terms of a linear PDE. Overall, the
parametric changes that deforms the mesh (such as geometry shape) are captured
in $\mathbf{A}$ whereas the rest (such as free-stream boundary conditions) are
captured in $\mathbf{b_{a}}$. Note that (3) is linear but under-determined
system and hence for uniqueness of the solution, constraints are added as
discussed in section 2.1. Finally, we re-write (3) by modifying the notation
as
$[g_{1}(\mathbf{u})^{\top},g_{2}(\mathbf{u})^{\top},\ldots,g_{O}(\mathbf{u})^{\top}]^{\top}\rightarrow[\mathbf{y}_{1}^{\top},\ldots,\mathbf{y}_{O}^{\top}]^{\top}=\mathbf{y}$
and $-\mathbf{b_{a}}\rightarrow\mathbf{f}$, leading to

$\mathbf{A}\mathbf{y}=\mathbf{f}$ (4)

Figure 1: Graphical depiction of the present methodology. The original non-
linear static system is transformed to an under-determined linear system with
closure [11]

and we work with the transformed version of the FOM in (4) to develop the ROM.
Such a transformation enables us to extract $\mathbf{A}$ [1] and furthermore,
makes the overall methodology amenable to parametric interpolation as will be
illustrated in section 2.2. The overall idea behind the lifting transformation
to the FOM is depicted in Figure 1 and the model reduction is performed on the
transformed equations (the right hand side of the figure), which is explained
in the following sub-section.

###  2.1 Model Order Reduction

The total number of observables, $O$ is essentially inifinite if one were
seeking a closed linear transformation of the non-linear FOM [10]. However we
seek to find a finite $O$ which might result in an un-closed linear
transformation which are then closed with a set of constraints. In this work,
$O$ represents the total number of terms in the FOM that are functions of the
state and are operated by a linear differential operator; each of the
$g_{i}(\mathbf{u})$’s is picked from knowledge of the FOM in its _continuous_
PDE form, as will be illustrated in section 3. For a FOM that is a system of
$S$ coupled PDEs, note that $O\geq S$ always and $O>S$ for a non-linear
system. Therefore the observables can be written

$\mathbf{y}=[\mathbf{y}^{\top}_{1},...,\mathbf{y}^{\top}_{S},\mathbf{y}^{\top}_{S+1},...,\mathbf{y}^{\top}_{O}]^{\top}$
(5)

Finally, to close the under-determined transformed system, we add algebraic
equations that establish the non-linear consistency relationship between the
observables and the state thereby providing closure. These constraints are of
the form

$h_{i}:=\mathbf{y}_{S+i}-f_{i}(\mathbf{y}_{1},...,\mathbf{y}_{S})=0,~{}~{}i=1,\ldots,O-S$
(6)

where, $f_{i}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is a function that
operates on the observables. Note that for a system of $S$ coupled PDEs, $O-S$
constraints are needed to be specified to achieve closure. Although, there is
no unique way of specifying these constraints, we provide some guidelines in
section 2.1. Equation (6) along with (4) together form a closed system upon
which model reduction is performed. To perform the projection, the truncated
basis set for each observable $\mathbf{y}_{i}$ is extracted by performing POD
on the snapshot matrix of $\mathbf{y}_{i}$ (generated by running the FOM at
$M$ distinct parameter points) and are denoted $\Phi_{i}\in\mathbb{R}^{N\times
k_{i}}$, i.e.

$\mathbf{Y}_{i}=\begin{bmatrix}\vdots&\vdots&\vdots&&\vdots\\\
\mathbf{y}_{i}^{(1)}&\mathbf{y}_{i}^{(2)}&\mathbf{y}_{i}^{(3)}&\ldots&\mathbf{y}_{i}^{(M)}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ \end{bmatrix}\underset{\text{thin-
svd}}{=}\mathbf{V}_{i}\mathbf{\Sigma}_{i}\mathbf{W}_{i}^{\top}$ (7)

where $\mathbf{y}_{i}^{(j)}$ is the $jth$ snapshot of observable
$\mathbf{y}_{i}$, $\Phi_{i}$ is the first $k_{i}$ columns of $\mathbf{V}_{i}$,
$M$ is the total number of snapshots used for model training and the POD basis
vectors are obtained from the thin-svd decomposition of $\mathbf{Y}_{i}$. This
leads to the trial basis matrix for the overall system defined as a block-
diagonal matrix of all the $O$ POD basis set given below

$\mathbf{\Phi}=\text{blkdiag}\\{\Phi_{1},\ldots,\Phi_{O}\\}~{}\in\mathbb{R}^{ON\times
k}$ (8)

where $k=k_{1}+\ldots+k_{O}$. The _reduced_ observable is then given by
$\tilde{\mathbf{y}}_{i}\approx\Phi_{i}^{\top}\mathbf{y}_{i}$. Recall that
$\mathbf{A}$ is non-square since it is $\in\mathbb{R}^{N\times ON}$, and hence
a suitable choice for the test basis for projection is
$\mathbf{\Psi}=\mathbf{A}\mathbf{\Phi}$. Note that this choice of the test
basis is equivalent to a galerkin projection ($\mathbf{\Psi}=\mathbf{\Phi}$)
on the normal equations. i.e. on
$\mathbf{A}^{\top}\mathbf{A}\mathbf{y}=\mathbf{A}^{\top}\mathbf{f}$. Let
$\mathbf{B}=\mathbf{A}^{\top}\mathbf{A}\in\mathbb{R}^{ON\times ON}$; then the
projection leads to

$\mathbf{\Phi}^{\top}\mathbf{B}\mathbf{\Phi}\tilde{\mathbf{y}}=\mathbf{\Phi}^{\top}\mathbf{A}^{\top}\mathbf{f}$
(9)

Setting
$\tilde{\mathbf{f}}=\mathbf{\Phi}^{\top}\mathbf{A}^{\top}\mathbf{f}\in\mathbb{R}^{k}$
and
$\tilde{\mathbf{B}}=\mathbf{\Phi}^{\top}\mathbf{B}\mathbf{\Phi}\in\mathbb{R}^{k\times
k}$, this leads to the reduced order model

$\tilde{\mathbf{B}}\tilde{\mathbf{y}}=\tilde{\mathbf{f}}$ (10)

The ROM given by Equation 10 is now a $k\times k$ system where $k<<N$ and is
solved along with the constraints presented in Equation 6, posed as a non-
linear program as shown below

$\displaystyle\underset{\tilde{\mathbf{y}}}{\text{minimize}}$
$\displaystyle\frac{1}{2}\|\tilde{\mathbf{B}}\tilde{\mathbf{y}}-\tilde{\mathbf{f}}\|_{2}^{2}$
(11) s.t.
$\displaystyle\mathbf{\Phi}^{\top}h(\mathbf{\Phi}\tilde{\mathbf{y}})=0$

The main hypothesis of this work is that the ROM given by Equation 11 still
approximately satisfies the governing equations and this is verified in the
Section 4. The optimization problem in Equation 11 needs special treatment to
handle the non-linear constraints which still depend on the full state of
observables, and is efficiently done using the DEIM [13]; see Appendix 7.3 for
details on implementation for a specific example. The ROM in (11) is solved
via Sequential Quadratic Programming (SQP) [14] with the objective function
and constraint tolerances set to $10^{-6}$ and the number of function
evaluations bounded by $4\times 10^{6}$. The initial guess to the solution of
(11) is given as the nearest flow snapshot to the test parameter at which
prediction is sought.

###  2.2 ROM Interpolation

The proposed approach that leads to the ROM in the form of (10) corresponds to
one parameter snapshot since $\tilde{\mathbf{B}}$ and $\tilde{\mathbf{f}}$ are
parameter dependent. Therefore, the approach generates a database of ROMs for
a pre-determined set of parameter snapshots, which are later interpolated to
predict the state at a new parameter. The interpolation is carried out in a
manner that retains the inherent structure and properties of the matrix
$\mathbf{\tilde{B}}$ post-interpolation. The general principle that is
followed is to map the matrices to a plane that is locally tangent to the
manifold in which they are originally embedded. The anchor point on the
manifold, which is the point of tangency is chosen to be one of the matrices
themselves. While this choice is arbitrary, in this work we use the matrix
that corresponds to the nearest (in the _standardized_ Euclidean sense)
parameter snapshot to the test parameter. The traditional Euclidean space
interpolation (where typical vector operations are valid) is then carried out
in the tangent plane after which they are mapped back to the manifold. The
mapping to and from the tangent plane are carried out via logarithmic and
exponential relationships as depicted in Figure 2.

Figure 2: A graphical representation of a manifold $\mathcal{M}$ and the
embedding of parametric matrices $\tilde{\mathbf{B}}(\boldsymbol{\theta})$. A
direct element-wise interpolation of $\tilde{\mathbf{B}}$ at
$\hat{\boldsymbol{\theta}}$ may not necessarily result in a matrix
$\in\mathcal{M}$, and is carried out after mapping to the tangent plane,
$\mathcal{T}(\mathcal{M})$

The matrix
$\tilde{\mathbf{B}}=\mathbf{\Phi}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{\Phi}$
in Equation 11 is symmetric positive definite (SPD) for the following reasons.
Firstly, $\mathbf{A}^{\top}\mathbf{A}$ is a co-variance matrix and hence is
symmetric positive semi-definite (see [15], sec. 5.3 ). Furthermore,
multiplication by orthogonal matrix $\mathbf{\Phi}$ of rank $k$ where
$k<rank(A)$ ensures $\tilde{\mathbf{B}}$ is SPD. SPD matrices of size $k\times
k$ form a special group called the $SPD(k)$ [16, 17] and the manifold they are
embedded in is denoted as $\mathcal{M}$. Also, for the set of all SPD matrices
$\mathbf{B}~{}\in\mathcal{M}$, the tangent plane is the set of all symmetric
matrices, $\mathbf{B^{\prime}}$ [17]. Any metric ($\mathcal{M}_{f}$) defined
on $SPD(k)$ for any two matrices uses the following functional relationship

$\mathcal{M}_{f}(\mathbf{B}_{1},\mathbf{B}_{2})=\mathbf{B}_{1}^{1/2}f\left(\mathbf{B}_{1}^{-1/2}\mathbf{B}_{2}\mathbf{B}_{1}^{-1/2}\right)\mathbf{B}_{1}^{1/2}$
(12)

which leads to the following results for the exponential and logarithmic
mapping for $SPD(k)$ [16] where, $\tilde{\mathbf{B}}_{0}$ is the anchor point
and $\tilde{\mathbf{B}^{\prime}}$ is the point whose mapping is desired. The
exponential mapping of $\tilde{\mathbf{B}}^{\prime}$ from tangent plane to
$\mathcal{M}$ at $\tilde{\mathbf{B}}_{0}\in\mathcal{M}$, to $\mathcal{M}$ is
given by

$\mathbf{Exp}_{\tilde{\mathbf{B}}_{0}}\tilde{\mathbf{B}}=\tilde{\mathbf{B}}_{0}^{1/2}\left(\tilde{\mathbf{B}}_{0}^{-1/2}exp(\tilde{\mathbf{B}}^{\prime})\tilde{\mathbf{B}}_{0}^{-1/2}\right)\tilde{\mathbf{B}}_{0}^{1/2}$
(13)

and the logarithmic mapping of $\tilde{\mathbf{B}}\in\mathcal{M}$ to tangent
plane to $\mathcal{M}$ at $\tilde{\mathbf{B}}_{0}\in\mathcal{M}$:

$\mathbf{Log}_{\tilde{\mathbf{B}}_{0}}\tilde{\mathbf{B}^{\prime}}=\tilde{\mathbf{B}}_{0}^{1/2}log\left(\tilde{\mathbf{B}}_{0}^{-1/2}\tilde{\mathbf{B}}\tilde{\mathbf{B}}_{0}^{-1/2}\right)\tilde{\mathbf{B}}_{0}^{1/2}$
(14)

The results presented in Equations 13 and 14 are used in this work to perform
mapping to and from the tangent space. Once on the tangent space the matrices
are interpolated element-wise using multivariate Lagrange polynomials of
$2^{nd}$ order. The overall method is summarized in Figure 3.

Figure 3: Flowchart representation of the overall framework

##  3 Governing Equations & Experimental Setup

We demonstrate the methodology on the inviscid flow past airfoils, for which
the NACA0012 and RAE2822 are chosen as baselines at subsonic and transonic
flow regimes respectively. We begin by first illustrating the method outlined
in section 2 on the compressible Euler equations and then discussing the
parametrization behind the chosen test cases.

Table 2: Free-stream conditions for the NACA and RAE test cases | NACA | RAE
---|---|---
$p_{\infty}$ | 101,325 | 28,745 $Pa$
$\rho_{\infty}$ | 1.225 | 0.44 $kg/m^{3}$
$a_{\infty}$ | 340.296 | 301.86 $m/s$
$\mu_{\infty}$ | 1.78E-5 | 1.49E-5 $Pa-s$
$\mathbb{M}$ | 0.60 | 0.734
$\alpha$ | 2.0 | 2.79 $deg.$

###  3.1 Governing Equations

The Euler equations governing the 2D, compressible, inviscid flow past an
airfoil are the governing equations on which we perform model reduction. This
coupled non-linear system is solved via the commercial black-box CFD solver,
STARCCM+ [18]. The equations in conservation form is provided in Eq.15 below,

$\nabla_{x}\mathbf{F}+\nabla_{y}\mathbf{G}=0$ (15)

where

$\displaystyle\mathbf{F}=\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho uv\\\
\rho uH\end{bmatrix},~{}\mathbf{G}=\begin{bmatrix}\rho v\\\ \rho uv\\\ \rho
v^{2}+p\\\ \rho vH\end{bmatrix}$ $\displaystyle H=E+\frac{p}{\rho}$
$\displaystyle\rho E=\frac{1}{2}\rho(u^{2}+v^{2})+\frac{p}{\gamma-1}$

and $\nabla_{x}$ and $\nabla_{y}$ are the $x$ and $y$ components of the
gradient operator $\nabla$ respectively. The following transformation is then
performed

$[\rho u,\rho v,\rho uv,p,\rho u^{2},\rho v^{2},\rho uH,\rho
vH]^{\top}\rightarrow[y_{1},y_{2},y_{3},y_{4},y_{5},y_{6},y_{7},y_{8}]^{\top}$

from the state variables to observables, leading to the lifted model.

$\begin{bmatrix}\nabla_{x}&\nabla_{y}&&&&&&\\\
&&\nabla_{y}&\nabla_{x}&\nabla_{x}&&&\\\
&&\nabla_{x}&\nabla_{y}&&\nabla_{y}&&\\\ &&&&&&\nabla_{x}&\nabla_{y}\\\
\end{bmatrix}\begin{bmatrix}y_{1}\\\ y_{2}\\\ y_{3}\\\ y_{4}\\\ y_{5}\\\
y_{6}\\\ y_{7}\\\ y_{8}\end{bmatrix}=\mathbf{0}$ (16)

where in the above equation, empty spaces in the matrix denote zeros. The
equation upon discretization leads to

$\underbrace{\begin{bmatrix}\mathbf{G}_{x}&\mathbf{G}_{y}&&&&&&\\\
&&\mathbf{G}_{y}&\mathbf{G}_{x}&\mathbf{G}_{x}&&&\\\
&&\mathbf{G}_{x}&\mathbf{G}_{y}&&\mathbf{G}_{y}&&\\\
&&&&&&\mathbf{G}_{x}&\mathbf{G}_{y}\\\
\end{bmatrix}}_{\mathbf{A}}\underbrace{\begin{bmatrix}\mathbf{y}_{1}\\\
\mathbf{y}_{2}\\\ \mathbf{y}_{3}\\\ \mathbf{y}_{4}\\\ \mathbf{y}_{5}\\\
\mathbf{y}_{6}\\\ \mathbf{y}_{7}\\\
\mathbf{y}_{8}\end{bmatrix}}_{\mathbf{y}}=\underbrace{-\begin{bmatrix}\mathbf{b}_{a1}\\\
\mathbf{b}_{a2}\\\ \mathbf{b}_{a3}\\\ \mathbf{b}_{a4}\\\ \mathbf{b}_{a5}\\\
\mathbf{b}_{a6}\\\ \mathbf{b}_{a7}\\\
\mathbf{b}_{a8}\end{bmatrix}}_{\mathbf{f}}$ (17)

where, $\mathbf{G}_{x}$ and $\mathbf{G}_{y}$ represents the discrete version
of the gradient operators $\nabla_{x}$ and $\nabla_{y}$ respectively and the
empty spaces denote block matrices of zeros. The parameter-dependent
$\mathbf{A}$ matrix is obtained directly by discretizing the linear terms
$\mathbf{G}_{x}$ and $\mathbf{G}_{y}$ via finite volume method. The grid is
exported in the CFD General Notation System (CGNS) [19] for this purpose, as
mentioned in Figure 3. The snapshots $\mathbf{y}$ are applied to the
$\mathbf{A}$ matrix and the RHS $\mathbf{f}$ is extracted for each parameter.
With the FOM reduced to the $\mathbf{A}\mathbf{y}=\mathbf{f}$ form and
$\mathbf{A}\in\mathbb{R}^{4N\times
8N},~{}\mathbf{y},\mathbf{f}\in\mathbb{R}^{8N}$, Equation 17 represents an
under-determined system. Therefore they are closed using non-linear
constraints given by Equation 18. Notice that the constraints express the
relationship between the first $S=4$ observables ($y_{1}$ through $y_{4}$) and
the remaining $O-S;~{}(O=8)$ observables ($y_{5}$ through $y_{8}$). As
mentioned in the previous section, the choice of the first $S$ observables and
hence the $O-S$ constraints is non-unique. However, from experience trying out
different choices in this work, it is found that the following heuristics
ensure a stable transformation from the observables back to the state: (i) the
terms starting from the lowest order are picked as the first $S$ observables
($\rho u,\rho v,\rho uv,p$ in this case) and (ii) one of the first $S$
observables is set to be a state variable ($y_{4}=p$ in this case). It should
be noted that all the observables that are in excess of the number of
equations can be expressed as some function of the rest because the number of
independent observables are only as many as the number of PDEs in the FOM
($S=4$). The constraints are expressed in terms of the continuous form of the
state and observable as follows

$\displaystyle h_{1}=\rho u^{2}-\frac{(\rho u)(\rho uv)}{\rho v}\equiv
y_{5}-\frac{y_{1}y_{3}}{y_{2}}=0$ (18) $\displaystyle h_{2}=\rho
v^{2}-\frac{(\rho v)(\rho uv)}{\rho u}\equiv y_{6}-\frac{y_{2}y_{3}}{y_{1}}=0$
$\displaystyle h_{3}=\rho uH-\rho u\left(E+\frac{p}{\rho}\right)\equiv
y_{7}-y_{1}\left(E+\frac{y_{4}y_{3}}{y_{1}y_{2}}\right)=0$ $\displaystyle
h_{4}=\rho vH-\rho v\left(E+\frac{p}{\rho}\right)\equiv
y_{8}-y_{2}\left(E+\frac{y_{4}y_{3}}{y_{1}y_{2}}\right)=0$

With discrete observables ($\mathbf{y}_{i}$) all the operations in (15) are
performed element-wise. In all the results presented in the following
sections, the POD modes constituting $99.99\%$ of the energy (cumulative
fraction sum of the singular values in $\mathbf{\Sigma}_{i}$) are retained.
The error metrics used for all the results are defined as follows

$\begin{split}C_{P}~{}\text{Error}&=\frac{\left\|C_{P}^{FOM}-C_{P}^{ROM}\right\|_{\infty}}{\left\|C_{P}^{FOM}\right\|_{\infty}}\times
100\\\
C_{d}~{}\text{Error}&=\frac{\left|C_{d}^{FOM}-C_{d}^{ROM}\right|}{C_{d}^{FOM}}\times
100\\\
C_{l}~{}\text{Error}&=\frac{\left|C_{l}^{FOM}-C_{l}^{ROM}\right|}{C_{l}^{FOM}}\times
100\end{split}$ (19)

###  3.2 Test cases and parametrization

The baseline shapes are parameterized using Class Shape Transformation (CST)
[20, 21]. The CST model of parametrization defines a class function $c$ and a
shape function $s$ and the curve being parameterized is specified as their
product. The main idea is that the class function serves to define a general
class of geometry such as airfoils, missiles or sears-haack body, while the
shape function serves to define the unique shape within a particular class of
shapes (such as a NACA0012 vs RAE2822 airfoil). The class function, $c(\psi)$
is more generally defined as

$c_{n_{1}}^{n_{2}}(\psi):=\psi^{n_{1}}(1-\psi)^{n_{2}}$ (20)

where the variable $\psi$ represents the non-dimensional chord-wise distance.
$n_{1}$ and $n_{2}$ define the specific class; for instance
$n_{1}=0.5,~{}n_{2}=1$ and hence $\sqrt{\psi}(1-\psi)$ defines airfoils with
rounded leading edge and a sharp trailing edge [20]. The unique shape of an
airfoil is driven by the shape function, specified as follows

$s(\psi)=\sum_{i=0}^{n}A_{i}\psi^{i}$ (21)

where $A_{i}$ are the coefficients which are also the shape parameters. The
NACA0012 and RAE2822 are parameterized using 6 and 8 variables respectively,
whose values are given by $A_{NACA0012}$ and $A_{RAE2822}$ where the top and
bottom rows correspond to the upper and lower surfaces of the airfoils.
Further details of the parameterization are provided in Appendix 7.1.

$\displaystyle A_{NACA0012}$
$\displaystyle=\begin{bmatrix}0.1689&0.2699&0.1387\\\
-0.1689&-0.2699&-0.1387\end{bmatrix}$ (22) $\displaystyle A_{RAE2822}$
$\displaystyle=\begin{bmatrix}0.1268&0.4670&0.5834&0.2103\\\
-0.1268&-0.5425&-0.5096&0.0581\end{bmatrix}$

(a) RAE2822

(b) NACA0012

Figure 4: Comparison of the CST approximation against the true curve

##  4 Model Validation

### NACA0012

The CST coefficients representing the NACA0012 baseline is perturbed $\pm$ 30%
to generate new airfoil shapes, a sample of which is shown in Figure 5. A
total of 170 such points were generated using a Latin Hypercube design,
$M=$160 of which was used in model building while the remaining was used to
validate the model. The system matrix $\tilde{\mathbf{B}}$ is interpolated in
the tangent space to the manifold of symmetric positive definite matrices, as
explained in Section 2.2, while the RHS $\tilde{\mathbf{f}}$ is interpolated
in the Euclidean space. In both cases a multi-variate polynomial in the
Lagrange form is used for interpolation as mentioned previously (see
Renganathan(2018) [11], Algorithm 1 for details).

Figure 5: Family of airfoils generated by perturbing (by $\pm 30\%$) the CST
coefficients of the NACA0012 baseline

The comparison of the ROM predicted pressure coefficient on the airfoil
surface against the FOM solution, for a select two of the validation set is
shown in Figure 6. Across all the 10 validation cases, the maximum and average
errors in $C_{P}$ are $\approx$ 5% and $\approx$ 2% respectively, whereas the
computational speedup is $\approx$ 100x. Therefore, the effectiveness of
proposed approach is well established. Similar accuracy was observed with the
lift coefficient, which are tabulated in Appendix 7.2. The overall flow-field
comparison in terms of overlaid pressure and mach number contours is also
shown in Figure 6 which further re-iterate the low prediction errors observed.

---
(a) Validation Case-1 (b) Validation Case-4 (c) Validation Case-1 (d)
Validation Case-4

Figure 6: Validation of the subsonic NACA0012 test case. Top row shows airfoil
shapes and $C_{P}$ distributions, bottom row shows Mach number and absolute
pressure contours.

### RAE2822

Now we demonstrate the method in the transonic regime using the RAE2822 test
case. Similar to the NACA test case, CST coefficients representing the
baseline are perturbed $\pm$ 30 % to generate new airfoil shapes, a sample of
which is shown in Figure 7. A total $M=$160 snapshots were again used in model
building. However, the freestream mach number for this case is set as
$\mathbb{M}=0.734$ which leads to a shock whose strength and location are
affected by perturbing the shape CST coefficients. The rest of the freestream
conditions are summarized in Table 2

Figure 7: Family of airfoils generated by perturbing (by $\pm 30\%$) the CST
coefficients of the RAE2822 baseline

We begin by comparing the ROM predictions against FOM solutions in terms of
the coefficient of pressure distributions in Figure 8. Two specific cases are
showcased in the figure to illustrate the strength and weakness of the
approach in transonic regime. Overall, from all the validation cases, it is
observed that the ROM does capture the shock location within a maximum error
of 5%. However, this translates in to an average $L_{2}$ error of 8.5% in
$C_{P}$ predictions and about 16% in the prediction of $C_{d}$. The average
error in $C_{l}$ predictions were much lower at 4%. Such a trend where the
lift coefficient were predicted with greater accuracy than drag was also
previously observed in [1]. The main reason for this is expected to be that in
the absence of viscous effects, the computation of the drag force shows
greater sensitivity to discrepancy in pressure. This is because unless at high
angles of attack, the projected area in the direction normal to the freestream
velocity is very small leading to larger sensitivities in $C_{d}$
computations. However, as will be demonstrated in section 5.3, the proposed
approach demonstrates more robustness in the drag prediction than competing
methods.

---
(a) Val. Case-5 (b) Val. Case-7 (c) Val. Case-5 (d) Val. Case-7

Figure 8: Comparison of pressure coefficient $C_{P}$ predicted by the ROM with
the true solution due to the FOM for various airfoil shapes that represent the
validation cases

### Discussion

Under subsonic shock-free conditions, the predictions are consistently under
5% whereas under transonic conditions, the predictive capability drops. In the
presence of discontinuities in the flow field such as moving shocks, a POD-
based method is unable to accurately predict shock location and strength. Such
a limitation has been reported by others in the literature who propose a
domain-decomposition method to isolate the shock-prone regions of the flow-
field where the FOM is solved whereas the ROM is used to approximate the rest
of the flow field [22, 23]. However such methods are intrusive in the sense
that they might require access to source code of the FOM to effectively
manipulate the flow domains which is not feasible with black-box codes.
Another method reported in the literature to address this problem is to do a
space-transformation of the snapshots containing moving discontinuities [24].
The goal of the proposed approach are first steps towards enabling projection-
based MOR for black-box static parametric systems although such enhancements
can be added to it to address specific problems. Additionally, the solution to
the ROM in (10) shows sensitivity to initial guess and hence can lead to a
local optimum which also contributes to the performance penalty of the
approach in transonic regime. Despite the limitations in the transonic regime,
it has been observed that the ROM still predicts the shock within 5% chord-
length variability. Furthermore, the $C_{P}$, $C_{d}$ and $C_{l}$ are
predicted with an average error of $14\%$, $16\%$ and $4\%$ respectively which
are still useful in the conceptual and preliminary stages of the design
process where some accuracy can be traded for larger gain in computational
costs for faster decision making.

##  5 Application: Many-Query Problems

We pick two problems which are pertinent to aerospace design and require
several queries to the model which is not practical if one had only the FOM.
The first problem is a non-convex optimization problem which greatly benefits
from derivative-free optimization techniques such as the Genetic Algorithm
(GA) [25] that efficiently balances exploration and exploitation to determine
the global optimum. The second problem concerns the approximation of
probability densities which similarly require several queries to the model.

As a baseline for comparison of the proposed approach, we use a more commonly
used non-intrusive technique which replaces the projection step in MOR with a
direct interpolation of the reduced state. As mentioned previously, such an
approach no longer guarantees that the ROM still satisfies the governing
equations of the FOM, but their accuracy has been found useful in multiple
studies [5, 26]. Here we replace the projection step with the Kriging [27, 28]
interpolator and use that for comparison against the proposed method. In the
rest of the section the proposed projection-based approach is referred as
_POD-Proj._ whereas the interpolation-based approach is referred as _POD-
Krig._

###  5.1 POD + Kriging

Recall from section 2 that the full and reduced observable are related via the
relationship
$\tilde{\mathbf{y}}_{i}\approx\mathbf{\Phi}^{\top}\mathbf{y}_{i}$. Here each
element of the reduced observable $\tilde{\mathbf{y}}_{i}$ is assumed to be a
smooth function of the parameters $\boldsymbol{\theta}$ and is interpolated in
the $\boldsymbol{\theta}$-space via Kriging. The Kriging approach assumes that
the true response is a relaization of a Gaussian process (GP) which is
completely specified by a _mean_ and _covariance function_ and estimates the
hyperparameters of the GP from observations at finite number of locations.
Kriging interpolates noise-free data, regresses noisy data and the prediction
at unknown sites is given by a Gaussian probability distribution as opposed to
a deterministic value. However it is customary to use the expected value of
the predictive distribution as a surrogate for the true function.

The $jth$ element of the $ith$ reduced observable is assumed to take the form

$\tilde{y}_{ij}=\mu+z(\boldsymbol{\theta})$ (23)

where $z(\boldsymbol{\theta})\sim\mathcal{N}(\mathbf{0},\mathbf{C})$ with
$\mathbf{0}$ being a vector of zeros of appropriate length, $\mathbf{C}$ the
covariance matrix and $\mu$ the unknown mean. The GP is assumed to be
statistically _stationary_ [29] with homoscedastic variance which simplifies
the covariance matrix as $\mathbf{C}_{ij}=\sigma^{2}\mathbf{R}_{ij}$, where
$\mathbf{R}$ is the _correlation_ matrix. The assumption of stationarity also
means that
$\mathbf{R}_{ij}=\mathbf{R}(\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|)$
where $\|\cdot\|$ denotes the Euclidean distance. The elements of the
correlation matrix are defined by a parametrized _kernel_ as
$\mathbf{R}_{ij}=r(\boldsymbol{\theta}_{i},\boldsymbol{\theta}_{j}~{};\ell)$.
The hyperparameters of the Kriging model are therefore
$\\{\ell,\mu,\sigma^{2}\\}$. In this work, we assume the stationary and
isotropic _squared-exponential_ kernel are given by the following equation

$\mathbf{R}_{ij}=r(\boldsymbol{\theta}_{i},\boldsymbol{\theta}_{j};\ell)=\text{exp}\left(-\frac{\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|_{2}^{2}}{2\ell^{2}}\right)$
(24)

The training data for each $\tilde{y}_{ij}$ is the $jth$ row of the product
$\mathbf{\Sigma}_{i}\times\mathbf{W}_{i}^{\top}$ in the thin-svd step in (7);
let that be denoted as $\tilde{\mathbf{y}}^{j}$. Each $\tilde{\mathbf{y}}^{j}$
is a vector of length $M$ corresponding to the parameter snapshots
$\boldsymbol{\Theta}=[\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{M}]$.
Then the Kriging model assumes that the prediction at some unknown
$\boldsymbol{\theta}^{*}$ forms a joint normal distribution with
$\tilde{\mathbf{y}}^{j}$. That is

$\begin{bmatrix}\tilde{\mathbf{y}}^{j}\\\
\tilde{y}_{ij}^{*}\end{bmatrix}=\sim\mathcal{N}\left(\mathbf{0},\sigma^{2}\begin{bmatrix}\mathbf{R}&r(\boldsymbol{\Theta},\boldsymbol{\theta}^{*})\\\
r(\boldsymbol{\theta}^{*},\boldsymbol{\Theta})&1\end{bmatrix}\right)$ (25)

Conditioning the distribution of $\tilde{y}_{ij}^{*}$ in (25) on the
observations one then obtains the posterior predictive distribution [28] given
by

$\tilde{y}_{ij}*|\boldsymbol{\theta}^{*},\boldsymbol{\Theta},\tilde{\mathbf{y}}^{j}\sim\mathcal{N}\left(r(\boldsymbol{\theta}^{*},\boldsymbol{\Theta})\mathbf{R}^{-1}\tilde{\mathbf{y}}^{j},\sigma^{2}\left[1-r(\boldsymbol{\theta}^{*},\boldsymbol{\Theta})\mathbf{R}^{-1}r(\boldsymbol{\Theta},\boldsymbol{\theta}^{*})\right]\right)$
(26)

where
$r(\boldsymbol{\theta}^{*},\boldsymbol{\Theta})=r(\boldsymbol{\Theta},\boldsymbol{\theta}^{*})^{\top}$
is a column vector of the correlation between $\boldsymbol{\theta}^{*}$ and
$\boldsymbol{\Theta}$ and the unknown mean
$\mu=r(\boldsymbol{\theta}^{*},\boldsymbol{\Theta})\mathbf{R}^{-1}\tilde{\mathbf{y}}^{j}$
(which is also our Kriging predictor). The other two hyperparameters
($\ell,\sigma^{2}$) are determined by maximizing the log marginal likelihood
of $\tilde{y}_{ij}*$

$\text{log}~{}p\left(\tilde{\mathbf{y}}^{j}|\boldsymbol{\Theta},\ell,\sigma^{2}\right)=-\frac{1}{2}\tilde{\mathbf{y}}^{j\top}\mathbf{R}^{-1}\tilde{\mathbf{y}}^{j}-\frac{1}{2}\text{log}|\mathbf{R}|-\frac{M}{2}\text{log}~{}2\pi$
(27)

###  5.2 Inverse Design

In aerodynamic design, a specific aerodynamic load distribution about an
aerodynamic object is of interest. For instance, under incompressible flow
assmptions for a finite wing, an elliptic lift distribution along the wing is
known to produce the least induced drag [30]. Similarly, in the preliminary
design of propellers [31] and turbines[32], a certain lift distribution along
the blade is an input to the design process. In such cases the actual design
(shape) parameters that produce such a load distribution is of interest. We
call such a problem the _inverse design_. Here we fix the desired pressure
coefficient distribution as our _target_ and search the design space for the
shape parameters that would best approximate the target. Then the predicted
airfoil shape is compared against the true shape. Therefore we are interested
in solving the following optimization problem

$\displaystyle\underset{\boldsymbol{\theta}}{arg\text{min}}$
$\displaystyle~{}~{}\frac{1}{2}\|C_{P}(\mathbf{u},\boldsymbol{\theta})-C^{*}_{P}(\mathbf{u},\boldsymbol{\theta})\|_{2}^{2}$
(28) subject to: $\displaystyle R(\mathbf{u},\boldsymbol{\theta})=0$
$\displaystyle\boldsymbol{\theta}_{l}\leq\boldsymbol{\theta}\leq\boldsymbol{\theta}_{u}$

where $C^{*}_{P}$ is the target pressure distribution and $R$ is the residual
operator. Naturally we want to replace the full-order governing equations with
the ROM and hence we solve the modified problem

$\displaystyle\underset{\boldsymbol{\theta}}{arg\text{min}}$
$\displaystyle~{}~{}\frac{1}{2}\|C_{P}(\mathbf{u},\boldsymbol{\theta})-C^{*}_{P}(\mathbf{u},\boldsymbol{\theta})\|_{2}^{2}$
(29) subject to:
$\displaystyle\mathbf{\Psi}^{\top}R(\mathbf{\Phi}\tilde{\mathbf{u}},\boldsymbol{\theta})=0$
$\displaystyle\boldsymbol{\theta}_{l}\leq\boldsymbol{\theta}\leq\boldsymbol{\theta}_{u}$

Table 3: Free-stream conditions for inverse design $P_{\infty}$ | 101325 Pa
---|---
$T_{\infty}$ | 288 K
$\rho_{\infty}$ | 1.225 $kg/m^{3}$
$a_{\infty}$ | 340.296 m/s
$M_{\infty}$ | 0.6
$\alpha$ | 2 ∘

The free-stream conditions used for this test case are summarized in Table 3.
The optimum shape was searched using a GA based optimizer with a population
size of 30 per generation and a total of 60 generations. The constraint and
function convergence tolerance were set to $10^{-5}$ and $10^{-3}$
respectively and the optimization required a total of 1830 function
evaluations to determine the final design. The final design is shown in Figure
9. The predictions by the proposed approach appears quite similar to the POD-
Krig. However, the difference is more noticeable in the GA optimizer
convergence history shown in Figure 10 where the POD-Krig. leads to a sub-
optimal design compared to the proposed approach. In this specific example the
proposed approach only marginally outperforms the baseline approach, but this
is put to test further in predicting the lift and drag coefficients section
5.3. It was earlier shown in [1] that in the transonic regime a POD-Krig. like
approach is prone to predicting non-physical shock patterns mainly because the
reduced state is not necessarily smoothly varying. However, such an approach
is still effective if the flow does not contain parameter-dependent
discontinuities as demonstrated in this example. The main benefit of the
proposed approach is the gain in computational times. It required
approximately 3.7 hrs of wall-clock time for the inverse design problem
whereas the equivalent FOM wall-clock time for the same number of function
evaluations is expected to take roughly 300 hrs.

(a) Airfoil shape comparison

(b) $C_{P}$ distribution comparison

(c) Airfoil shape comparison

(d) $C_{P}$ distribution comparison

Figure 9: Comparison of predicted-target design with the target. Top row: POD-
Proj., bottom row: POD-Krig. Figure 10: GA optimizer convergence for inverse
design problem

###  5.3 Uncertainty Quantification

Next, we demonstrate the present methodology towards the uncertainty
quantification problem. Uncertainty in aerospace design is inevitable and
manifests itself either due to lack of knowledge (ex: biased models) or
variability inherent in the system (ex: noisy manufacturing process). The
aerospace design process should therefore account for the uncertainties in
order to make reliable decisions early on in the design process. Specifically,
we are interested in quantifying the uncertainty in the airfoil lift and drag
coefficients due to the manufacturing process induced variations in the
airfoil shape. We use the same $\pm~{}30\%$ variation in the airfoil CST
coefficients and uniformly sample this design space. A Monte Carlo (MC) [33]
simulation is carried out to propagate the input uncertainty into the model to
quantify the uncertainties in the quantities of interest - namely, the lift
and drag coefficients. We use the NACA test case under operating conditions
defined in Table 3.

A total of 4000 samples were uniformly sampled from the $\pm~{}30\%$ range on
the CST coefficients on which the proposed approach (POD-Proj.) and the
baseline approach (POD-Krig.) are compared. The FOM is also evaluated at the
same set of points in order to provide a reference for comparison. The
probability density function is approximated from the 4000 samples via the
kernel density estimation (KDE) method [34, 35] with 100 equally-spaced points
and a band-width of $0.001$ and $0.017$ for $C_{d}$ and $C_{l}$ respectively.
The density plots are shown in Figure 11. In predicting the density plot for
$C_{l}$ the two methods perform quite equally. However for $C_{d}$, the
proposed approach clearly outperforms POD-Krig. As discussed briefly in
section 7.2, this is mainly because the $C_{d}$ shows greater sensitivity to
$C_{P}$ distributions compared to $C_{l}$ particularly in 2-D inviscid
simulations and hence even modest error in predicting $C_{P}$ can translate
into larger errors in $C_{d}$. Secondly, the POD-Krig. approach does not
account for the physics of the system unlike the proposed approach and hence
performs poorly when the Kriging interpolator does not learn the behavior of
the POD coefficients in the parameter space very well. It should be noted that
there are other ways to improve upon the Kriging prediction, for instance via
using a _non-stationary_ kernel [36, 28] to capture the underlying correlation
in the data. Such models might come with a higher dimensional parametrization
but have the potential to emulate rapidly varying non-linear responses very
well. However, in the present work the goal is only to make a very general
comparison of the proposed approach with a very basic alternative method that
is applicable for non-intrusive ROM methods. The statistical measures that
quantify the shape of the density plots is summarized in Table 4. As revealed
by the table, the statistics for the proposed POD-Proj. approach matches the
stats predicted from the FOM much better than the POD-Krig approach.

Another important aspect for discussion are the computational times required
for each of the methods. All surrogate model computations were run in serial
mode on a desktop computer with 4 cores and 16GB RAM. The ROM computations for
4000 samples consumed approximately 7.2 hrs of wall-clock time whereas the
POD-Krig. approach consumed approximately 6.5 hrs. An equivalent budget of
high-fidelity simulations run in serial for the same level of convergence as
the ROM would have cost $\approx$ 667 hrs of wall-clock time; therefore the
computational speed-up with the ROM is $\sim~{}100\times$. However the FOM was
run on a high-performance computer taking advantage of parallel computing in
order to keep the wall-clock times tractable.

Overall, the projection-based ROM is able to capture general trends such as
the range of output quantities of interest and shapes of their distributions
with much better accuracy than an approach that circumvents the projection
step (in this case via Kriging). This emphasizes the power of projection-based
ROMs since they operate on the actual governing equations of the system as
opposed to other data-driven surrogate modeling techniques. This allows them
to perform more robustly with parameter variation compared to competing
methods. Overall, the presented results provide evidence that the projection-
based ROM has capability to potentially supplement the expensive high-fidelity
models while offering compelling computational speed-ups without sacrificing
the accuracy significantly.

Table 4: Comparison of output statistics. Bold-face entries are significantly worse compared to rest | $C_{l}$ | $C_{d}$
---|---|---
| POD-Proj. | POD-Krig. | FOM | POD-Proj. | POD-Krig. | FOM
Mean | 0.2996 | 0.2968 | 0.2980 | 0.0076 | 0.0334 | 0.0075
Median | 0.2989 | 0.2967 | 0.2969 | 0.0076 | 0.0340 | 0.0074
Std. Dev | 0.0587 | 0.0531 | 0.0598 | 0.0018 | 0.0068 | 0.0007
Skewness | 0.0743 | -0.0144 | 0.0431 | 0.2294 | -0.2765 | 0.3887
Kurtosis | 2.4794 | 2.4978 | 2.7697 | 4.5908 | 2.5047 | 2.5210

(a) Drag coefficient

(b) Lift coefficient

Figure 11: Kernel density estimate from 4000 monte-carlo samples.

###  5.4 Computational Costs

The overall computational cost is dominated by the offline phase where the
model is built. Here, we provide an estimate of the computational cost in
terms of Floating Point OPerationS (FLOPS) necessary to build the ROM (POD-
Proj.) as a function of the grid size ($N$) and the number of FOM snapshots
($M$). The cost of the online phase (including ROM interpolation) is trivial
comparatively and the wall-clock time is more relevant in this scenario. The
off-line phase includes 4 major steps whose computational cost are summarized
in Table 5; see [1] for details. Note that only the dominating factors of the
cost are provided in the table. In the same table, the cost of an intrusive
projection-based MOR method is also provided for comparison. It can be seen
that the most expensive steps of the method are the POD and the projection
which scale as $\sim~{}NM^{2}$ and $\sim~{}N^{2}M$ respectively; as $N$
increases the cost of these steps increases quite rapidly. The finite volume
discretization is relatively a cheaper step that scales linearly with grid
size. Compared to a conventional intrusive MOR method, the current method
incurs an additional cost due to the finite-volume discretization which can be
considered a penalty paid due to the lack of access to the source code of the
FOM. Whereas POD-Krig. is computationally cheaper than projection-based
methods, their limitation in predicting complex flow fields have been
demonstrated in this study.

Table 5: Summary of offline computational cost Operation | (non-intrusive) POD-Proj. | POD-Krig. | (intrusive) POD-Proj.
---|---|---|---
Snapshot Scaling | $\mathcal{O}(N)$ | $\mathcal{O}(N)$ | $\mathcal{O}(N)$
POD | $\mathcal{O}(NM^{2})$ | $\mathcal{O}(NM^{2})$ | $\mathcal{O}(NM^{2})$
Finite Vol. Discret. | $\mathcal{O}(MN)$ | - | -
Projection | $\mathcal{O}(N^{2}M)$ | - | $\mathcal{O}(N^{2}M)$

The wall-clock times for online evaluation are summarized in Table 6, where
the computational gains of the proposed approach are compelling and therefore
given their demonstrated accuracy in predicting non-linear flows, are viable
candidates during early and middle stages of aerospace design. While the POD-
Krig. approach requires a comparable amount of wall-clock time as the proposed
approach, they are based on prediction for only one variable ($p$) whereas the
proposed approach solves for $O=8$ observables. Therefore for a fair
comparison, the POD-Krig. wall-clock times are multiplied by $O$ in Table 6 to
provide the perspective. Having said that, when surrogates of only one field
variable (such as pressure) is of interest, this study suggests that a simpler
method such as the POD-Krig. could be sufficient.

Table 6: Comparison of computational (wall-clock) times | FOM | POD-Krig. | POD-Proj.
---|---|---|---
Evaluation | 600s | 5.85s$(\times O)$ | 6.48s
Inverse Design | 300hrs | 3.4hrs$(\times O)$ | 3.7hrs
Uncert. Quant. | 667hrs | 6.5hrs$(\times O)$ | 7.2hrs

##  6 Conclusion

We extend the earlier work by [1] to systems with parametric geometry in
addition to parametric boundary conditions, which are more commonly
encountered in the aerospace design process. The methodology is validated
under subsonic and transonic compressible inviscid flow. Under subsonic
conditions, exceptional accuracy (<5%) was observed in approximating the flow-
field as well as ouputs $C_{l}$, $C_{d}$ and $C_{P}$, whereas in the transonic
regime the average error is higher (in the 15-20% range). The method suffers
from typical limitations of a POD-based methods to capture highly non-linear
flow such as ones with moving shocks [37, 22]. A common trend observed in this
work is that $C_{P}$ and $C_{l}$ are predicted with much better accuracy than
$C_{d}$, albeit showing better accuracy and robustness than a competing non-
intrusive method.

Furthermore to demonstrate its computational efficiency, the method is applied
to 2 specific applications in the many-query context: (i) inverse design and
(ii) uncertainty quantification via Monte Carlo sampling. In both situations,
the ROM is executed $\mathcal{O}(1000)$ times at a wall-clock time of 2-8 hrs,
while the equivalent budget of FOM would have consumed 200-800 hrs. Therefore
first and foremost, these results demonstrate the utility of the methodology
towards real-time decision making. The inverse design test case is mainly the
test of the ROM to satisfy the physics of the problem, even in its
approximated form. The predictions of the ROM is within $5\%$ of the target
solution while again, achieving it at superior computational efficiency
compared to the FOM. Lastly, a Monte Carlo analysis with 4000 uniformly
sampled points from the input space was used to approximate the probability
distributions of the two main outputs considered in this work: the $C_{d}$ and
$C_{l}$. Overall, the $C_{l}$ showed better match with the FOM results, while
the $C_{d}$ predicted the statistics with relatively higher discrepancy.
Regardless, the predictions of the ROM turn out to capture general trends at a
fraction of the computational cost of the FOM, while again showing more
robustness than a competing non-intrusive method.

The present work is a first step towards performing projection-based model
reduction with black-box models. Therefore to keep the exposition simple, the
2-dimensional compressible euler equations are used as the test case. However,
the methodology naturally extends to 3-dimensional systems with $S>4$ without
modification. The only requirement is complete knowledge of the governing
equations in continuous PDE form, which is always available in the form of
code documentation with black-box commercial codes. The author is currently
investigating purely data-driven methods where this requirement can be waived.
Another direction for future work is the adaptive construction of the model,
where snapshots are sequentially generated based on certain _greedy_
criterion. Finally, state-of-the art methods to enhance the performance of the
method for flows with discontinuities is also being investigated.

Overall, the present methodology establishes feasibility for projection based
MOR for static parametric systems available as a black-box in addition to
demonstrating their utility towards aerospace design. Comparison against the
best known non-intrusive MOR technique (given the same constraints as the
present study), namely POD+Interpolation suggested that the projection-based
approach leads to more accurate predictions in the presence of parameter-
dependent discontinuities. However, the author acknowledges that projection-
based ROMs, in addition to incurring a greater off-line cost, can also lead to
numerical stability issues (specifically in time-dependent systems) which
needs to be addressed. Therefore the choice between either of these methods is
dependent upon the flow regime under consideration and hence the domain
knowledge of the engineers has to be leveraged.

## Acknowledgments

## References

  * Renganathan et al. [2018] Renganathan, S. A., Liu, Y., and Mavris, D. N., “Koopman-Based Approach to Nonintrusive Projection-Based Reduced-Order Modeling with Black-Box High-Fidelity Models,” _AIAA Journal_ , Vol. 56, No. 10, 2018, pp. 4087–4111.
  * Johnson et al. [2005] Johnson, F. T., Tinoco, E. N., and Yu, N. J., “Thirty years of development and application of CFD at Boeing Commercial Airplanes, Seattle,” _Computers and Fluids_ , Vol. 34, No. 10, 2005, pp. 1115–1151.  10.1016/j.compfluid.2004.06.005.
  * Holmes, Philip., Lumley, John L., Berkooz, Gahl and Rowley [1998] Holmes, Philip., Lumley, John L., Berkooz, Gahl and Rowley, C. W., _Turbulence, Coherent Structures, Dynamical Systems and Symmetry_ , Vol. 36, 1998.  10.2514/2.399.
  * Xiao et al. [2015] Xiao, D., Fang, F., Buchan, A. G., Pain, C. C., Navon, I. M., and Muggeridge, A., “Non-intrusive reduced order modelling of the Navier Stokes equations,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 293, 2015, pp. 522–541.  10.1016/j.cma.2015.05.015, URL http://linkinghub.elsevier.com/retrieve/pii/S0045782515001887.
  * Christophe Audouze and Nair [2013] Christophe Audouze, F. D. V., and Nair, P. B., “Nonintrusive Reduced-Order Modeling of Parametrized Time-Dependent Partial Differential Equations,” _Numerical Methods for Partial Differential Equasion_ , Vol. 29, No. 5, 2013, pp. 1587–1628.  10.1002/num, URL http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve{\&}db=PubMed{\&}dopt=Citation{\&}list{\\_}uids=20376194.
  * Audouze C. and B. [2009] Audouze C., D. V. F., and B., N. P., “Reduced-order modeling of parameterized PDEs using time-space parameter principal component analysis,” _INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING_ , Vol. 80, 2009, pp. 1025–1057.  10.1002/nme.
  * Bui-Thanh et al. [2004] Bui-Thanh, T., Damodaran, M., and Willcox, K. E., “Aerodynamic Data Reconstruction and Inverse Design Using Proper Orthogonal Decomposition,” _AIAA Journal_ , Vol. 42, No. 8, 2004, pp. 1505–1516.  10.2514/1.2159.
  * Brunton et al. [2016] Brunton, S. L., Proctor, J. L., and Kutz, J. N., “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , 2016, p. 201517384.
  * 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.
  * Koopman [1931] Koopman, B. O., “Hamiltonian systems and transformation in Hilbert space,” _Proceedings of the National Academy of Sciences_ , Vol. 17, No. 5, 1931, pp. 315–318.
  * Renganathan [2018] Renganathan, S. A., “A Methodology for Non-Intrusive projection-based model reduction of expensive black-box PDE-based systems and application in the many-query context,” Ph.D. thesis, Georgia Institute of Technology, 2018.
  * Kutz et al. [2016] Kutz, J. N., Proctor, J. L., and Brunton, S. L., “Generalizing Koopman theory to allow for inputs and control,” _arXiv:1602.07647v1_ , 2016, pp. 1–21.
  * Chaturantabut and Sorensen [2010] Chaturantabut, S., and Sorensen, D. C., “Nonlinear Model Reduction via Discrete Empirical Interpolation,” _SIAM Journal on Scientific Computing_ , Vol. 32, No. 5, 2010, pp. 2737–2764.  10.1137/090766498.
  * Schittkowski [1986] Schittkowski, K., “NLPQL: A FORTRAN subroutine solving constrained nonlinear programming problems,” _Annals of operations research_ , Vol. 5, No. 2, 1986, pp. 485–500.
  * Golub and Van Loan [2012] Golub, G. H., and Van Loan, C. F., _Matrix computations_ , Vol. 3, JHU Press, 2012.
  * Rahman et al. [2005] Rahman, I. U., Drori, I., Stodden, V. C., Donoho, D. L., and Schröder, P., “Multiscale Representations for Manifold-Valued Data,” _Multiscale Modeling & Simulation_, Vol. 4, No. 4, 2005, pp. 1201–1232.  10.1137/050622729.
  * Barachant et al. [2010] Barachant, A., Bonnet, S., Congedo, M., and Jutten, C., “Riemannian geometry applied to BCI classification,” _International Conference on Latent Variable Analysis and Signal Separation_ , Springer, 2010, pp. 629–636.
  * STA [2017] “STARCCM+ url: http://mdx.plm.automation.siemens.com/star-ccm-plus,” , 2017\.  URL http://mdx.plm.automation.siemens.com/star-ccm-plus.
  * Poirier et al. [1998] Poirier, D., Allmaras, S. R., McCarthy, D. R., Smith, M. F., and Enomoto, F. Y., “The CGNS system,” _AIAA paper_ , , No. 98-3007, 1998.
  * Kulfan and Bussoletti [2006] Kulfan, B., and Bussoletti, J., “Fundamental Parameteric Geometry Representations for Aircraft Component Shapes,” _11th AIAA/ISSMO multidisciplinary analysis and optimization conference_ , 2006, p. 6948.
  * Kulfan [2008] Kulfan, B. M., “Universal parametric geometry representation method,” _Journal of Aircraft_ , Vol. 45, No. 1, 2008, pp. 142–158.
  * LeGresley [2006] LeGresley, P. A., “Application of proper orthogonal decomposition (POD) to design decomposition methods,” Ph.D. thesis, 2006.
  * Lucia et al. [2002] Lucia, D. J., King, P. I., and Beran, P. S., “Domain Decomposition for Rediced-Order Modeling of a Flow with Moving Shocks,” _AIAA Journal_ , Vol. 40, No. 11, 2002, pp. 2360–2363.
  * Welper [2017] Welper, G., “Interpolation of functions with parameter dependent jumps by transformed snapshots,” _SIAM Journal on Scientific Computing_ , Vol. 39, No. 4, 2017, pp. A1225–A1250.
  * Goldberg [1989] Goldberg, D. E., _Genetic Algorithms in Search, Optimization and Machine Learning_ , 1st ed., Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1989.
  * D. Xiao and Hu [2015] D. Xiao, C. P., F. Fang, and Hu, G., “Non-intrusive reduced-order modelling of the Navier–Stokes equations based on RBF interpolation,” _INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS_ , Vol. 79, 2015, pp. 580–595.
  * Stein [2012] Stein, M. L., _Interpolation of spatial data: some theory for kriging_ , Springer Science & Business Media, 2012.
  * Rasmussen [2003] Rasmussen, C. E., “Gaussian processes in machine learning,” _Summer School on Machine Learning_ , Springer, 2003, pp. 63–71.
  * Santner et al. [2003] Santner, T. J., Williams, B. J., Notz, W., and Williams, B. J., _The design and analysis of computer experiments_ , Vol. 1, Springer, 2003.
  * Anderson Jr [2010] Anderson Jr, J. D., _Fundamentals of aerodynamics_ , Tata McGraw-Hill Education, 2010.
  * Adkins and Liebeck [1994] Adkins, C. N., and Liebeck, R. H., “Design of optimum propellers,” _Journal of Propulsion and Power_ , Vol. 10, No. 5, 1994, pp. 676–682.
  * [32] Renganathan, A., Denney, R. K., Duquerrois, A., and Mavris, D. N., _Validation and Assesment of Lower Order Aerodynamics Based Design of Ram Air Turbines_ , ????  10.2514/6.2014-3463, URL https://arc.aiaa.org/doi/abs/10.2514/6.2014-3463.
  * Mooney [1997] Mooney, C. Z., _Monte carlo simulation_ , Vol. 116, Sage Publications, 1997\.
  * Epanechnikov [1969] Epanechnikov, V. A., “Non-parametric estimation of a multivariate probability density,” _Theory of Probability & Its Applications_, Vol. 14, No. 1, 1969, pp. 153–158.
  * Rosenblatt [1956] Rosenblatt, M., “Remarks on some nonparametric estimates of a density function,” _The Annals of Mathematical Statistics_ , 1956, pp. 832–837.
  * Paciorek and Schervish [2004] Paciorek, C. J., and Schervish, M. J., “Nonstationary covariance functions for Gaussian process regression,” _Advances in neural information processing systems_ , 2004, pp. 273–280.
  * Beran and Silva [2001] Beran, P. S., and Silva, W., “Reduced Order Modeling: New Approaches for Computational Physics,” _42rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference_ , 2001.

##  7 Appendix

###  7.1 Shape parameterization

The CST model of parameterization defines a class function $C$ and a shape
function $S$ and the curve being parameterized is specified as their product.
The class function, $C(\psi)$ is more generally defined as

$C_{n_{1}}^{n_{2}}(\psi):=\psi^{n_{1}}(1-\psi)^{n_{2}}$

where the variable $\psi$ represents the non-dimensional chord-wise distance,
$n_{1}$ and $n_{2}$ define the specific class. The unique shape of an airfoil
is driven by the shape function, specified as follows

$S(\psi)=\sum_{i=0}^{n}A_{i}\psi^{i}$ (30)

It is particularly useful to define a unit shape function, i.e. $S(\psi)=1$
such that the individual coefficients $A_{i}$ 333 $A_{i}$’s are denoted by
$K_{i,n}$ for unit shape functions can be obtained as generic constants. For
instance for $n=1$ the simplest decomposition one could get for the shape
function is $S(\psi)=S_{0}(\psi)+S_{1}(\psi)$ where $S_{1}(\psi)=\psi$ and
$S_{2}(\psi)=1-\psi$ where the coefficients $A_{0}=1$ and $A_{1}=1$.
Similarly, for the general $nth$ order shape function, the decomposition of
the unit shape function can be done using Bernstein polynomials

$S(\psi)=\sum_{i=0}^{n}K_{i,n}\psi^{i}(1-\psi)^{n-i}$ (31)

where the coefficients are the binomial coefficients given by

$K_{i,n}={n\choose i}=\frac{n!}{i!~{}(n-i)!}$

The final shape of the airfoil shape is then given by

$\mathbf{y}(\psi)=C(\psi)S(\psi)$ (32)

The unit shape functions and the corresponding airfoil geometries are
illustrated in the Figure 12. It can be seen that such a parametrization
results in each component shape function peak being equally distributed
between the leading and trailing edges leading to the same effect in the
component airfoils. It is now a matter of scaling up or down, the binomial
coefficients of the Bernstein polynomials in order to approximate the unique
airfoil shape of interest.

Figure 12: Examples of the shape function decomposition (into Bernstein
polynomials) for various values of the order $n$ (left) and the resulting
component airfoils (right). The coefficients $A_{i}$ correspond to unit shape
function, which can be scaled up/down to obtain a specific airfoil shape.

The coefficients $A_{i}$ represent the actual parameters of the shape, given
$n$ the order of the Bernstein polynomials. An $nth$ order CST parametrization
has $n+1$ parameters. If separate parametrizations are sought for the upper
and lower surfaces of the airfoil, then the CST parametrization leads to
$2(n+1)$ parameters to specify the whole shape of the airfoil, where the $n$
needs to be determined for a specific geometry under consideration. However,
typically $n=3-5$ are observed to be adequate to parametrize the airfoil
shapes considered in this work. One way to determine $n$ and the associated
polynomial coefficients is to find the values that minimizes certain error
between the true shape of the airfoil and the resulting approximation via CST.
In this work, the parameters for a given airfoil shape are determined by
solving the following minimization problem

$\underbrace{\text{minimize}}_{A_{j},~{}n}~{}\left\|\mathbf{y}(\tilde{\psi})-C(\tilde{\psi})S(\tilde{\psi},~{}A_{j})\right\|_{2}^{2}$
(33)

where $\tilde{\psi}\in\mathbb{R}^{n+1}$ are $n+1$ equally spaced points
sampled from $\psi$ spanning $[0,1]$ 444Note that picking equally spaced
points works well for the current airfoil geometries under consideration in
this work and may not work for any arbitrary geometry. More generically, a
least-squares fit considering all the points might be more suitable. This way,
the smallest possible $n$ and their corresponding Bernstein coefficients are
determined. For the RAE2822 airfoil shape, the following parameterization was
obtained ($n=3$):

$A_{RAE2822}=\begin{bmatrix}0.1268&0.4670&0.5834&0.2103\\\
-0.1268&-0.5425&-0.5096&0.0581\end{bmatrix}$

where the first and second rows represent the parameterization of the upper &
lower surfaces of the airfoil; the comparison of the CST curve and the actual
RAE2822 shape is shown in Figure 4(a). It can be seen that the CST
parametrization gives an adequate approximation to the true curve with only 8
parameters. The coefficients may now be perturbed to modify the baseline
airfoil shape.

Similarly, the NACA0012 airfoil shape is approximated via CST and is also
shown in Figure 4(b). In this case, due to the lack of camber, the CST gives
very good approximation with $n=2$. Additionally, due to the symmetry of the
airfoil about the chord, the parameters (given below) are equal in magnitude
and opposite in sign. However, all the $2(n+1)$ degrees of freedom are
considered in this work for the NACA0012 airfoil.

$A_{NACA0012}=\begin{bmatrix}0.1689&0.2699&0.1387\\\
-0.1689&-0.2699&-0.1387\end{bmatrix}$

###  7.2 Model validation data

The values of the output quantities of interest and their associated errors
(defined by Eq. 19) is presented in Tables 7 and 8 below.

Table 7: Comparison of $C_{P}$, $C_{l}$ and $C_{d}$ between ROM & FOM for the NACA0012 test case Case |  $C_{P}$ Error % |  $C_{l}$ (ROM) |  $C_{l}$ (FOM) | Error %
---|---|---|---|---
1 | 1.29 | 0.1889 | 0.1912 | 1.20
2 | 0.74 | 0.2018 | 0.2070 | 2.50
3 | 0.80 | 0.2932 | 0.2943 | 0.37
4 | 1.86 | 0.2795 | 0.2865 | 2.44
5 | 1.36 | 0.3550 | 0.3621 | 1.96
6 | 0.62 | 0.3691 | 0.3664 | 0.73
7 | 0.46 | 0.3298 | 0.3272 | 0.79
8 | 2.79 | 0.3229 | 0.3312 | 2.50
9 | 0.76 | 0.3109 | 0.3137 | 0.89
10 | 5.43 | 0.2710 | 0.3065 | 11.58
Table 8: Comparison of $C_{P}$, $C_{l}$ and $C_{d}$ between ROM & FOM for the RAE2822 test case Case |  $C_{P}$ Error % |  $C_{d}$ (ROM) |  $C_{d}$ (FOM) | Error % |  $C_{l}$ (ROM) |  $C_{l}$ (FOM) | Error %
---|---|---|---|---|---|---|---
1 | 7.68 | 0.0161 | 0.0174 | 7.47 | 0.9174 | 0.9608 | 4.52
2 | 12.93 | 0.0336 | 0.0302 | 11.26 | 1.0669 | 0.9825 | 8.59
3 | 2.84 | 0.0257 | 0.0262 | 1.91 | 1.1140 | 1.1446 | 2.67
4 | 8.14 | 0.0321 | 0.0264 | 21.59 | 1.1004 | 1.0901 | 0.95
5 | 13.77 | 0.0194 | 0.0224 | 13.39 | 0.8965 | 1.0345 | 13.34
6 | 12.78 | 0.0298 | 0.0484 | 38.43 | 1.0288 | 1.0082 | 2.04
7 | 5.91 | 0.0306 | 0.0279 | 9.68 | 0.9403 | 0.9286 | 1.26
8 | 4.99 | 0.0245 | 0.0245 | 0 | 0.9124 | 0.8914 | 2.36
9 | 5.63 | 0.0292 | 0.0315 | 7.30 | 0.9507 | 0.9779 | 2.78
10 | 9.96 | 0.0326 | 0.0217 | 50.23 | 0.9750 | 0.9793 | 0.44

###  7.3 Discrete Empirical Interpolation Method (DEIM)

The Discrete Empirical Interpolation Method (DEIM) is briefly reviewed here
and as an illustration one of the non-linear constraints used in Eq. 18 is
evaluated. For a non-linear function $\mathbf{f}(\theta)\in\mathbb{R}^{N}$ the
DEIM approximates $\mathbf{f}$ by projecting it onto a subspace spanned by
$\\{\mathbf{x}_{1},...,\mathbf{x}_{q}\\}\subset\mathbb{R}^{N}$ as

$\mathbf{f}(\theta)\approx\mathbf{X}c(\theta)$ (34)

where $\mathbf{X}=[\mathbf{x_{1}},...,\mathbf{x_{q}}]\in\mathbb{R}^{N\times
q},~{}~{}q<<N$ is determined via a POD of the snapshots of $\mathbf{f}$ and is
assumed to be globally valid in the design space that bounds the design
parameters $\theta$ and $\mathbf{c}(\theta)\in\mathbb{R}^{q}$ are the
coefficients of the basis expansion. Then the approximation of $\mathbf{f}$
requires only the determination of $\mathbf{c}(\theta)$ which requires only
$q$ equations. The DEIM gives a distinguished set of $q$ points from the over-
determined system $\mathbf{f}(\theta)=\mathbf{X}\mathbf{c}(\theta)$. Given a
permutation matrix $\mathbf{P}$ that would give $q$ such distinguished rows of
a matrix when pre-multiplied, then the $q\times q$ system necessary to solve
for the coefficients is given by

$\mathbf{P}^{\top}\mathbf{f}(\theta)=(\mathbf{P}^{\top}\mathbf{X})\mathbf{c}(\theta)$
(35)

So the approximation of $\mathbf{f}(\theta)$ is then given by

$\mathbf{f}(\theta)\approx\mathbf{X}(\mathbf{P}^{\top}\mathbf{X})^{-1}\mathbf{P}^{\top}\mathbf{f}(\theta)$
(36)

If the $q$ row-indices (that are extracted by pre-multiplying with
$\mathbf{P}^{\top}$) are represented by a vector, $\mathbf{\varrho}$, then in
the above equation, $\mathbf{P}^{\top}\mathbf{f}(\theta)$ is equivalent to
extracting the $\varrho$ rows of $\mathbf{f}$. Therefore the approximation of
$\mathbf{f}(\theta)$ requires only $q$ computations which is efficient because
$q<<N$. Similarly, a non-linear function that depends on the state,
$\mathbf{f}(\mathbf{u})$ can be approximated as

$\mathbf{f}(\mathbf{u})\approx\mathbf{X}(\mathbf{P}^{\top}\mathbf{X})^{-1}\mathbf{P}^{\top}\mathbf{f}(\mathbf{u})$
(37)

Since $\mathbf{u}=\Phi_{k}^{\top}\tilde{\mathbf{u}}$ and setting
$\tilde{\mathbf{f}}=\Phi_{k}^{\top}\mathbf{f}(\mathbf{u})$,
$\tilde{\mathbf{f}}$ can be approximated as

$\tilde{\mathbf{f}}=\Phi_{k}^{\top}\mathbf{X}(\mathbf{P^{\top}}\mathbf{X})^{-1}\mathbf{f}(\mathbf{P}^{\top}\Phi_{k}\tilde{\mathbf{u}})$
(38)

In the above equation, the term
$\Phi_{k}^{\top}\mathbf{X}(\mathbf{P^{\top}}\mathbf{X})^{-1}$ is independent
of the state and hence can be pre-computed and $\mathbf{P}^{\top}\Phi_{k}$ is
just extraction of the $\varrho$ rows of $\Phi_{k}$. Therefore using the DEIM,
the non-linear term can be expressed in terms of the reduced state,
$\tilde{\mathbf{u}}$ and hence can be efficiently computed.

Now the DEIM is illustrated on evaluating the first constraint of Equation 18
which in discretized form is given below

$\mathbf{h}_{1}=\mathbf{y}_{5}-\frac{\mathbf{y}_{1}\mathbf{y}_{3}}{\mathbf{y}_{2}}$
(39)

Let $\varrho_{5}$ be the vector containing the $q$ row-indices returned by
DEIM via snapshots of the non-linear term $\mathbf{y}_{5}$ and $\Phi_{1}$,
$\Phi_{2}$, $\Phi_{3}$, $\Phi_{5}$ be the projection matrix of
$\mathbf{y}_{1}$, $\mathbf{y}_{2}$, $\mathbf{y}_{3}$ and $\mathbf{y}_{5}$
respectively. Then

$\tilde{\mathbf{h}}_{1}=\tilde{\mathbf{y}_{5}}-\Phi_{5}^{\top}\mathbf{X}~{}[\mathbf{X}(\varrho_{5},:)]^{-1}\left\\{\frac{\Phi_{1}(\varrho_{5},:)\tilde{\mathbf{y}_{1}}~{}~{}\Phi_{3}(\varrho_{5},:)\tilde{\mathbf{y}_{3}}}{\Phi_{2}(\varrho_{5},:)\tilde{\mathbf{y}_{2}}}\right\\}$
(40)

In the above equation, the term outside of the braces can be pre-computed.
Additionally since
$\mathbf{y}_{5}=\frac{\mathbf{y}_{1}\mathbf{y}_{3}}{\mathbf{y}_{2}}$,
$\mathbf{X}=\Phi_{5}$ and hence the term reduces to
$[\mathbf{X}(\varrho_{5},:)]^{-1}$ which is $q\times q$ and hence can be
cheaply computed. Therefore using the DEIM, the non-linear constraints are
evaluated in terms of the reduced state variables which makes it
computationally cheap.