Transported snapshot model order reduction approach for parametric, steady-state fluid flows containing parameter dependent shocks
Nirmal J. Nair Maciej Balajewicz
physics.flu-dyn




1]Nirmal J. Nair*

1]Maciej Balajewicz

NAIR AND BALAJEWICZ


Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA

*Nirmal J. Nair, Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA. 

 njn2@illinois.edu

Air Force Office of Scientific Research, Grant/Award Number: FA9550-17-1-0203

[Summary]A new model order reduction approach is proposed for parametric steady-state
nonlinear fluid flows characterized by shocks and discontinuities whose spatial
locations and orientations are strongly parameter dependent. In this method,
solutions in the predictive regime are approximated using a linear
superposition of parameter dependent basis. The sought after parametric
reduced-basis are obtained by transporting the snapshots in a spatially and parametrically dependent transport field. Key to the
proposed approach is the observation that the transport fields are
typically smooth and continuous, despite the solution themselves not being so.
As a result, the transport fields can be accurately expressed using a
low-order polynomial expansion. Similar to traditional projection-based model order reduction
approaches, the proposed method is formulated mathematically as a residual
minimization problem for the generalized coordinates. The proposed approach is also integrated with well-known hyper-reduction strategies to obtain significant computational speed-ups.
The method is successfully applied to the reduction of a parametric 1-D flow in a
converging-diverging nozzle, a parametric 2-D supersonic flow over a forward
facing step and a parametric 2-D jet diffusion flame in a combustor.




,
 
Transported snapshot model order reduction approach for parametric, steady-state fluid flows containing parameter dependent shocks. Int J Numer Methods Eng. 2018;1-29.
<https://doi.org/10.1002/nme.5998>

Transported snapshot model order reduction approach for parametric,
steady-state fluid flows containing parameter dependent shocks
    Balajewicz M.
    December 30, 2023  File: 
==================================================================================================================================




 


§ INTRODUCTION


Computational models of high-dimensional systems arise in a rich variety of
engineering and scientific contexts. Computational Fluid Dynamics (CFD) for
example has become an indispensable tool for many engineering applications
across a wide range of industries.  Unfortunately, high-fidelity CFD simulations
are often so computationally prohibitive that they cannot be used as often as
needed or used only in special circumstances rather than routinely.
Consequently, the impact of CFD on parametric and time-critical applications
such as design, optimization, and control has not yet been fully
realized. Model Order Reduction (MOR) is a serious contender for bridging this
gap. 

Most existing MOR approaches are based on projection. In projection-based MOR,
the state variables are approximated in a low-dimensional subspace. Over the
years, a number of approaches for calculating a reduced-order basis (ROB) have
been developed, e.g., proper orthogonal decomposition
(POD) <cit.>, dynamic mode decomposition
(DMD) <cit.>, balanced POD
(BPOD) <cit.>, balanced
truncation <cit.> and the reduced-basis
method <cit.>. 

For a reduced-order model (ROM) to be truly useful, it must be capable of generating accurate
predictions for parameter values that are different from those sampled for the
purpose of constructing a ROB. Generating ROBs and ROMs that are robust to
parameter variations is an active area of research <cit.>. The
choice of parameter sample points is critical to any method used to generate
the basis. For problems with smaller number of parameters, a simple approach
like random sampling using the Latin hypercube method is often
sufficient <cit.>.  For problems with large number of
parameters, more sophisticated sampling methods are usually required.  In the
standard greedy sampling
approach <cit.>,
the sample points are chosen one-by-one in an adaptive manner. At every
iteration, the goal is to find the parameter value for which the error between
the ROM and the full order model (FOM) is largest. The FOM is sampled at this point and the new
information is used to generate a new reduced-basis.

Another critical issue involves the choice between a global or local basis.
Although global basis have been shown to perform adequately for many
applications <cit.>, particularly challenging problems often necessitate the
use of multiple local reduced
basis <cit.>. In these
cases, several local basis are constructed and linked to particular regions in
the parameter or state
space <cit.>.
The price of this additional flexibility is the switching algorithm required
to switch between the local basis. Finally, to improve the ROM performance in
predictive regimes, it is also possible to interpolate the basis or the ROM
matrices <cit.>.

Achieving parametric robustness is particularly challenging when the sought
after solutions contain sharp gradients, discontinuities or shocks. These
situations arise in a wide range of important engineering applications,
for example, high-speed fluid
flows <cit.>, multi-phase flows with
evolving material interfaces <cit.>, computational
finance <cit.> and structural contact problem with
evolving contact regions <cit.>. 

Over the years, a large variety of discontinuity-aware MOR techniques have been
developed. In the first class of such methods, the problem of modeling
discontinuities is avoided entirely by exploiting symmetry and transport
reversal properties of certain hyperbolic
PDEs <cit.>.  Other
methods involve decomposition into global and advection modes governed by
optimal mass transfer <cit.> or more direct modeling of
discontinuities using basis splitting <cit.> or snapshot
transformation <cit.>. For unsteady solutions with
shocks, accurate low-rank solutions can be obtained using a Lagrangian framework
 <cit.> where both position and state of Lagrangian
particle variables are approximated by their respective ROBs.  Other methods
decompose the solution into a variable separable form consisting of an evolution
term to capture moving shocks and a diffusion term to capture the changing
shapes <cit.>. Finally, other methods avoid the problem of
modeling discontinuities by domain decomposition where direct numerical
simulation or reconstruction methods are used in regions containing the
discontinuities
 <cit.>. 

In this manuscript we outline our new proposed approach for parametric model
reduction of solutions containing moving shocks and discontinuities. In our
proposed approach, solutions in the predictive regime are approximated using parameter dependent basis. Snapshots are transported in a spatially and parametrically dependent transport field to yield transported snapshots. These transported snapshots form the parameter dependent basis for the proposed MOR. Key to our proposed approach is the
observation that the transport fields are smooth and thus, can be
themselves approximated using a low-order expansion. Our method may be
interpreted as a data-driven generalization of previous works
of Rowley and Marsden <cit.> and Iollo and
Lombardi <cit.>.

The remainder of the paper is organized as follows.  In  <ref>, the problem of interest and the traditional projection-based model
order reduction approach is recapitulated.  The proposed methodology for
approximating unsampled solutions via transported snapshots for one-dimensional parameter variational problems is detailed  in
 <ref> and its extension to multi-dimensions is detailed in  <ref>. The details for numerical implementation of our new approach in the discrete framework are provided in  <ref>.
In  <ref>, the proposed approach is integrated with a hyper-reduction algorithm. In  <ref>, the novelty of the proposed approach is compared with similar existing methods in literature. In
 <ref>, the performance of the proposed method is
evaluated on several simple but representative fluid flow models. Limitations of the proposed approach are discussed in  <ref>. Finally, in
 <ref>, conclusions are offered and prospects for future work
are summarized.





§ PROBLEM STATEMENT





 §.§ Full order model


We consider full order models consisting of
hyperbolic or convection-dominated parabolic PDEs such as the Euler or
high-Reynolds-number Navier-Stokes equations:


    ∂u/∂ t+f(u,x;μ)=0


where the state variable u=u(x,t;μ) depends on space x∈Ω, Ω being the flow domain, time t ∈ [0,t_max] and a vector of N_d parameters μ∈𝒟⊂ℝ^N_d
(where 𝒟 is a bounded domain). f is the nonlinear function that  contains the convective and diffusive fluxes.

The steady state equation for this system can
be written as:


    R(u;μ) f(u,x;μ)=0


where R(u;μ) is the steady state residual.  To obtain steady state
solutions, Eq. (<ref>) can be discretized in space by a standard finite difference/volume/element method. The resulting set of equations can be solved directly by an iterative method or a
time-stepping method can be used to advance the semi-discretized form of unsteady Eq.
(<ref>) to a steady state.





 §.§ Traditional nonlinear model reduction


In traditional projection-based MOR, the state variable u(x;μ) is approximated in a global low-dimensional trial subspace as follows


    u(x;μ) 
        ≈u_r = ∑_n=1^N_ku_n(x) a_n(μ)


where u_n is the basis of this subspace, N_k is the number of basis, and
a(μ) ∈ℝ^N_k denotes the generalized coordinates in
this basis. Substituting the approximation (<ref>) into the residual equation (<ref>) yields
R(u_r;μ)=0.
Consequently, the generalized coordinates are chosen to minimize the
residual of the Galerkin expansion.


    (v_n, R(u_r;μ))_Ω


where v_n is the basis of the test subspace and the inner product is defined as


    (u,v)_Ω:=∫_Ωu.v dx


If the test basis v_n ≠u_n, then the projection in Eq. (<ref>) is called as Petrov-Galerkin projection with the specific case v_n =Ju_n known as least-squares Petrov-Galerkin (LSPG) projection, where J=∂R/∂a denotes the Jacobian of the residual. For nonlinear, non-self adjoining problems such as those represented in
this case by the set of ODEs, this approach is more robust than a Galerkin
projection, where v_n = u_n. 

Solving the minimization problem in equation (<ref>) requires the evaluation of the residual of the
governing equations of the high-dimensional state variables. The complexity of
this computation scales with the size of the FOM. Therefore, while MOR approximates the FOM in a low-dimensional subspace, part of the computational cost still scales with the size of the FOM.
For general nonlinear systems, an additional level of approximation – sometimes
called “hyper-reduction” such as the discrete empirical interpolation method (DEIM) <cit.>, Gauss-Newton with approximated tensors (GNAT) method <cit.>, energy-conserving sampling and weighting (ECSW) method <cit.>
– is therefore required.




 §.§ Drawbacks



In data-driven projection-based MOR, the reduced-order bases are usually
constructed offline by collecting solution snapshots u(x;μ_s) of  
problem (<ref>) for different instances μ_s for s=1,⋯,N_s, of the parameter
vector μ.  
The reduced-order bases are then
formed by selecting a small subset of the snapshots or via compression using,
for example, POD.

For a ROM to be useful however, it must be capable of providing solutions at
parameters μ^* not sampled during the offline basis construction
phase, μ^* ≠μ_s. Although parameter robustness is an active
area of research, it is particularly challenging when the sought after solutions
contain discontinuities or sharp gradients whose spatial orientations are
strongly parameter dependent. 

To illustrate, consider the simple problem of a quasi 1-D supersonic flow in a
converging-diverging nozzle governed by 1-D Euler equations. Area profile, A(x), of the
nozzle is parabolic with equal inlet and outlet area, A(0)=A(L), and the
throat is located at L/2, where L is the length of the nozzle. For this
problem, the throat area μ=A(L/2) is the varying parameter of interest. Refer to
 <ref> for details of this problem. Four snapshots at parameters μ_s=[0.5, 0.875, 1.25, 1.625]
are generated and the corresponding steady density solutions are shown in
Fig. <ref>.




Optimal construction of a new solution at an unsampled parameter μ^* is
given by a superposition of the sampled snapshots:


    u_r(x;μ^*)=∑_s=1^N_su(x;μ_s)a_n(μ^*)


where the coordinates, a(μ^*), are
obtained by projecting the FOM solution onto the snapshots. The FOM and optimally constructed solutions at μ^*=1.0 (μ_2<μ^*<μ_3) are plotted in Fig. <ref>. It can be observed that the optimally constructed solution is dominated by staircase shock type error. Moreover, the optimal coordinates are a(μ^*)=[-0.006, 0.333, 0.688, -0.016] where the coordinates corresponding to the first and fourth snapshots, (a_1 and a_4, respectively), are significantly lower than the coordinates corresponding to the second and third snapshots, (a_2 and a_3, respectively). This implies that the optimal construction is significantly dominated by local nearby snapshots.

This toy problem demonstrates that: (a) optimal constructions
of solutions characterized by parameter dependent shocks and discontinuities are
typically local and sparse in the sense that only two bases were required for the construction, and (b) the optimal construction provides a “staircase” approximation of the true solution. In summary, neither global nor local reduced-basis can be expected to yield
efficient approximations in the predictive regime of solutions characterized by
shocks, discontinuities, and sharp gradients whose physical locations and
orientations are parameter dependent.




§ TRANSPORTED SNAPSHOT MODEL ORDER REDUCTION (TSMOR)



In this manuscript, we introduce and summarize our new MOR approach for parametric
and steady nonlinear fluid flows characterized by moving shocks, discontinuities
and sharp gradients. In this section, the proposed method is initially developed and detailed for one-dimensional parameter variational problems, while the extension of this method for multi-dimensional parameter variations are later detailed in  <ref>.

Our proposed approach is motivated by the observation that for many problems of
interest, such as the motivational problem illustrated in Fig. <ref>,
snapshots local in the parameter space are transported in physical space. More
precisely, if u(x;μ_j) and u(x;μ_j+k) are the
sampled solutions of the FOM at parameters μ_j and μ_j+k,
respectively, then there exists continuous and smooth transport field
c_j(x,Δμ) such that


    u(x;μ_j+k) ≈u(x + c_j(x,Δμ);μ_j)


where Δμ = μ_j+k - μ_j is the parameter variation
between the original snapshot u(x,μ_j) and the target snapshot
u(x,μ_j+k). 

Approximating the solution at an unsampled parameter  μ^* ≠μ_s,
for s = 1,…,N_s proceeds similarly to the traditional projection-based
MOR. More specifically, we assume that the solution can be approximated as a
linear superposition of parameter dependent basis functions


    u(x;μ^*) ≈u_r = ∑_n=1^N_ku_n(x;μ^*) a_n(μ^*)


where the reduced-basis functions u_n correspond to
the transported local snapshots


    u_n = u(x + c_k_n(x,Δμ);μ_k_n)   for   n=1,…,N_k


where the transport field c_k_n(x;Δμ) is evaluated at the unsampled parameter variation, Δμ = μ^* - μ_k_n and the snapshots, u(x,μ_k_n) for
n=1,⋯,N_k, are a subset of the solution snapshots computed offline, μ_k_n∈μ_s.

Finally, the generalized coordinates, a(μ^*), are
identified by minimizing the residual of the Petrov-Galerkin expansion. 


    (v_n, R(u_r;μ^*))_Ω


where (·,·)_Ω is the standard inner product as defined in Eq. (<ref>).

In summary, a new solution at an unsampled parameter using our proposed approach is given by a linear superposition of transported snapshots. This approach can be decomposed into the standard 
offline-online strategy where the transport fields are identified offline while the residual of the Petrov-Galerkin projection is minimized in the online stage.


Remark 1. The transport fields in Eq. (<ref>) are identified offline via a training procedure. Clearly, any direct
identification procedure for these fields would be intractable for large-scale
systems and, very likely, yield an over-determined and ill-conditioned system.
These computational issues can be avoided by adding smoothness constraints to the transport fields. Details of this procedure are
outlined in  <ref>.

Remark 2. For the sake of brevity, we have summarized the approach for 
one-dimensional parameter variations. However, the extension of this method to
multi-dimensional parameter variations are straightforward and the details of
this extension are provided in  <ref>.

Remark 3. In practice, FOMs are typically derived in a semi-discrete form by discretizing a system of PDEs in space. Hence, usually, the solution
snapshots are available only as vectors and not as functions.
Consequently, it is not possible to simply evaluate the transported snapshots
for arbitrary transport field c_k_n(x,Δμ) as shown in Eq. (<ref>). For discrete
models, this step must be performed using a numerical interpolation procedure.
Computational details of this procedure are outlined in  <ref>. 

Remark 4. Since the reduced-bases are parameter dependent, for
every new prediction at an unsampled parameter, an entire set of new reduced-bases must be
generated. However, the computational costs of generating these bases can be
expected to be proportional to the size of the FOM. Furthermore, the computational cost associated with the computation of the nonlinear residual function also scales with the size of the FOM. This expensive cost of residual evaluation can be mitigated by precomputation of certain terms that contain polynomial nonlinearities. However, such a precomputation procedure may not be viable for other classes of nonlinear functions. Hence, to gain significant amounts of computational speed-up,
a hyper-reduction strategy must be utilized. Details of this procedure are
outlined in  <ref>.





 §.§ Offline stage



Given the set of snapshots at sampled parameters μ_s, the methodology for identifying the transport fields c_s(x;Δμ), for
s=1,⋯,N_s, where Δμ=μ-μ_s, for each snapshot is
explained in this section.

As mentioned in Remark 1 of  <ref>, to make the identification procedure for the transport field tractable for large systems, smoothness constraints are required. More precisely, since the transport fields are assumed to be smooth in space, they can be approximated using a low order polynomial expansion in space using spatial basis functions f_p(x). The basis functions f_p(x) are selected a priori such as Chebyshev polynomials or Fourier modes. Then, similar to polynomial fitting, the transport field is expressed as a polynomial expansion in parameter space as:


    c_s(x;Δμ)=  ∑_p=1^N_p∑_q=1^N_q c_m^sf_p(x)g_q(Δμ)   for   m=1,…,N_p N_q


where c_m^s are the coefficients of the expansion, N_p and N_q are the number of functions f_p(x) and g_q(Δμ), respectively, and g_q(Δμ) are the monomials given by:


    g_q(Δμ)= (Δμ)^q  for   q=1,…,N_q


where Δμ=μ-μ_s.
The choice of the basis functions f_p(x) is problem-dependent, hence it is detailed in  <ref> with the help of test problems.

The least-squares fitting procedure to identify the coefficients c^s:= c_m^s is now explained. Firstly, a subset of N_v solution snapshots u(x;μ_v_n) for n=1,⋯,N_v, where μ_v_n∈μ_s, is chosen. Here, the parameters μ_v_n are called as the training parameters for the parameter μ_s and the corresponding snapshots are called training snapshots. The procedure for selecting N_v training snapshots from N_s is described in  <ref>. The coefficients c^s are evaluated by minimizing the least-squares training error between the transported snapshot u(x+c_s(x;Δμ);μ_s) and training snapshots u(x;μ_v_n):


    c^smin∑_n=1^N_v( ϵ_n, ϵ_n )_Ω


where ϵ_n(x)=u(x+c_s(x;Δμ);μ_s)-u(x;μ_v_n); Δμ=μ_v_n-μ_s and (·,·)_Ω is the standard inner product as defined in Eq. (<ref>).
The optimization problem (<ref>) is solved N_s times for the coefficients
c^s for s=1,⋯,N_s snapshots. 

In summary, the identification of the transport field in Eq. (<ref>) is posed as a least-squares fitting problem for the coefficients of the polynomial fitting function.




  §.§.§ Choice of training snapshots



According to Eq. (<ref>), a snapshot at μ_j+k can be approximated by transporting its neighboring snapshot at μ_j. The accuracy of this approximation tends to decrease as the absolute value, |k|, increases. Hence, N_v training snapshots are chosen which correspond to the N_v neighboring snapshots to u(x;μ_s) in the parameter space.

For instance, consider the sampled parameters μ_1<…<μ_5 in Fig. <ref>. For the evaluation of the transport field for the snapshot at μ_3, we have μ_s=μ_3. The snapshot at μ_3 has N_v=2 nearest neighboring snapshots at μ_2 and μ_4 which form the training snapshots μ_v_1 and μ_v_2, respectively. It is noted that the offline training procedure does not restrict itself to N_v=2 and allows for more neighboring snapshots to be incorporated. Finally, a biased stencil can be used at the boundary of the parameter space for the selection of the training snapshots.








 §.§ Online stage

In this stage, parameter dependent reduced-bases u_n(x;μ^*) at an unsampled parameter μ^* (≠μ_s) are constructed and generalized coordinates a(μ^*) are identified.

For a parameter μ^* in the predictive regime, local reduced-bases u_n are evaluated by transporting local snapshots u(x;μ_k_n) for
n=1,⋯,N_k, as shown in Eq. (<ref>), where μ_k_n∈μ_s is a subset of snapshot solutions. The procedure for selecting the subset of snapshots u(x;μ_k_n) is described in  <ref>. For the predictive regime, the transport field c_k_n(x;Δμ) in Eq. (<ref>), identified in the offline stage, is evaluated at the unsampled parameter, Δμ = μ^* - μ_k_n.

Finally, the generalized coordinates are obtained by solving the minimization problem (<ref>). The initial guesses for the generalized coordinates at the unsampled parameter a(μ^*)^(0) are obtained by a linear combination of the generalized coordinates at the known snapshots, a(μ_k_n). The generalized coordinates at the known snapshots are given by a(μ_k_n)=e_k_n, where e_k_n∈ℝ^N_k is the k_nth canonical unit vector. The weights of the linear combination are obtained such that they are inversely proportional to | Δμ| and the sum of all weights equals one.





  §.§.§ Choice of basis
 


In traditional projection-based MOR applied to smooth elliptic problems, ROM
performance is expected to improve by increasing the number of bases N_k.
Unfortunately, this is usually not the case when the solutions of interest are
characterized by strong shocks and discontinuities. For
example, Abgrall et al. <cit.> demonstrated that optimal constructions of
such solutions are typically  sparse in the generalized coordinates and local in
the parameter space. This property was also demonstrated in the nozzle problem considered in  <ref>.
Therefore, in this work, we use only a small number of local bases corresponding to N_k nearest neighboring snapshots
selected from the set of N_s snapshots.

For instance, consider the same sampled parameters μ_1<…<μ_5 in Fig. <ref>. For the prediction at a new unsampled parameter μ^*, N_k=2 nearest neighboring snapshots are given by μ_3 and μ_4 which form the snapshots μ_k_1 and μ_k_2, respectively, required for constructing the transported snapshots or local bases. It is noted that the proposed TSMOR approach does not restrict itself to N_k=2 bases and allows for more snapshots to be incorporated. Finally, for μ^* lying outside of the sampled parameter space (i.e. μ^*<μ_1 or μ^*>μ_5), a biased stencil can be used for the selection of the neighboring snapshots.








§ EXTENSION OF TSMOR TO MULTI-DIMENSIONAL PARAMETER VARIATIONS



In this section, the proposed TSMOR approach is detailed for multi-dimensional parameter variations where the snapshots are sampled on a Cartesian grid in parameter space, μ∈𝒟⊂ℝ^N_d.
Without loss of generality, the parameter vector μ∈ℝ^N_d is normalized to fit in a hypercube such that every element of the parameter vector satisfies: -1 ≤μ_i ≤ 1.

TSMOR for multi-dimensional case proceeds exactly in the same manner as that for the one-dimensional case but with the difference that the one-dimensional parameter μ and parameter variation Δμ are replaced by the corresponding vectors μ and Δμ, respectively. Consequently, solutions at an unsampled
parameter μ^* ≠μ_s are approximated as a
linear superposition of parameter dependent basis functions


    u(x;μ^*) ≈u_r = ∑_n=1^N_ku_n(x;μ^*) a_n(μ^*)


where the reduced-basis functions, u_n, correspond to
the transported local snapshots


    u_n = u(x + c_k_n(x,Δμ);μ_k_n)   for   n=1,…,N_k


where the transport field c_k_n(x;Δμ) is evaluated at the unsampled parameter variation, Δμ=μ^*-μ_k_n and the snapshots, u(x,μ_k_n) for
n=1,⋯,N_k, are a subset of the solution snapshots computed offline, μ_k_n∈μ_s. 

Finally, the generalized coordinates a(μ^*) are
identified by minimizing the residual of the Petrov-Galerkin expansion. 


    (v_n, R(u_r;μ^*))_Ω


where (·,·)_Ω is the standard inner product as defined in Eq. (<ref>). 

The proposed TSMOR approach for the multi-dimensional case is also decomposed into the standard 
offline-online strategy where the transport fields are identified offline while the residual of the Petrov-Galerkin projection
is minimized in the online stage.




 §.§ Offline stage



Given the set of snapshots at sampled parameters μ_s, the methodology for identifying the transport fields c_s(x;Δμ), for
s=1,⋯,N_s, where Δμ=μ-μ_s, for each snapshot is
explained in this section.

Similar to TSMOR for one-dimensional parameter variations, the convection field is expressed as a low order expansion in space using spatial basis functions f_p(x). The basis functions f_p(x) are selected a priori such as Chebyshev polynomials
or Fourier modes. For the multi-dimensional case, the dependence of the transport field on parameter variation is modeled as a multi-variate polynomial expansion in parameter space:


    c_s(x;Δμ)=  ∑_p=1^N_p∑_q=1^N_q c_m^sf_p(x)g_q(Δμ)   for   m=1,…,N_p N_q


where c_m^s are the coefficients of the expansion, N_p and N_q are the number of functions f_p(x) and g_q(Δμ), respectively, and g_q(Δμ) is given by:


    g_q(Δμ)= ∏_j=1^N_d(Δμ_j)^h_q,j  for  h_q,j∈ℤ and  q = 1,…,N_q


where Δμ_j for j=1,…,N_d are the vectorial
elements of Δμ. Few terms in Eq. (<ref>) are expanded using Eq. (<ref>) for clarity as:


    c_s(x;Δμ) = ∑_p=1^N_p c_m^s f_p(x)Δμ_1  +  c_m+1^s f_p(x) Δμ_2  +  …  + c_r^s f_p(x) Δμ_1^2  +  c_r+1^s f_p(x) Δμ_1 μ_2 +  …


Since the choice of f_p(x) is problem-dependent, it is detailed in  <ref> with the help of test problems.

The least-squares fitting procedure to identify the coefficients c^s:= c_m^s is similar to the one-dimensional case. Firstly, a subset of solution snapshots u(x;μ_v_n) for n=1,⋯,N_v, where μ_v_n∈μ_s, is chosen. The procedure for selecting these N_v training snapshots, u(x;μ_v_n), from N_s is described in  <ref>. The coefficients c^s are evaluated by minimizing the least-squares training error between the transported snapshot u(x+c_s(x;Δμ);μ_s) and training snapshots u(x;μ_v_n):


    c^smin∑_n=1^N_v( ϵ_n, ϵ_n )_Ω


where ϵ_n(x)=u(x+c_s(x;Δμ);μ_s)-u(x;μ_v_n); Δμ=μ_v_n-μ_s and (·,·)_Ω is the standard inner product as defined in Eq. (<ref>).
The optimization problem (<ref>) is solved N_s times for the coefficients
c^s for s=1,⋯,N_s snapshots. 

In summary, the transport field in Eq. (<ref>) is identified by solving a least-squares fitting problem for the coefficients of the polynomial fitting function.




  §.§.§ Choice of training snapshots



Similar to one-dimensional parameter variations, N_v training snapshots are chosen which correspond to the N_v nearest neighboring snapshots to u(x;μ_s). 

For instance, consider the sampled parameters in a two-dimensional parameter space denoted by `×', `∘' and `∙' as shown in Fig. <ref>. For the evaluation of the transport field for the snapshot at μ_s, N_v=8 nearest neighboring snapshots denoted by `×' are chosen which form the training snapshots μ_v_1 through μ_v_8. The three-dimensional analogue would use N_v=26 neighboring snapshots for the offline training procedure and extends similarly for higher dimensions. Finally, a biased stencil can be used at the boundary of the parameter space for the selection of the training snapshots.









 §.§ Online stage

In this stage, parameter dependent reduced-bases u_n(x;μ^*) at an unsampled parameter μ^* ≠μ_s are constructed and generalized coordinates a(μ^*) are identified.

For a parameter μ^* in the predictive regime, local reduced-bases u_n are evaluated by transporting local snapshots u(x;μ_k_n) for
n=1,⋯,N_k, as shown in Eq. (<ref>), where μ_k_n∈μ_s is a subset of snapshot solutions. The procedure for selecting the snapshots u(x;μ_k_n) is described in  <ref>. For the predictive regime, the transport field c_k_n(x;Δμ) in Eq. (<ref>), identified in the offline stage, is evaluated at the unsampled parameter, Δμ = μ^* - μ_k_n.

Finally, the generalized coordinates are obtained by solving the minimization problem (<ref>). The initial guesses for the generalized coordinates at the unsampled parameter a(μ^*)^(0) are obtained by a linear combination of the generalized coordinates at the known snapshots, a(μ_k_n). The generalized coordinates at the known snapshots are given by a(μ_k_n)=e_k_n, where e_k_n∈ℝ^N_k is the k_nth canonical unit vector. The weights of the linear combination are obtained such that they are inversely proportional to Δμ_2 and the sum of all weights equals one.



  §.§.§ Choice of basis
 


Similar to the one-dimensional case, we use only a small number of local bases corresponding to N_k nearest neighboring snapshots
selected from the set of N_s snapshots. 

For instance, consider the same sampled parameters in a two-dimensional parameter space denoted by `×' and `∘' as shown in Fig. <ref>. For the prediction at a new unsampled parameter μ^*, N_k=4 nearest neighboring snapshots denoted by `×' are chosen which form the snapshots μ_k_1 through μ_k_4 required for constructing the transported snapshots or local bases. It is noted that the proposed TSMOR approach does not restrict itself to N_k=4 basis and allows for more snapshots to be incorporated for this two-dimensional case. The three-dimensional analogue would use N_v=8 neighboring snapshots for basis construction and extends similarly for higher dimensions. Finally, for μ^* lying outside of the sampled parameter space, a biased stencil can be used for the selection of the neighboring snapshots.









§ IMPLEMENTATION OF TSMOR IN DISCRETE FRAMEWORK



As mentioned in remark 3 of  <ref>, FOM solutions are typically available as vectors of state variables specified at various spatial locations. Thus, this section discusses the implementation strategy for applying the proposed TSMOR approach in the discrete framework.

Let w(μ_s) ∈ℝ^N_w denote the discrete FOM solution of the snapshot solution u(x;μ_s). In the continuous form, the transported snapshot u(x+c_s(x;Δμ);μ_s), where Δμ=μ-μ_s, is directly evaluated by computing the solution at the new spatial locations x+c_s(x;Δμ). However, in the discrete form, an interpolation step is required to evaluate the transported snapshots. More specifically, if w(μ_s) is the vector of state variables on a
computational grid with Cartesian coordinates x_i ∈ℝ^N_w
for i =1,2,3, then the transported
solution snapshot is evaluated by interpolating w(μ_s) from x_i+c_s(x_i;Δμ) to the original grid x_i. In the discrete form, c_s(x_i;Δμ)∈ℝ^N_w denotes the transport field in Eq. (<ref>) evaluated at x_i grid points and the transported snapshot from μ_s to μ is denoted by w'(μ_s,μ) ∈ℝ^N_w.

Consequently, the solution at the new parameter using the TSMOR approach is
given by


    w(μ^*) 
        ≈w_r(μ^*)
        =       U(μ^*)a(μ^*)


where U(μ^*) ∈ℝ^N_w× N_k is a matrix whose
columns contain the transported solution snapshots
w'(μ_k_n,μ^*) of the corresponding snapshots
w(μ_k_n). More precisely, U(μ^*)_:,n := w'(μ_k_n,μ^*), where the subscript (:,n) denotes nth column of a matrix. The generalized
coordinates a(μ^*) are
chosen to minimize the residual at the new parameter μ^*:


    amin‖R(U(μ^*)a(μ^*);μ^*)‖ _p


where ‖·‖ _p is the standard ℓ_p-norm. For p=2, this residual minimization (<ref>) is equivalent to a least-squares Petrov-Galerkin projection of the residual equation in the continuous framework given by (<ref>) with test basis v_n=Ju_n where J=∂R/∂a denotes the Jacobian of the residual.

Similarly, the least-squares minimization problem (<ref>) in discrete form can be expressed as:


    c^smin∑_n=1^N_vw'(μ_s,μ_v_n)-w(μ_v_n)^2 _2




Remark 1: Evaluation of the transported snapshot in the discrete form faces a drawback of the possibility of negative volume elements. Since it is usually not possible to evaluate
the solution residual and perform
interpolation on a computational grid with negative volume elements, the training error minimization problem (<ref>) can be augmented with inequality constraints for the volume elements:


    c^smin∑_n=1^N_vw'(μ_s,μ_v_n)-w(μ_v_n)^2 _2 
       subject to   V_j>δ_v,  for  j=1,…,N_V


where V_j are the N_V element volumes and δ_v>0 is the minimal positive
volume. 

For a structured Cartesian grid, this constraint simplifies to
 (x'_i)_j+1 -(x'_i)_j>δ_v for i=1,2,3 and j =
 1,…,N_V, where x'_i=x_i+c_s(x_i;Δμ).

Remark 2: Since the transported snapshot w'(μ_s,μ_v_n) is an approximation of the snapshot w(μ_v_n), w'(μ_s,μ_v_n) should  satisfy the boundary conditions corresponding to the parameter μ_v_n. Generally, for a large class of problems, linear extrapolation at the boundaries of the computational domain (during the interpolation of w(μ_s) from x_i+c_s(x_i;Δμ) to x_i) are found to be a reasonable approximation of the boundary conditions for low to moderate parameter variations. However, approximation of boundary conditions by extrapolation would be a poor choice if the desired boundary conditions have sharp gradients or discontinuities. To tackle this issue, the transport of the boundary nodes can be restricted such that the need for extrapolation is avoided. The transport of the boundary nodes can be restricted by enforcing appropriate conditions on the coefficients c^s of the transport c_s(x_i;Δμ). Since this issue is problem dependent, it is discussed in detail in  <ref> with the help of test problems.

Remark 3: Choice of norms, ℓ_p:
In traditional projection-based MOR, the generalized coordinates are usually
selected to minimize the ℓ_2-norm of the residual. Although this approach
has been demonstrated to work adequately for many applications, in the case when
the FOM is comprised of a system of hyperbolic conservation laws, minimizing the
ℓ_1-norm has been shown to be preferable <cit.>. For our proposed approach, ℓ_1 norm for the residual minimization (<ref>) is found to perform better than other choices of norms.
Unfortunately, the optimal choice of norm remains an open problem.

The classical approach for solving ℓ_1-norm minimization problems involves
recasting the problem as a linear program or, alternatively, solving iteratively
using, for example, Iteratively Reweighted Least Squares
(IRLS) <cit.>. Methodology for
minimizing the ℓ_1-norm is also explained in <cit.>. 





 §.§ Summary of TSMOR


The offline and online stages of the proposed TSMOR approach in the discrete framework are summarized in Algorithms <ref> and <ref>, respectively.









§ HYPER-REDUCTION



In the online stage of TSMOR, evaluation of the residual
R(U(μ^*)a(μ^*)) in Eq. (<ref>) scales with
the size of the FOM, N_w. Furthermore, for every new prediction at μ^*, a new parameter dependent basis U(μ^*) has to be constructed, the computation of which also scales with N_w. Hyper-reduction can significantly reduce these
computational complexities. A review of the state of art hyper-reduction
techniques, such as DEIM <cit.>,
GNAT <cit.> and ECSW <cit.> is provided
in this section. Furthermore, the methodology to equip the proposed TSMOR
approach with hyper-reduction strategies is outlined. 




 §.§ Review of hyper-reduction techniques


In hyper-reduction for traditional projection-based MOR described in  <ref>, the residual R(Ua(μ^*)) is evaluated only at a small subset of n_w
interpolation entries ε⊂{1,⋯,N_w}. The interpolation
matrix Z∈ℝ^N_w × n_w is thus defined as


    Z=[e_ε_1,e_ε_2,⋯,e_ε_n_w]


where e_ε∈ℝ^N_w is the εth canonical unit vector.
R(Ua(μ^*)) is then approximated in a low-dimensional subspace
U_R ∈ℝ^N_w × N_R


    R(Ua(μ^*)) ≈U_R (Z^T U_R)^+ Z^T R(Ua(μ^*))


Finally, equation (<ref>) is projected in a low-dimensional
subspace using state basis U. 
Since the computation of Z^T R(Ua(μ^*)) involves the evaluation of the residual at only n_w grid locations, the resulting minimization problem is independent of the size of FOM, N_w.

Another class of hyper-reduction techniques, for instance ECSW, involves
minimization of the weighted residuals computed only at ε points.
Thus, all the computations are performed on these collocation points and
interpolation of the residuals using basis functions is avoided.




 §.§ Hyper-reduction applied to TSMOR



In this work, we adopt a collocation-based  hyper-reduction approach similar to ECSW. 
More specifically, the residual
minimization problem (<ref>) in this hyper-reduction framework is written
as

    amin‖Z^T R(U(μ^*)a(μ^*);μ^*)‖ _p



Generally, in most CFD problems, the Jacobian matrix is sparse. Hence the
computation of the residuals is dependent only on a few subset of n̂_w
entries ε̂⊂{1,⋯,N_w}. The corresponding
interpolation matrix is denoted as P∈ℝ^N_w ×n̂_w. Thus U(μ^*) a(μ^*) can be
replaced by PP^T U(μ^*) a(μ^*) resulting in the cheap computation of the collocated reduced-basis, P^T U(μ^*), whose columns contain the collocated transported snapshots, i.e., (P^T U(μ^*))_:,n=P^T w'(μ_k_n,μ^*). 

The collocated transported snapshot at the collocation points ε̂ is
computed by interpolating the collocated snapshot, P^T w(μ_k_n), from x̂_i+ĉ_k_n(x̂_i;Δμ) to the original grid x̂_i. 
Here ĉ_k_n(x̂_i;Δμ)=P^T c_k_n(x_i;Δμ) ∈ℝ^n̂_w are the collocated transports and x̂_i= P^T x_i ∈ℝ^n̂_w are the Cartesian coordinates of the original grid at the collocation points.  

To summarize, equation (<ref>) involves the computation of Z^T R(P P^T U(μ^*)a(μ^*)) which necessitates the computation of R(U(μ^*)a(μ^*)) and U(μ^*) only at ε and ε̂ indices respectively, resulting in a reduction of the computational complexity from N_w to n̂_w.




  §.§.§ Identification of collocation points


Various hyper-reduction techniques in literature employ different strategies for
identifying the collocation points or interpolation entries ε.
Generally, these approaches are specific to their respective hyper-reduction
procedure. For instance, these algorithms are based on minimization of the error
in the interpolated snapshots <cit.>, greedy approach
to minimize error associated with gappy-POD projection of
residual <cit.> and solving a sparse non-negative least-squares
(NNLS) problem <cit.>. 
In this work, we employ the standard DEIM Algorithm <ref> to identify
the collocation points.

The DEIM algorithm uses basis functions U_R of the
nonlinear residual to identify the interpolation entries ε. In our
implementation of the algorithm, U_R are the POD basis
of the snapshots of residuals. The snapshots of residuals are collected at each iteration while solving the FOM Eq. (<ref>) during the offline stage.
For boundary value problems, in addition to
DEIM indices, it is important to include inlet/outlet grid points into
ε since these boundary conditions contain vital information about
the dynamics of the problem. Details about inclusion of these boundary points are explained in  <ref> where this topic is covered for each flow problem. Finally, corresponding interpolation entries
ε̂ for computing the residuals can be related to ε
depending on the type of finite-difference/volume/element scheme.








§ COMPARISON WITH PREVIOUS MOR APPROACHES
 


Our proposed approach shares similarities with previous parametric MOR
techniques in the literature. Similar to Reiss et al. <cit.>, Iollo and Lombardi <cit.> and Welper <cit.>, we approximate the sought after solutions
as a superposition of transported snapshots or their basis. However, the methodology to identify the transport fields and corresponding transported snapshots is different where our approach provides the following advantages. 

Shifted-POD method <cit.>, developed for unsteady flow problems, involves the identification of shift velocities based on known or data-driven unsteady transport phenomena. However, this approach has not yet been realized and demonstrated for the reduction of steady flow problems. 

The transport obtained by solving an optimal mass transfer problem <cit.> is optimal in the sense that the displacement of the computational domain to transform one snapshot to the other is minimal. The formulation of the optimal mass transfer problem, however, has the following issues. Firstly, the optimal transports obtained between various pairs of snapshots may not provide an intermediate transported snapshot in the predictive regime that approximately satisfies the physics of the underlying problem. In other words, we seek a physics-based transport field instead of an optimal transport (though an optimal transport may show similarities with physics based transport for a subset of flow problems). In our proposed approach, the transport field is evaluated by a least-squares fitting procedure which minimizes the approximation error of multiple local transported snapshots. This ensures that the transport field captures the physics-based dynamics of snapshot transformation from all the neighboring snapshots. Secondly, it has been demonstrated for image processing applications that for images with sharp features, optimal transport can lead to numerical artifacts at these sharp edges <cit.>, while a Lagrangian approach for solving the optimal mass transfer problem can lead to diffusion of sharp boundaries <cit.>. In the fluid mechanics community, these sharp features can be interpreted as shocks, implying that optimal transport may not be easily extendable for flow problems containing shocks.

In transformed snapshot interpolation (TSI) <cit.>, the transports for approximating the solution in the predictive regime are obtained by interpolating the transports evaluated between various pairs of known snapshots. Our proposed method for identifying the transport fields can be considered as a generalization of TSI, with the difference being that our method is essentially a least-squares fit in contrast to interpolation in TSI. The advantage of fitting as compared to interpolation is that the fitting polynomial can be chosen to contain fewer terms than a corresponding interpolation polynomial. This property allows for easy extension of our method to multi-dimensional parameter variations by mitigating the curse of dimensionality faced by TSI. Furthermore, our residual minimization approach is projection-based in contrast to the interpolation-based approach in TSI. One of the potential advantages of projection-based methods is that they retain the
underlying structure of the dynamical system and thus provide, in principle,
more robust predictive capabilities.    

Finally, similar
to Abgrall et al. <cit.>, we minimize the system residual in the
ℓ_1 norm which has been shown to be beneficial for problems containing
shocks and discontinuities.





§ NUMERICAL EXPERIMENTS



In this section, TSMOR is applied to the steady Euler equations modeling
supersonic flow inside a quasi 1-D nozzle and 2-D flow over a forward facing
step and a nonlinear advection-diffusion equation modeling a jet diffusion flame in a
combustor. These problems are chosen because the steady flow solutions contain
shocks or flame fronts whose spatial locations and orientations are parameter
dependent.  Similar numerical experiments have been studied
by Lucia et al. <cit.>, Mojgani and Balajewicz <cit.>, Zahr et al. <cit.>; Welper <cit.>; and Galbally et al. <cit.> respectively. 

In this
section, the results generated by the proposed TSMOR approach are compared
with traditional projection-based MOR techniques such as LSPG with test subspace v_n =
Ju_n as described in  <ref> where u_n is the POD basis of the snapshot matrix and J is the Jacobian of the residual function. Furthermore, comparison is also made with recent parametric MOR techniques such as L_1-dictionary approach <cit.> where the solution is given by a linear combination of local snapshots or dictionary elements and L_1-norm of the residual is minimized. All results considered in this section are predictive, that is, the predicted
solutions all lie in parameter regions not sampled during the offline training
phase. The performance of these MOR techniques are analyzed by computing the relative error between the predicted and FOM solutions where the error is defined as:


    Error(%)=|| w(μ^*)-w_r(μ^*)  ||_2/||w(μ^*)||_2× 100


where w(μ^*) and w_r(μ^*) are the FOM and predicted solutions using the above-mentioned MOR methods, respectively. All the computations were done in Matlab.




 §.§ Quasi 1-D flow in a converging-diverging nozzle




  §.§.§ Problem description


The 1-D Euler equations in a quasi 1-D converging-diverging nozzle
are considered:

    1/A∂ AF/∂ x=Q x∈[0,L]

where A=A(x) is the area profile and 
 
    w(μ)=[   ρ; ρ u; ρ E ],   F=[      ρ u;  ρ u^2+p; (ρ E+p)u ],   Q=[          0; p/A∂ A/∂ x;          0 ]


with homogeneous Dirichlet boundary conditions ρ(0;μ)=1,  p(0;μ)=1 and
p(L;μ)=0.7.

The boundary conditions are chosen such that a shock is formed in the diverging
section of the nozzle. Length of the nozzle L is 10 units. Area profile of the
converging-diverging nozzle is parabolic with equal inlet and outlet area, A(0)=A(L)=3, and
the throat is located at L/2. For this problem, the throat area μ=A(L/2)
is the parameter of interest. Steady state solutions are obtained by
discretizing the corresponding governing equations in space using a central
second-order finite difference scheme on a uniform grid which is divided into 1000 grid points with grid spacing
Δ x=0.01. A first-order accurate artificial viscosity scheme using
ν=Δ x / 2 is used to stabilize the solution. The resulting nonlinear
system of algebraic equations is solved in Matlab using the built-in 
algorithm. Fig. <ref> shows the steady density solutions for different
values of throat area μ.




  §.§.§ Implementation of TSMOR


A snapshot matrix M containing 4 snapshots at parameters μ_s=[0.5, 0.875, 1.25, 1.625] is generated. The coefficients of the polynomial expansion (<ref>) c^s for each snapshot are computed offline by solving the training error minimization (<ref>). The basis f_p(x) are chosen to be Fourier sine series with m modes:


    f(x)=[1, sin(πx/L), …, sin((m-1)πx/L) ]


For this simple problem, g_q(Δμ)  is given by: 

    g(Δμ)=Δμ


The interpolation from the
transported grid to the original grid for calculating the transported snapshots
was performed using  algorithm. The training error minimization (<ref>) is solved using  algorithm.

First, convergence of the proposed TSMOR approach with respect to the number of Fourier modes m is studied by predicting new solutions in the predictive regime μ^*. Fig. <ref> shows the mean and maximum relative error in the TSMOR solutions predicted at two uniformly distributed parameters in every interval of μ_s for different number of Fourier modes m. It can be seen that as the number of Fourier modes increases, the error converges to a low value of 0.27%. The TSMOR convergence plot is compared to the convergence of LSPG solutions with respect to the number of POD basis k of the snapshot matrix M. Similar to Fig. <ref>, Fig. <ref> displays the relative errors in the LSPG solutions predicted at the same set of parameters for different number of bases. It can be observed that, the error converges only to 7.16% even though all 4 POD bases were used for prediction. Hereafter, all the TSMOR predicted results presented for this problem are produced with 3 Fourier modes.







Next, the performance of the proposed TSMOR approach is compared to several existing MOR techniques. Fig. <ref> illustrates the predictive capabilities of several MOR
approaches for the parameters μ^*=1.0 and μ^*=1.5. The FOM density solution is given by the gray
lines while the new proposed TSMOR approach using 2 local bases corresponding to two nearby snapshots is given by the red lines. The blue lines correspond to the solution obtained by LSPG method using 4 POD modes of the snapshot matrix M. The green lines correspond to L_1-dictionary approach using 2 local snapshots or dictionary elements. The proposed TSMOR
approach reproduces the solution remarkably well. In contrast, LSPG does not predict the correct shock
location while L_1-dictionary solutions are dominated by staircase shock type errors.

Next, the TSMOR approach is equipped with the hyper-reduction strategy mentioned in  <ref>. 
First, 30 collocation points are
obtained by employing the DEIM algorithm. Second, these points are augmented
with inlet and outlet points, i=1, and i=N_w, respectively, if not already included in those 30 DEIM collocation points.
Finally, an additional n̂_w≈ n_w × 2=64
points ε̂ are included to enable the evaluation of the residuals via
the central finite difference scheme. 
Fig. <ref> illustrates the predictive capabilities of hyper-reduced TSMOR (TSMOR+HR)
for the same parameters μ^*=1.0 and μ^*=1.5. 
The hyper-reduced TSMOR density solutions are compared with non-hyper-reduced LSPG and L_1-dictionary approaches. As before, excellent agreement with the FOM solution is demonstrated by hyper-reduced TSMOR.







In Fig. <ref>, relative errors between the FOM solutions and
predicted solutions using hyper-reduced TSMOR, non-hyper-reduced LSPG and L_1-dictionary approaches are illustrated across the entire parameter range of interest.  For
this case, predictions are made at two uniformly distributed parameters in every interval of μ_s. It can be observed that the solutions predicted using hyper-reduced TSMOR have an average error of only 0.27% as compared to 7.22% in LSPG and 5.88% in L_1-dictionary approach. Thus, for all parameters considered, the TSMOR approach significantly
outperforms LSPG and L_1-dictionary methods.  
Finally, wall-times and speed-ups for the FOM and the online stage of hyper-reduced TSMOR are illustrated in
Fig. <ref> and Fig. <ref>, respectively. Here, speed-up is defined as the ratio of wall-times of FOM to the online stage of hyper-reduced TSMOR. TSMOR+HR
delivers a speed-up of approximately four orders of magnitude across the entire
parameter range while maintaining a high level of accuracy. 














 §.§ Supersonic flow over a forward facing step



  §.§.§ Problem description


This problem consists of a supersonic flow over a 2-D forward facing step in a
wind tunnel setup with walls on top and bottom as described
by <cit.> and also shown in Fig. <ref>. Length (L)
and height (H) of the wind tunnel are 3 units and 1 unit respectively. The step
has a height of 0.2 units and is located at the bottom wall starting at 0.6
units from the left-end of the tunnel. 2-D Euler equations governing the
supersonic flow over a forward facing step are:

    ∂F/∂ x+∂G/∂ y=0,  x∈[0,L],  y∈[0,H]

where

    w(μ)=[   ρ; ρ u; ρ v; ρ E ];
    F=[      ρ u;  ρ u^2+p;     ρ uv; (ρ E+p)u ];
    G=[      ρ v;     ρ uv;  ρ v^2+p; (ρ E+p)v ]


with homogeneous Dirichlet inlet boundary conditions ρ(0,y;μ)=1.4, p(0,y;μ)=1,  u(0,y;μ)=μ and v(0,y;μ) = 0;
where μ is the inlet supersonic Mach number which is taken to be the varying
parameter of interest. No penetration solid wall boundary conditions are imposed
on the top, bottom and step wall surfaces. The equations are discretized in
space using a second-order, central finite difference scheme on a uniform
Cartesian grid which is divided into 0.48 million grid points with Δ x=Δ
y=0.025. The solutions are stabilized using a first-order artificial viscosity
scheme where the artificial viscosity is set to be ν=Δ x / 0.8. The
resulting equations are solved by marching to steady state using a second-order
Strong Stability Preserving (SSP) Runge–Kutta scheme <cit.>.
Fig. <ref> shows the steady state density contour for inlet Mach number μ=3.3,
while the corresponding 1-D plots at y=0.7 for inlet Mach numbers μ=3.3, 3.6, and 3.9 are shown in Fig. <ref>. 







  §.§.§ Implementation of TSMOR






A snapshot matrix M containing 5 snapshots at parameters μ_s=[3.3, 3.45, 3.6, 3.75, 3.9] is generated. 
The coefficients of the polynomial expansion (<ref>) c^s for each snapshot are computed offline by solving the training error minimization (<ref>). The bases of the polynomial expansion, f_x_p(x) and f_y_p(y), for the x and y transport fields, c_s_x(x,Δμ) and c_s_y(y,Δμ), respectively, are chosen to be Fourier sine series with m modes each:


    f_x(x)=[1, sin(πx/L), …, sin((m-1)πx/L) ]
    
    f_y(y)=[1, sin(πy/H), …, sin((m-1)πy/H) ]


and g_q(Δμ) is given by:


    g(Δμ)=[Δμ, Δμ^2 ]


The interpolation from the
transported grid (x+c_s_x(x;Δμ),y+c_s_y(y;Δμ)) to the original grid (x,y) for calculating the transported snapshots
was performed using  algorithm in Matlab. The training error minimization (<ref>) is solved using the  algorithm.

First, convergence of the proposed TSMOR approach with respect to the number of Fourier modes m is studied by predicting new solutions in the predictive regime μ^*. Fig. <ref> shows the mean and maximum relative error in the TSMOR solutions predicted at two uniformly distributed parameters in every interval of μ_s for different number of Fourier modes m. It can be seen that as the number of Fourier modes increases, the error converges to a low value of 0.54%. The TSMOR convergence plot is compared to the convergence of LSPG approach with respect to the number of POD basis k of the snapshot matrix M. Similar to Fig. <ref>, Fig. <ref> displays the relative errors in the LSPG solutions predicted at the same set of parameters for different number of bases. It can be observed that, the error converges only to 5.24% even though all 5 POD basis were used for prediction. Hereafter, all the TSMOR predicted results presented for this problem are produced with 9 Fourier modes.

Next, the performance of the proposed TSMOR approach is compared to several existing MOR techniques. Fig. <ref> illustrates the predictive capabilities of several MOR
approaches for the parameter μ^*=3.4. The FOM density solution is
shown in Fig. <ref> while the new proposed TSMOR solution using 2 local bases corresponding to two nearby snapshots is shown in
Fig. <ref>. Fig. <ref> corresponds to the solution obtained by LSPG using 4 POD modes of the snapshot matrix M. Fig. <ref> corresponds to L_1-dictionary approach using 2 local bases or dictionary elements. To compare these solutions qualitatively, the predicted density distributions at various y-locations for the FOM, TSMOR, LSPG and L_1-dictionary approaches are shown in
Fig. <ref>. It can be observed that the proposed TSMOR approach significantly outperforms
LSPG and L_1-dictionary methods.










Finally, the relative solution error between the FOM solution and
predicted solution using TSMOR, LSPG and L_1-dictionary approaches across the entire parameter range of interest is
given in Fig. <ref>. For
this case, predictions are made at two uniformly distributed parameters in every interval of μ_s. It can be observed that the solutions predicted using TSMOR have an average error of only 0.6% as compared to 5.9% in LSPG and 5.0% in L_1-dictionary approach.









 §.§ Jet diffusion flame in a combustor



  §.§.§ Problem description


This problem consists of jets of fuel and oxidizer injected into a combustion chamber as shown in Fig. <ref>. Length (L) and height (H) of the chamber are 18 mm and 9 mm respectively. The width of fuel and oxidizer inlets are denoted by H_f and H_o respectively. Inside the chamber, the fuel and oxidizer diffuse to form a diffusion flame where the combustion reaction is governed by an advection-diffusion type governing equation:

    ∇·(Ww(μ))-∇·(ν∇w(μ))+F(w(μ))=0  x∈[0,L],  y∈[0,H]


where the state variable w(μ) represents the concentration of fuel in the chamber, W=W_xî + W_yĵ is the velocity field and ν is the diffusion coefficient. The nonlinear reaction term F(w(μ)) is of Arrhenius type given by:

    F(w(μ))=Aw(μ)(c-w(μ))e^-E/(d-w(μ))

where c, d, A and E are constants and (c-w(μ)) represents the oxidizer concentration. Dirichlet boundary condition is prescribed at the inlet:

    w(0,y;μ)=0  y∈[0,H_o)


    w(0,y;μ)=c  y∈[H_o,H_o+H_f]


    w(0,y;μ)=0  y∈(H_o+H_f,H]


whereas homogeneous Neumann boundary conditions are prescribed on other boundaries.



In this problem, we consider parameter variations in three dimensions where μ=[W_y, ln(A), r] consists of the velocity field in y-direction W_y, Arrhenius parameter ln(A) and ratio of fuel to oxidizer inlet widths r=H_f/H_o. W_y, ln(A) and r influence the direction, length and width of the flame, respectively. The values of the remaining constants are: W_x=0.17 m/s, ν=5× 10^-6m^2/s, c=0.2, d=0.24 and E=0.1091. The equations are discretized in
space using a second-order, central finite difference scheme on a uniform
Cartesian grid which is divided into 1 million grid points with Δ x=Δ
y=10^-5mm. A fine discretization is chosen to demonstrate the speedup associated with the hyper-reduction approach. The
resulting equations are solved using Newton's method until convergence of 10 orders of magnitude. Fig. <ref> shows the fuel concentration contours computed at eight different corners of the 3-D parameter space. Clear and distinct flame fronts having different directions, lengths and widths can be observed in these plots.







  §.§.§ Implementation of TSMOR

A snapshot matrix M containing 48 snapshots is generated on a 4× 4 × 3 grid in parameter space 𝒟 at parameters W_y× ln(A) × r ≡ [-0.02,0,0.02,0.04]× [7.0,7.2,7.4,7.6] × [0.467, 0.846, 1.364]. The coefficients of the polynomial expansion (<ref>) c^s for each snapshot are computed offline by solving the training error minimization (<ref>). Since the parameter W_y changes the angle of the flame, a strong coupling between x and y coordinates is required. As a result, the basis of the polynomial expansion, f_x_p(x,y) and f_y_p(x,y), for the x and y transport fields, c_s_x(x,y;Δμ) and c_s_y(x,y;Δμ), respectively, are chosen to be a combination of coupled and decoupled Fourier sine series with m modes. Here we take f_x_p(x,y)=f_y_p(x,y)=f_p(x,y):


    m=1   f(x,y)=sin(πx/2L)
       m=2   f(x,y)=[f(x,y)|_m=1, sin(πx/L), sin(πy/H) ]
       m=3   f(x,y)=[f(x,y)|_m=2, sin(2πx/L), sin(πx/L)sin(πy/H), sin(2πy/H) ]
       m=4   f(x,y)=[f(x,y)|_m=3, sin(3πx/L), sin(2πx/L)sin(πy/H), sin(πx/L)sin(2πy/H), sin(3πy/H) ]


Note that the basis functions in Eq. (<ref>) is cumulative in the sense that for m>1, the basis includes the basis corresponding to the previous modes as well. 
Since the inlet boundary condition for this problem is discontinuous, boundary conditions on the transported snapshots are enforced by enforcing appropriate conditions on the transports, as described in remark 2 of  <ref>.
More specifically, to avoid extrapolation in the x-direction, the following condition is imposed on the x-transports at the inlet: c_s_x(x=0,y;Δμ)=0. This condition is enforced by setting the coefficients of decoupled basis sin(i π y/H)=0 for i=1,2,3 in the expansion of c_s_x(x,y;Δμ). Note that this condition is not enforced on y-transports c_s_y(x,y;Δμ).
For this multi-dimensional parameter problem, g_q(Δμ)  is given by:


    g(Δμ)=[Δμ_1, Δμ_2, Δμ_3, Δμ_1Δμ_2, Δμ_2Δμ_3, Δμ_3Δμ_1, Δμ_1Δμ_2Δμ_3]



Transformation of the structured Cartesian grid based on the coupled
transports (<ref>) leads to a non-tensor grid. Hence,
interpolation from the non-tensor transported grid to the original grid
must be performed using the  algorithm in Matlab. 
Unfortunately, this interpolation scheme is computationally prohibitive for large
grids. Therefore, during the offline stage, the snapshots are uniformly
downsampled by a factor Δ s in x and y coordinates leading to a coarser grid with N_Δ s points. A grid convergence study is conducted to identify Δ s and N_Δ s.
Fig. <ref> plots the relative error in the TSMOR solutions predicted at μ^*=[-0.01,7.1,0.644] for increasing N_Δ s. It can be observed that the reduction in error is minimal after 3000 grid points which corresponds to a downsampling factor of 25. Hence, the snapshots are uniformly
downsampled by a factor of 25. Note that this downsampling process is conducted only during the offline stage.









Next, convergence of the proposed TSMOR approach with respect to the number of Fourier modes m is studied by predicting new solutions in the predictive regime μ^*. Fig. <ref> shows the mean and maximum relative error in the TSMOR solutions predicted at 18 uniformly distributed parameters in the parameter space 𝒟 for different numbers of Fourier modes m. It can be observed that as the number of Fourier modes increases, the error converges to a low value of 1.91%. The TSMOR convergence plot is compared to the convergence of LSPG approach with respect to the number of POD basis k of the snapshot matrix M. Similar to Fig. <ref>, Fig. <ref> displays the relative errors in the LSPG solutions predicted at the same set of parameters for different numbers of bases. It can be observed that, the error converges only to 30% even though all 48 POD basis were used for prediction. Hereafter, all the TSMOR predicted results presented for this problem are produced with 4 Fourier modes.

Next, the TSMOR approach is equipped with the hyper-reduction strategy mentioned in  <ref>. First, 70 collocation points are obtained by employing the DEIM algorithm. Second, these points are augmented with 30 uniformly distributed inlet points. Finally, an additional n̂_w≈ n_w × 4=400 points ε̂ are included to enable the evaluation of the residuals via the central finite difference scheme.

To demonstrate the performance of the proposed TSMOR approach, solutions are predicted at four different predictive regimes in the parameter space 𝒟 and compared with FOM and several other MOR techniques. Parameters corresponding to these test cases are tabulated in Table. <ref>. Figs. <ref>- <ref> illustrate the predictive capabilities of several MOR approaches for these test cases. Contour plots at contour levels w(μ^*)=0.018 and w(μ^*)=0.15 are displayed in these figures. The FOM is given by the gray lines while the new proposed hyper-reduced TSMOR approach using 8 local basis is given by the red lines. The blue lines correspond to the solution obtained by LSPG using 48 POD modes of the snapshot matrix M. The green lines correspond to L_1-dictionary approach using 8 local basis or dictionary elements. The proposed TSMOR approach predicts a solution which accurately matches the FOM solutions in all the four cases. In contrast, both LSPG and L_1-dictionary methods fail to capture the flame-front. 





The ROM errors for the various MOR approaches are summarized in
Table <ref>. Solutions predicted using TSMOR
have an average error of only 1.92% as compared to 29.02% in LSPG and 17.53% in
the L_1-dictionary approach. 
The table also provides the error in the POD-projected solution which is
obtained by projecting the FOM onto 48 POD basis of the snapshot matrix M.
The average error for the POD-projected solution is 8.65%. Thus, TSMOR significantly
outperforms LSPG and L_1-dictionary approaches and it is 3-4 times better than
POD-projected solutions. 









Finally, wall-times and speed-ups for the FOM and the online stage of hyper-reduced
TSMOR are given in Table. <ref>. Here, speed-up is defined as the ratio of wall-times of FOM to the online stage of hyper-reduced TSMOR. TSMOR+HR delivers a speed-up of
two orders of magnitude.


















§ LIMITATIONS OF TSMOR



TSMOR can be expected to provide significant improvements in accuracy over traditional projection-based ROMs for a large class of problems characterized by shocks and sharp gradients whose spatial locations and orientations are strongly parameter dependent. There are, however, a large class of important and very challenging problems where TSMOR – and indeed most other related methodologies in the literature – can be expected to fail.

For example, similar to most other model reduction approaches for parameter variations in the literature, TSMOR requires that the topology of the underlying solution varies smoothly with respect to the parameters. In other words, TSMOR cannot be expected to provide accurate predictions of bifurcations. In the case when the changes in solutions with respect to parameter variations are sharp yet still smooth and continuous, an adaptive parameter sampling <cit.> strategy during the offline stage can be expected to mitigate this issue.

Finally, since both the online and offline stages of TSMOR require solutions of non-convex optimization problems (Eq. (<ref>) and (<ref>), respectively), convergence to local, sub-optimal minima can be a real problem. However, in principle, this issue can usually be avoided by implementing a global optimization scheme, or – more typically – refining and/or adaptively sampling the parameter space during the offline stage.




§ CONCLUSIONS



In this manuscript, a transported snapshot model order reduction (TSMOR) method
for predicting new parametric steady state solutions containing moving shocks
and discontinuities is presented. In this method, the solution is approximated by a linear
combination of spatially transported snapshots. The transports are assumed to be smooth in parameter as well as physical space, and hence approximated as a low-order polynomial expansion. The coefficients of the polynomial expansion are obtained by solving a training error minimization problem in the offline stage. The generalized coordinates are derived by solving a residual minimization problem in the online stage. TSMOR is also integrated with hyper-reduction methods to reduce the computational complexity of evaluating the nonlinear residual and parameter dependent basis in the online stage. 

Numerical experiments consist of a 1-D converging-diverging
nozzle problem with throat area as the parameter, a supersonic flow over a
forward facing step with inlet Mach number as the parameter and a multi-dimensional parametric combustion problem with three parameters influencing the length, direction and width of the diffusion flame. For all parameters
considered, TSMOR is demonstrated to significantly outperform traditional
approaches such as those based on linear compression
schemes e.g. LSPG and more recent local basis approaches such as L_1-dictionary. Furthermore, speed-up of two and four orders of magnitude for the combustion and nozzle problems are achieved, respectively.




§ ACKNOWLEDGEMENT

This material is based upon the work supported by the Air Force Office of Scientific Research
under Grant No. FA9550-17-1-0203. 


WileyNJD-AMA



