## 1 Introduction

Rayleigh–Taylor instability is a classical hydrodynamic instability occuring at the interface between two fluids of different densities. A denser fluid is initially supported on top by a less dense fluid with a counterbalancing pressure gradient under the effects of a gravitational field. It was first introduced by Rayleigh [rayleigh1882investigation] in 1882 and further investigated by Taylor [taylor1950instability]. In recent years, it has been applied to various disciplines of science and engineering, including astrophyics [livio1980rayleigh, yamada1991rayleigh, blondin20011rayleigh, ribeyre2004compressible], atmospheric sciences and climate [keskinen1981nonlinear, huang1993nonlinear], geophysics [wilcock1991rayleigh, plag1995rayleigh, conrad1997growth], and fusion energy [betti1998growth, remington2019rayleigh]. The readers are referred to [chandrasekhar1961hydrodynamic, sharp1984overview, casner2019icf] for a comprehensive overview of Rayleigh–Taylor instability.

In recent years, numerous research efforts have been devoted to modeling the growth of Rayleigh–Taylor instability mathematically and numerically, so as to provide better understanding of the nonlinear dynamics. Linear stability theory can be used to analytically derive the exponential growth in the amplitude in the subwavelength regime in the short-time dynamics [taylor1950instability]. However, as the nonlinearity dominates and the dynamics become chaotic in the later stages, numerical simulation plays an important role in the studies of the long-time dynamics [tryggvason1990computations, youngs1991three, dimonte2004comparative, lee2013numerical, zanella2020two]

. Depending on the nature of the fluids, the physical process is governed by a set of conservation laws and mathematically modeled by a system of partial differential equations. In this work, we consider the scenario of a single-mode two-dimensional Rayleigh–Taylor instability between two compressible and isotropic ideal gases with different densities, in a parametric setting. The Euler equation is used to model the compressible gas dynamics in a complex multimaterial setting. Moreover, we study the effects of Atwood number as a problem parameter on the physical quantities. The governing equations are numerically solved in a moving Lagrangian frame, where the computational mesh is moved along with the fluid velocity. Our work is built on

[dobrev2012high], where a general framework of high-order curvilinear finite elements and adaptive time stepping of explicit time integrators is proposed for numerical discretization of the Lagrangian hydrodynamics problem over general unstructured two-dimensional and three-dimensional computational domains. The method shows great capability and lots of advantages, including high-order accuracy, total energy conservation, and the Rankine-Hugoniot jump conditions fulfilled at shock boundaries. However, forward simulations of Lagrangian hydrodynamics can be computationally very expensive, especially in long-time simulation of highly nonlinear problems like Rayleigh–Taylor instability. It is therefore desirable to develop efficient techniques for accelerating the computationally expensive simulations. In this work, we propose new reduced order model (ROM) techniques for Rayleigh–Taylor instability in compressible Euler equations.A reduced order model has nowadays become a popular and actively researched computational technique to reduce the computational cost of simulations while minimizing the error introduced in the reduction process. Many physics-based models are complex and nonlinear, and are formulated on large spatio-temporal domains, in which the computational cost can be prohibitively high. It may take a long time to run one forward simulation even with high performance computing. In decision-making applications where multiple forward simulations are needed, such as parameter study, design optimization [wang2007large, de2020three, de2018adaptive, white2020dual], optimal control [choi2015practical, choi2012simultaneous], uncertainty quantification [smith2013uncertainty, biegler2011large], and inverse problems [galbally2010non, biegler2011large], the computationally expensive simulations are not desirable. To this end, a reduced order model can be useful to obtain sufficiently accurate approximate solutions with considerable speed-up compared to a corresponding full order model (FOM).

Many model reduction schemes have been developed to reduce the computational cost of simulations while minimizing the error introduced in the reduction process. Most of these approaches seek to extract an intrinsic solution subspace for condensed solution representation by a linear combination of reduced basis vectors. The reduced basis vectors are extracted from performing proper orthogonal decomposition (POD) on the snapshot data of the FOM simulations. The number of degrees of freedom is then reduced by substituting the ROM solution representation into the (semi-)discretized governing equation. These approaches take advantage of both the

known governing equation and the solution data generated from the corresponding FOM simulations to form linear subspace reduced order models (LS-ROM). Example applications include, but are not limited to, the nonlinear diffusion equations [hoang2021domain, fritzen2018algorithmic], the Burgers equation and the Euler equations in small-scale [choi2019space, choi2020sns, carlberg2018conservative], the convection–diffusion equations [mojgani2017lagrangian, kim2021efficientII], the Navier–Stokes equations [xiao2014non, burkardt2006pod], rocket nozzle shape design [amsallem2015design], flutter avoidance wing shape optimization [choi2020gradient], topology optimization of wind turbine blades [choi2019accelerating], lattice structure design [mcbane2021component], porous media flow/reservoir simulations [ghasemi2015localized, jiang2019implementation, yang2016fast, wang2020generalized], computational electro-cardiology [yang2017efficient], inverse problems [fu2018pod], shallow water equations [zhao2014pod, cstefuanescu2013pod], Boltzmann transport problems [choi2021space], computing electromyography [mordhorst2017pod], spatio-temporal dynamics of a predator–prey system [dimitriu2013application], acoustic wave-driven microfluidic biochips [antil2012reduced], and Schrödinger equation [cheng2016reduced]. Survey papers for the projection-based LS-ROM techniques can be found in [gugercin2004survey, benner2015survey].In spite of successes of the classical LS-ROM in many applications, these approaches are limited to the assumption that the intrinsic solution space falls into a subspace with a small dimension, i.e., the solution space with a Kolmogorov -width decaying fast. This assumption is violated in advection-dominated problems, due to features such as sharp gradients, moving shock fronts, and turbulence, which hinder these model reduction schemes from being practical. Our goal in this paper is to develop an efficient reduced order model for hydrodynamics simulation with advection-dominated solutions. Some reduced order model techniques for hydrodynamics or turbulence models in the literature include [mou2020data, parish2017non, gadalla2020comparison, bergmann2009enablers, osth2014need, baiges2015reduced, san2018extreme], which are mostly built on the Eulerian formulation, i.e., the computational mesh is stationary with respect to the fluid motion. In contrast, numerical methods in the Lagrangian formulation, which are characterized by a computational mesh that moves along with the fluid velocity, are developed for better capturing the shocks and preserving the conserved quantities in advection-dominated problems. It therefore becomes natural to develop Lagrangian-based reduced order models to overcome the challenges posed by advection-dominated problems. Some existing work in this research direction include [mojgani2017lagrangian, lu2020lagrangian], where a Lagrangian POD and dynamic mode decomposition (DMD) reduced order model are introduced respectively for the one-dimensional nonlinear advection-diffusion equation. We remark that there are some similarities and differences between our work and [mojgani2017lagrangian] in using POD for developing Lagrangian-based reduced order models. It is important to note that our work is based on the more complicated and challenging two-dimensional compressible Euler equations. The reduced bases are built independently for each state variable in our work, while a single basis is built for the whole state in [mojgani2017lagrangian]

. Furthermore, we introduce the solution manifold decomposition concept so as to ensure adequate ROM speed-up by classifying solution snapshot samples for constructing ROM spaces with low dimension and assigning the appropriate ROM to be used in the time marching.

Recently, there have been many attempts to develop efficient ROMs for the advection-dominated or sharp gradient problems. The attempts can be divided mainly into two categories. The first category enhances the solution representability of the linear subspace by introducing some special treatments and adaptive schemes. A dictionary-based model reduction method for the approximation of nonlinear hyperbolic equations is developed in [abgrall2016robust], where the reduced approximation is obtained from the minimization of the residual in the norm for the reduced linear subspace. In [carlberg2015adaptive], a fail-safe -adaptive algorithm is developed. The algorithm enables ROMs to be incrementally refined to capture the shock phenomena which are unobserved in the original reduced basis through a-posteriori online enrichment of the reduced-basis space by decomposing a given basis vector into several vectors with disjoint support. The windowed least-squares Petrov–Galerkin model reduction for dynamical systems with implicit time integrators is introduced in [parish2019windowed, shimizu2020windowed], which can overcome the challenges arising from the advection-dominated problems by representing only a small time window with a local ROM. Our previous work [copeland2022reduced] adopts a similar approach for Lagrangian hydrodynamics. Another active research direction is to exploit the sharp gradients and represent spatially local features in ROM, such as the online adaptivity bases and adaptive sampling approach [peherstorfer2018model] and the shock reconstruction surrogate approach [constantine2012reduced]. In [taddei2020space], an adaptive space-time registration-based model reduction is used to align local features of parameterized hyperbolic PDEs in a fixed one-dimensional reference domain. Some new approaches have been developed for aligning the sharp gradients by using a superposition of snapshots with shifts or transforms. In [reiss2018shifted], the shifted proper orthogonal decomposition (sPOD) introduces time-dependent shifts of the snapshot matrix in POD in an attempt to separate different transport velocities in advection-dominated problems. The practicality of this approach relies heavily on accurate determination of shifted velocities. In [rim2018transport], an iterative transport reversal algorithm is proposed to decompose the snapshot matrix into multiple shifting profiles. In [welper2020transformed], inspired by the template fitting [kirby1992reconstructing]

, a high resolution transformed snapshot interpolation with an appropriate behavior near singularities is considered.

The second category replaces the linear subspace solution representation with the nonlinear manifold, which is a very active research direction. Recently, a neural network-based reduced order model is developed in

[lee2020model] and extended to preserve the conserved quantities in the physical conservation laws [lee2019deep]. In these approaches, the weights and biases in the neural network are determined in the training phase, and existing numerical methods, such as finite difference and finite element methods, are utilized. However, since the nonlinear terms need to be updated every time step or Newton step, and the computation of the nonlinear terms still scale with the FOM size, these approaches do not achieve any speed-up with respect to the corresponding FOM. Recently, Kim, et al., have achieved a considerable speed-up with the nonlinear manifold reduced order model [kim2021fast, kim2020efficient], but it was only applied to small problems. Manifold approximations via transported subspaces in [rim2019manifold] introduced a low-rank approximation to the transport dynamics by approximating the solution manifold with a transported subspace generated by low-rank transport modes. However, their approach is limited to a one-dimensional problem setting. In [rim2020depth], a depth separation approach for reduced deep networks in nonlinear model reduction is presented, in which the reduced order model is composed with hidden layers with low-rank representation.The method presented in this paper belongs to the first category. The underlying concept is to build small and accurate projection-based reduced-order models by decomposing the solution manifold into submanifolds. These reduced-order models are local in the sense that each of them will be valid only over a certain sub-interval of the temporal domain. The appropriate local reduced order model is chosen based on the current state of the system, and all the local reduced order models will cover the whole time marching in the online phase. The concept of a local reduced order model was introduced in [washabaugh2012nonlinear, amsallem2012nonlinear], where unsupervised clustering is used for the solution manifold decomposition. In our work, the solution manifold decomposition is based on a suitably defined physics-based indicator, which generalizes the time-windowing ROM approach for Lagrangian hydrodynamics simulations in our previous work [copeland2022reduced]. The idea of time-windowing ROM is to construct temporally-local ROM spaces which are small but accurate within a short period in advection-dominated problems. The method provided tremendous speed-up and accurate approximated solutions in various benchmark experiments including Sedov blast, Gresho vortices, Taylor-Green vortices, and triple-point problems. However, in order to approximate the solutions of Rayleigh–Taylor instability over a range of Atwood number, the time-windowing approach is insufficient to provide good approximations, since the penetration speed of the fluid interface varies with the Atwood number. In this work, we modify our ROM approach by introducing new techniques to capture the fast-moving shock boundaries, which provide significant improvements in Rayleigh–Taylor instability.

Our contribution in this paper is to propose an efficient model order reduction scheme for accelerating the simulation of Rayleigh–Taylor instability in compressible gas dynamics with varying Atwood number. Similar to [copeland2022reduced], the idea is to construct temporally-local ROM spaces which are small but accurate within a short period in advection-dominated problems to achieve a good speed-up and solution accuracy. Proper orthogonal decomposition (POD) is used to extract the dominant modes in solution representability, the solution nonlinear subspace (SNS) method is used to establish the subspace relations and construct the nonlinear term bases as in [choi2020sns], and an oversampling discrete empirical interpolation method (DEIM) serves as a hyper-reduction technique to reduce the complexity due to the nonlinear terms in the governing equations. The difference between this work and [copeland2022reduced] is the introduction of a new temporal domain partition scheme that allows a smaller dimension of reduced basis in each temporal subdomain and use of temporally-local ROMs at different Atwood numbers for better solution accuracy. In this work, we propose a general framework to decompose the solution manifold by a suitably defined indicator. The indicator is used to classify solution snapshot samples for constructing ROM spaces with low dimension and to assign the appropriate ROM to be used in the time marching. We consider two such indicators, namely the physical time and the penetration distance, and compare the performance in speed-up and solution accuracy. As we will see in our numerical experiments, the penetration distance is a good alternative indicator to resolve the deficiency of degenerating solution accuracy with respect to Atwood number observed in the time-windowing approach. We remark that, by taking the advantage of Lagrangian formulation of the Euler equations, the penetration distance is easily accessible and plays an important role in the newly proposed distance-windowing approach.

The rest of the paper is organized as follows. In Section 2, we introduce the governing equations and the numerical discretization which will be used as the full order model. Next, a projection-based ROM is described in Section 3. In Section 4, a general framework of solution manifold decomposition is introduced, and two practical examples will be discussed. Numerical results are presented in Section 5. Finally, conclusions are given in Section 6.

## 2 Numerical modeling

In this section, we describe a direct numerical simulation methodology for Rayleigh–Taylor instability in compressible gas dynamics, using a finite element method, which will serve as the full order model in this work. We first present the governing equations, initial and boundary conditions which give rise to Rayleigh–Taylor instability. We also present the key ingredients of the numerical solver.

### 2.1 Governing equations

We consider the two-dimensional mathematical model of Rayleigh–Taylor hydrodynamic instability caused by a gravitational field acting on compressible isotropic ideal gas with stratified densities. Suppose is the final time of the hydrodynamic process. For , let denote a continuous medium which deforms in time, with the initial configuration at being a rectangular domain with aspect ratio 1:4. The compressible gas dynamics is modeled by the Euler equations in a Lagrangian reference frame [harlow1971fluid], driven by a constant gravitational acceleration :

(2.1) | ||||||

Here, is the material derivative, denotes the density of the fluid,

denotes the deformation stress tensor,

denotes the Lagrangian position, denotes the velocity, and denotes the internal energy per unit mass. Each of these physical quantities depends on the time and the Eulerian coordinates of the particle . In gas dynamics, the stress tensor is isotropic, and we write , where denotes the thermodynamic pressure, and denotes the artificial viscosity stress. Due to the ideal gas assumption, the thermodynamic pressure is related to the density and the internal energy by the equation of state(2.2) |

where the adiabatic index is . For , the system is prescribed with a boundary condition , where is the outward normal unit vector on the domain boundary. In order to close the system, an initial condition needs to be imposed. We assume the gas is initially stratified at the interface , supported by a counterbalancing pressure gradient and perturbed by a velocity. In mathematical terms, at the time , for , we take

(2.3) |

where is the density ratio of the fluids. In the context of hydrodynamic instability, we define the dimensionless Atwood number by

(2.4) |

which increases mildly with the density ratio . The density ratio and the Atwood number play an important role in the perturbation of the fluid interface. In our work, the mathematical model is considered in the parametric setting by regarding the Atwood number as a problem parameter which lies in a parametric domain . In this way, the physical quantities are treated as functions not only of the time and the Eulerian coordinates of the particle , but also the problem parameter .

### 2.2 Finite element modeling

In [dobrev2012high], a general framework of high-order curvilinear finite elements and adaptive time stepping of explicit time integrators is proposed for numerical discretization of the Lagrangian hydrodynamics problem over general unstructured two-dimensional and three-dimensional computational domains. Using general high-order polynomial basis functions for approximating the state variables, and curvilinear meshes for capturing the geometry of the flow and maintaining robustness with respect to mesh motion, the method achieves high-order accuracy. On the other hand, a modification is made to the second-order Runge-Kutta method to compensate for the lack of total energy conservation in standard high-order time integration techniques. The introduction of an artificial viscosity tensor further generates the appropriate entropy and ensures the Rankine-Hugoniot jump conditions at shock boundaries.

Following [dobrev2012high], we adopt a spatial discretization for (2.1) using a kinematic space for approximating the position and the velocity, and a thermodynamic space for approximating the energy. The density can be eliminated, and the equation of mass conservation can be decoupled from (2.1). We assume high-order finite element (FE) discretization in space, and the finite dimensions and are the global numbers of FE degrees of freedom in the corresponding discrete FE spaces. For more details, see [dobrev2012high]. The FE coefficient vector functions for velocity and position are denoted as , and the coefficient vector function for energy is denoted as . The semidiscrete Lagrangian conservation laws can be expressed as a nonlinear system of differential equations in the coefficients with respect to the bases for the kinematic and thermodynamic spaces:

(2.5) | ||||||

where denotes the effects of the gravitation force in the discrete system.

Let , , be the hydrodynamic state vector. Then the semidiscrete conservation equation of (2.5) can be written in a compact form as

(2.6) |

where the nonlinear force term, , is defined as

(2.7) |

where and are nonlinear vector functions that are defined respectively as

(2.8) |

From now on, we drop the dependence on the parameter to simplify the notations where there is no ambiguity.

In order to obtain a fully discretized system of equations, one needs to apply a time integrator. We consider an explicit Runge-Kutta scheme called the RK2-average scheme, which is proved to conserve the discrete total energy exactly (see Proposition 7.1 of [dobrev2012high]). The temporal domain is discretized as , where

denotes a discrete moment in time with

, , and for , where . The computational domain at time is denoted as . We denote the quantities of interest defined on with a subscript . Starting with , , and , the discrete state is updated by(2.9) | ||||||

(2.10) | ||||||

(2.11) |

where the state is used to compute the updates

(2.12) |

in the first stage. Similarly, is used to compute the updates

(2.13) |

with in the second stage. Note that the RK2-average scheme is different from the midpoint RK2 scheme in the updates for energy and position. The RK2-average scheme uses the midpoint velocity and the average velocity to update energy and position in the first stage and the second stage respectively, while the midpoint RK2 uses the initial velocity and the midpoint velocity . Since an explicit Runge-Kutta method is used, we need to control the time step size in order to maintain the stability of the fully discrete scheme. We follow the automatic time step control algorithm described in Section 7.3 of [dobrev2012high]

, where the time step size is controlled by estimates at all quadrature points in the mesh used in the evaluation of the force matrix

.## 3 Reduced order model

In this section, we present the details of the projection-based reduced oder model for the semi-discrete Lagrangian conservation laws (2.5). The reduced order model is constructed in the offline phase and deployed in the online phase. In what follows, we first discuss all the essential ingredients of the reduced order model, and move on to discuss the construction.

### 3.1 ROM simulation

In the reduced order model, we restrict our solution space to a subspace spanned by a reduced basis for each field. That is, the subspaces for velocity, energy, and position fields are defined as

(3.1) |

with , , and . Using these subspaces, each discrete field is approximated in trial subspaces, , , and . For a generic problem parameter , we write

(3.2) | ||||

(3.3) | ||||

(3.4) |

where , , and denote the prescribed offset vectors for velocity, energy, and position fields respectively; the orthonormal basis matrices , , and are defined as

(3.5) |

and , , and denote the time-dependent generalized coordinates for velocity, energy, and position fields, respectively. One natural choice of the offset vectors is to use the initial values, i.e. , , and .

The nonlinear matrix function, , changes every time the state variables evolve. Additionally, it needs to be multiplied by the basis matrices whenever updates in the nonlinear term occur, which scales with the FOM size. Therefore, we cannot expect any speed-up without special treatment of the nonlinear terms. To overcome this issue, a hyper-reduction technique needs to be applied (cf. [chaturantabut2010nonlinear]), where and are approximated as

(3.6) |

That is, and are projected onto subspaces and , where , and , , denote the nonlinear term basis matrices, and and denote the generalized coordinates of the nonlinear terms. Now we show how the generalized coordinates, , can be determined by the following interpolation:

(3.7) |

where , , is the sampling matrix and is the

-th column vector of the identity matrix

. Note that Eq. (3.7) is an over-determined system. Thus, we solve the least-squares problem, i.e.,(3.8) |

The solution to the least-squares problem (3.8) is

(3.9) |

where the Moore–Penrose inverse of a matrix , , with full column rank is defined as . Instead of constructing the sampling matrix , for efficiency we simply store the sampling indices . More precisely, the reduced matrix can be precomputed and stored in the offline phase, and is multiplied to the sampled entries to obtain by (3.9) in the online phase. Similarly, following the same procedure to compute the generalized coordinates of the nonlinear term of the energy conservation equation, we have

(3.10) |

where we denote the sampling matrix by .

Furthermore, we apply the solution nonlinear subspace (SNS) method in [choi2020sns] to establish the subspace relations and . Using the SNS relation together with the above reduced approximations, the hyper-reduced system of (2.5) can be written as :

(3.11) | ||||

(3.12) | ||||

(3.13) |

Let , , be the reduced order hydrodynamic state vector. Then the semidiscrete hyper-reduced system (3.11) can be written in a compact form as

(3.14) |

where the nonlinear force term, , is defined as

(3.15) |

Applying the RK2-average scheme to the hyper-reduced system (3.11), the RK2-average fully discrete hyper-reduced system reads:

(3.16) | ||||||

(3.17) | ||||||

(3.18) |

where the lifted ROM approximation given by

(3.19) |

is used to compute the updates

(3.20) |

in the first stage. Similarly, is used to computed the updates

(3.21) |

with in the second stage. The lifting is computed only for the sampled degrees of freedom, avoiding full order computation. Again, the time step size is determined adaptively using the automatic time step control algorithm, with the state replaced by the lifted ROM approximation , and the time step size controlled by estimates at all quadrature points used in the evaluation of the force matrices in the hyper-reduction sample mesh.

### 3.2 Solution bases construction

Now we describe how to obtain the reduced basis matrices , , and for the solution variables. It suffices to describe how to construct the reduced basis for the energy field only, i.e., , because other bases will be constructed in the same way. Proper orthogonal decomposition (POD) is commonly used to construct a reduced basis. POD [berkooz1993proper] obtains

from a truncated singular value decomposition (SVD) approximation to a FOM solution snapshot matrix. It is related to principal component analysis in statistical analysis

[hotelling1933analysis] and Karhunen–Loève expansion [loeve1955] in stochastic analysis. In order to collect solution data for performing POD, we run FOM simulations on a set of problem parameters, namely . For , let and be the final time and the number of time steps in the training FOM simulation with the problem parameter . By choosing , a solution snapshot matrix is formed by assembling all the FOM solution data including the intermediate Runge-Kutta stages, i.e.(3.22) |

where is the energy state at th time step with problem parameter for computed from the FOM simulation, e.g. the fully discrete RK2-average scheme (2.9), and with being the number of Runge-Kutta stages in a time step. Then, POD computes its thin SVD:

(3.23) |

where and are orthogonal matrices, and

is the diagonal singular value matrix. Then POD chooses the leading

columns of to set , where is -th column vector of . The basis size, , is determined by the energy criteria, i.e., we find the minimum such that the following condition is satisfied:(3.24) |

where is a -th largest singular
value in the singular matrix, , and denotes a threshold.
^{1}^{1}1We use the default value unless stated otherwise.
The POD basis minimizes over all with orthonormal columns, where denotes the
Frobenius norm of a matrix , defined as
with being the element of . Since
the objective function does not change if is post-multiplied by
an arbitrary orthogonal matrix, the POD
procedure seeks the optimal -dimensional subspace that
captures the snapshots in the least-squares sense. For more details on POD, we
refer to [hinze2005proper, kunisch2002galerkin].
The same procedure can be used to construct the other solution bases and .
Finally, the reduced bases for the nonlinear terms are
constructed by the subspace relations and .

### 3.3 Sampling indices selection

It remains to describe how to obtain the sampling matrices, i.e. and . The discrete empirical interpolation method (DEIM) is a popular choice for nonlinear model reduction. It suffices to describe how to construct the sampling matrix for the momentum nonlinear term only, i.e., , as the other matrix will be constructed in the same way. The sampling matrix is characterized by the sampling indices , which can be found either by a row pivoted LU decomposition [chaturantabut2010nonlinear] or the strong column pivoted rank-revealing QR (sRRQR) decomposition [drmac2016new, drmac2018discrete]. Algorithm 1 of [chaturantabut2010nonlinear] uses the greedy algorithm to sequentially seek additional interpolating indices corresponding to the entry with the largest magnitude of the residual of projecting an active POD basis vector onto the preceding basis vectors at the preceding interpolating indices. The number of interpolating indices returned is the same as the number of basis vectors, i.e. . Algorithm 3 of [carlberg2013gnat] and Algorithm 5 of [carlberg2011efficient] use the greedy procedure to minimize the error in the gappy reconstruction of the POD basis vectors . These algorithms allow over-sampling, i.e. , and can be regarded as extensions of Algorithm 1 of [chaturantabut2010nonlinear]. Instead of using the greedy algorithm, Q-DEIM is introduced in [drmac2016new] as a new framework for constructing the DEIM projection operator via the QR factorization with column pivoting. Depending on the algorithm for selecting the sampling indices, the DEIM projection error bound is determined. For example, the row pivoted LU decomposition in [chaturantabut2010nonlinear] results in the following error bound:

(3.25) |

where denotes norm of a vector, and is the condition number of , bounded by

(3.26) |

On the other hand, the sRRQR factorization in [drmac2018discrete] reveals a tighter bound than (3.26):

(3.27) |

where is a tuning parameter in the sRRQR factorization (i.e., in Algorithm 4 of [gu1996efficient]). If SNS is used, the nonlinear term basis is not necessarily orthogonal, but we can still obtain similar error bounds. For example, when the sRRQR factorization is used, (3.25) becomes

(3.28) |

with satisfying (3.27) (c.f. Theorem 5.2 of [choi2020sns]). Once the sampling indices are determined, the reduced matrix can be precomputed and stored.

## 4 Decomposition of solution manifold

Although the reduced order model reported in Section 3 can serve as a powerful technique in short-time simulation of Rayleigh–Taylor instability, it is not applicable to long-time simulation. The nonlinear evolution leads to the formation of mushroom-shaped vortices, which is the main interest of numerical studies in the hydrodynamic instability. The main difficulty is that, in order to capture the formation of the vortex, the time step is adaptively reduced. In some numerical examples presented, the full order model simulation ends up with to time steps. This gives rise to a huge number of snapshot samples, which imposes a heavy burden in storage and computational cost for the SVD computations. On the other hand, the penetration of the fluid interface grows with time and characterizes the solution dynamics, i.e., advection-dominated. As a result, the linear dependence among the snapshots is weak, and therefore there is no intrinsic low-dimensional subspace which can universally approximate the solution manifold, comprised of all the solutions over the parameter domain and temporal domain. Therefore, it is impossible for the reduced order model approach to achieve any meaningful speed-up and good accuracy.

To this end, we propose a framework to overcome these difficulties by employing multiple reduced order models. In the offline phase, we construct each of these reduced order models from a small subset of the snapshot samples, to ensure low dimension. In the online phase, each of these reduced order models will be used in a subdomain of time and parameter where it is supposed to provide good approximation. This framework involves a decomposition of the solution manifold and relies on an indicator which will be used to classify the snapshot samples and assign the reduced-order models. The idea is to decompose the solution manifold into submanifolds where the Kolmogorov -width decaying fast with respect to the subspace dimension, within which we can collect snapshots with strong linear dependence. This enables us to build accurate multiple low-dimensional subspaces.

### 4.1 General framework

We describe the general framework of the decomposition of the solution manifold, from which we will derive two practical examples later in this section. Let be an indicator which maps the triplet to a real value. The range of the indicator will be partitioned into subintervals, i.e.

(4.1) |

This partition can be either prescribed or determined by the snapshot data collected in the offline phase. For simplicity, we assume that , and is increasing with time for any initial state and any parameter . These assumptions are valid in the concrete examples we present in this paper. In general, one can relax these assumptions.

With the partition of the indicator range, instead of directly assembling all the snapshot samples into a single huge snapshot matrix as in (3.22), the FOM states will be first classified into groups. Let be a group index. We denote the subset of paired indices of time and offline parameter whose corresponding snapshot belongs to the th group as

(4.2) |

Then the snapshot matrix of the energy variable in the th group will be formed by assembling the corresponding snapshots, i.e.

(4.3) |

and POD is used to construct the energy solution basis as described in Section 3.2. The same procedure can be used to construct the other solution bases and .

For a generic problem parameter , let be the final time in the ROM simulation with the problem parameter , where the temporal domain is discretized as . We remark that can be different from the final time used in the snapshot sampling, and that even with the same problem setting and same final time, it is very likely that the temporal discretization used in the hyper-reduced system is different from the full order model, since the temporal discretization is adaptively controlled by the states. The computation in the online phase will be performed using different reduced bases in subintervals of the temporal domain , i.e.

(4.4) |

with being the partition points of the temporal domain. We remark that is the number of subintervals that the temporal domain is decomposed into for the parameter , and does not exceed the number of groups . In general, the temporal domain partition is parameter-dependent, i.e. the number of subintervals and the end-point of each subinterval depend on the problem parameter , and they are assigned by the indicator iteratively in the time marching of the ROM simulation. More precisely, starting with , for