On-the-fly Reduced Order Modeling of Passive and Reactive Species via Time-Dependent Manifolds
OOn-the-fly Reduced Order Modeling of Passive and Reactive Species viaTime-Dependent Manifolds
Donya Ramezanian , Arash G. Nouri , and Hessam Babaee ∗ Department of Mechanical Engineering and Materials Science, University of Pittsburgh
Abstract
One of the principal barriers in developing accurate and tractable predictive models in turbulent flowswith a large number of species is to track every species by solving a separate transport equation, whichcan be computationally impracticable. In this paper, we present an on-the-fly reduced order modelingof reactive as well as passive transport equations to reduce the computational cost. The presented ap-proach seeks a low-rank decomposition of the species to three time-dependent components: (i) a set oforthonormal spatial modes, (ii) a low-rank factorization of the instantaneous species correlation matrix,and (iii) a set of orthonormal species modes, which represent a low-dimensional time-dependent manifold .Our approach bypasses the need to solve the full-dimensional species to generate high-fidelity data — asit is commonly performed in data-driven dimension reduction techniques such as the principle componentanalysis. Instead, the low-rank components are directly extracted from the species transport equation.The evolution equations for the three components are obtained from optimality conditions of a variationalprinciple. The time-dependence of the three components enables an on-the-fly adaptation of the low-rankdecomposition to transient changes in the species. Several demonstration cases of reduced order modelingof passive and reactive transport equations are presented.
Tractable numerical simulation of passive and reactive species transport in turbulent flows has been thesubject of intense research in the past few decades [1] due to its importance in a diverse range of engi-neering and scientific applications. See for example [2–9]. High fidelity numerical simulations using directnumerical/large-eddy simulation (DNS/LES) of turbulent reactive flows with a large number of species is costprohibitive [10]. Detailed kinetic models usually contain hundreds of species and thousands of reactions [11–14]. Each species requires solving a transport equation, i.e., a partial differential equation (PDE). The cost ofsolving a large number of transport PDEs on a DNS/LES grid can be prohibitive [15]. Moreover, the memoryand input/output (I/O) cost of storing a large number of species imposes severe limitation on the numberof species that can be simulated. For example, currently I/O limitations allow for the storage of the speciesat every 400 simulation time steps in a DNS of turbulent combustion, while intermittent phenomena such asignition kernel occurs in the order of 10 simulation time steps [16]. The cost of I/O and storage continuesto grow in comparison to the cost of floating point operations in the future high performance computingarchitectures and this trend will impose even more stringent I/O and memory restrictions [17].To reduce the computational cost of large scale complex turbulent reactive simulations, a vast array oftechniques have been developed with the aim of reducing the number of species and/or reactions. Skele-tal reduction, for example [1], finds an optimized subset of a detailed model by eliminating unimportant ∗ Corresponding author. Email:[email protected]. a r X i v : . [ c s . C E ] J a n pecies and reactions. The skeletal reduction is often performed for zero-dimensional reactors and basedon different criteria, e.g. sensitivity analysis [18–20], reaction flux analysis [21–23], etc. [24, 25]. Otherreduction techniques include lumping methods [26, 27] and time-scale analysis techniques e.g. quasi steadystate approximation [28], partial equilibrium approximation [29], rate controlled constrained equilibrium [30],computational singular perturbation [31, 32] and intrinsic low dimensional manifold [33].Reduction schemes based on principle component analysis (PCA) [34–37] is another strategy to enablerealistic high fidelity simulations of turbulent flows with many species. The PCA-based reductions belongto data-driven reduction techniques, in which a reduced composition space is obtained by performing PCAon a DNS/LES dataset and then solve transport equations for the reduced principal components. In thePCA workflow, one must at least solve one full-dimensional DNS/LES a priori and store the data. Thisstep is done offline and one of its implications is that the PCA-based reductions are fined-tuned for a targetproblem and in general there it is difficult to guarantee that the reduced composition space is accurate underif operating conditions such as boundary conditions, geometry, Reynolds number and Mach number change.This motivates for an on-the-fly reduction scheme.Recently, new model-driven dimension reduction techniques have been introduced. in which the low-rank structure is extracted from the model. Dynamically orthogonal (DO) [38], bi-orthogonal (BO) [39] anddynamically bi-orthonormal decompositions (DBO) [40] are model-driven stochastic reduced-order modelingtechniques. For linear deterministic systems, optimally time-dependent (OTD) reduction was introduced[41]. In all of these techniques, closed-form (partial) differential equations for the evolution of the low-rankstructures (modes) are derived. Leveraging the mathematical elegance of these techniques, for some cases,numerical and dynamical system analyses are used to shed light into their performance. See, for example,[42] for the equivalence of DO and BO; [40] for the equivalence of DBO with BO/DO; [43, 44] for a theoreticalerror bounds for the approximation error; [45] for the exponential convergence of r OTD modes to the r mostdominant Lyapunov vectors of a dynamical system; [45] for a variational principle for deriving data-drivenDO evolution equations.In this paper, we present a new DBO formulation for on-the-fly low-rank decomposition of species trans-port equation. In particular, we present a novel variational principle for the extraction of the componentsof the low-rank decomposition. The reduction is achieved by extracting instantaneous correlations betweendifferent species on the fly and directly from the species transport equation without the need to generate data.We demonstrate the performance of the presented method for three cases. The remainder of this paper isorganized as follows: In §
2, we present a variational principle for the determination of the DBO components.In §
3, we illustrate the performance of the presented method by solving passive species transport equationswith one thousand species. We also apply the presented method to solve incompressible and compressiblereactive flows. Finally, a brief summary of the present work is presented in § Our target problem is the passive/reactive transport equation given by: ∂ Φ ∂t + ( v · ∇ )Φ = ∇ · ( ∇ Φ α ) + S (Φ , ρ, T ) , (1)and augmented with appropriate initial and boundary conditions. Here, x ∈ D ⊂ R d is the spatial coordinate,where d = 1 , t is time. In Eq. (1), v is the velocity vector,Φ = [ φ ( x, t ) , φ ( x, t ) , ..., φ n s ( x, t )] is the species concentration, and S (Φ , ρ, T ) is the source term, where ρ is the density and T is the temperature. For clarity in the exposition, we use the quasimatrix notationas introduced in [46], in which one of the matrix dimensions, representing the continuous space, is shown2ith ∞ and the other dimension is defined with an integer. Following this notation Φ( x, t ) ∈ R ∞× n s and S (Φ , ρ, T ) ∈ R ∞× n s . Moreover, α = diag { /ReSc , /ReSc , . . . , /ReSc n s } is a diagonal matrix of size n s × n s at each physical point x , where, Re is the Reynolds number, and Sc i is the Schmidt number ofthe i th species. We derive the equations for the generic case where α = α ( T ( x, t )) is space-time variabledue to its dependence to temperature. The term ( v · ∇ )Φ is interpreted as a quasimatrix of size ∞ × n s :( v · ∇ )Φ = [( v · ∇ ) φ , ( v · ∇ ) φ , . . . , ( v · ∇ ) φ n s ]. In case of S (Φ , ρ, T ) = 0, Eq. (1) represents the passivetransport equation, and for incompressible flow the source term is only a function of species, i.e. S = S (Φ).We define an inner product in the spatial domain between two fields u ( x ) and v ( x ) as: (cid:10) u ( x ) , v ( x ) (cid:11) = (cid:90) D u ( x ) v ( x ) dx, (2)and the L norm induced by this inner product as: (cid:13)(cid:13) u ( x ) (cid:13)(cid:13) = (cid:10) u ( x ) , u ( x ) (cid:11) . (3)The Frobenius norm of a quasimatrix A ( x ) = [ a ( x ) , a ( x ) , . . . , a n ( x )] ∈ R ∞× n is defined as: (cid:13)(cid:13) A ( x ) (cid:13)(cid:13) F = n (cid:88) i =1 (cid:90) D a i ( x ) dx. (4)We also define the inner product between quasimatrices U ( x, t ) ∈ R ∞× m and v ( x, t ) ∈ R ∞× n as S ( t ) = (cid:10) U ( x, t ) , V ( x, t ) (cid:11) , where S ( t ) ∈ R m × n is a matrix with components S ij ( t ) = (cid:10) u i ( x, t ) , v j ( x, t ) (cid:11) , where v j ( x, t ) is the j th columnof V ( x, t ). Following the above definition, we observe that (cid:10) U R U , V (cid:11) = R TU (cid:10) U, V (cid:11) and (cid:10)
U, V R V (cid:11) = (cid:10) U, V (cid:11) R V for any R U ∈ R m × m and R V ∈ R n × n . Our goal is to solve for a low-rank decomposition of Φ( x, t ) instead of solving Eq. (1). To this end, weconsider the DBO decomposition for the full species as follows:Φ( x, t ) = r (cid:88) j =1 r (cid:88) i =1 u i ( x, t )Σ ij ( t ) y Tj ( t ) + E ( x, t ) , (5)where U ( x, t ) = [ u ( x, t ) | u ( x, t ) | ... | u r ( x, t )] , (6a) Y ( t ) = [ y ( t ) | y ( t ) | ... | y r ( t )] . (6b)Here, r < n s is the reduction size, U ( x, t ) ∈ R ∞× r is a quasimatrix whose columns are orthonormal spatialmodes, Σ( t ) ∈ R r × r is a reduced factorization of the correlation matrix, Y ( t ) ∈ R n s × r is the matrix oforthonormal species modes and E ( x, t ) ∈ R ∞× n s is the reduction error (Fig. 1). The key observation aboutthe above decomposition, given by Eq. (5) is that all three components are time dependent. As we explain,the time dependence of DBO components enables the decomposition to adapt on the fly to changes in Φ( x, t ).Orthonormality of spatial modes and species modes implies that: (cid:10) u i ( x, t ) , u j ( x, t ) (cid:11) = δ ij , (7a) y Ti ( t ) y j ( t ) = δ ij . (7b)3 … ≈ ∞ × n s ∞ × r r × r r × n s ( x, t ) U ( x, t ) ⌃( t ) Y ( t ) T Figure 1: Schematic of dynamically bi-orthonormal decomposition of species. The species quasimatrix Φ( x, t ) is decomposed to U ( x, t )Σ( t ) Y ( t ) T , where U ( x, t ) represents a set of orthonormal spatial modes, Σ( t ) is a factorization of the low-rank correlation matrix and Y ( t ) is the matrix orthonormal species modes and it represents a low-dimensional time-dependent manifold. In the above decomposition, n s represents the number of species, r is the rank of DBO and in the discrete representation ∞ will be replaced with the total number ofgrid points. The schematic of the DBO decomposition is shown in Figure 1. In the cases, where Φ( x, t ) representsmass fraction, n s − (cid:80) n s i =1 φ i ( x, t ) = 1. Moreover, it is possible to include temperature and pressure in DBO in a straightfor-ward manner. However, to maintain the generality of the presented approach for passive transport andincompressible reactive flow, we only focus on the transport of n s species.There is a duality in the DBO decomposition, in which U is a low-rank basis in the physical domainand analogously Y constitutes a low-rank basis in the species space. As we demonstrate, DBO closelyapproximates the dominant instantaneous correlated structures in the species. To formalize this connection,we define the species two-point correlation operator and its approximation by DBO as follows: C ( x, x (cid:48) , t ) = Φ( x, t )Φ( x (cid:48) , t ) T (cid:39) U ( x, t )Σ( t ) Y ( t ) T Y ( t )Σ( t ) T U ( x (cid:48) , t ) T = U ( x, t ) C ( t ) U ( x (cid:48) , t ) T , (8)where we have used the orthonormality condition Y ( t ) T Y ( t ) = I . In the above equation, C ( x, x (cid:48) , t ) ∈ R ∞×∞ is the species two-point correlation operator. In the discrete representation, C ( x, x (cid:48) , t ) is a matrix of size N by N , where N is the total number of grid points, and C ( t ) = Σ( t )Σ( t ) T ∈ R r × r is the low-rank correlationmatrix. This shows that Σ( t ) is a factorization of the symmetric positive matrix C ( t ). As we demonstrate,the eigenvalues of C ( t ) approximate the r largest eigenvalues of C ( x, x (cid:48) , t ). In this section we present a variational principle whose optimality conditions lead to closed form evolutionequations for the components of the DBO decomposition. The variational principle is given by: F ( ˙ U ( x, t ) , ˙Σ( t ) , ˙ Y ( t )) = (cid:13)(cid:13)(cid:13)(cid:13) ∂ ( U ( x, t )Σ( t ) Y ( t ) T ) ∂t − M (Φ) (cid:13)(cid:13)(cid:13)(cid:13) F , (9)subject to the orthonormality conditions of u i ( x, t ) and y i ( t ) modes given by Eqs. (7a)-(7b). In Eq. (9),˙( ) = ∂ ( ) /∂t and M (Φ) is the right hand side of species transport equation: M (Φ) = − ( v · ∇ )Φ + ∇ · ( ∇ Φ α ) + S (Φ , ρ, T ). The variational principle given in Eq. (9) seeks to minimize the difference between theright hand side of the species transport equation and the time derivative of the DBO decomposition subjectto the orthonormality constraints given by Eqs. (7a)-(7b). The control parameters are ˙ U , ˙ Y and ˙Σ. Theorthonormality constraints can be incorporated into the variational principle via Lagrange multipliers. To4his end, we first take a time derivative of the orthonormality constraints. This results in: (cid:10) ˙ u i ( x, t ) , u j ( x, t ) (cid:11) + (cid:10) u i ( x, t ) , ˙ u j ( x, t ) (cid:11) = 0 , (10)˙ y Ti ( t ) y j ( t ) + y Ti ( t ) ˙ y j ( t ) = 0 . (11)We denote ϕ ij ( t ) = (cid:10) u i ( x, t ) , ˙ u j ( x, t ) (cid:11) and θ ij ( t ) = y Ti ( t ) ˙ y j ( t ). From Eqs. (10) and (11), it is apparent thatboth ϕ ( t ) ∈ R r × r and θ ( t ) ∈ R r × r are skew-symmetric matrices, i.e. ϕ T ( t ) = − ϕ ( t ) and θ T ( t ) = − θ ( t ).By incorporating the orthonormality constraints into the variational principle via Lagrange multipliers, thefollowing unconstrained optimization problem is obtained: G ( ˙ U , ˙Σ , ˙ Y , λ, γ ) = (cid:13)(cid:13)(cid:13)(cid:13) ∂ ( U ( x, t )Σ( t ) Y ( t ) T ) ∂t − M (Φ) (cid:13)(cid:13)(cid:13)(cid:13) F + λ ij ( (cid:10) u i , ˙ u j (cid:11) − ϕ ij ) + γ ij ( y Ti ˙ y j − θ ij ) , (12)where λ ij and γ ij , i, j = 1 , ..., r are the Lagrange multipliers. In simple words, the variational principleseeks to optimally update the DBO components given the change in the right hand side of the speciestransport equation while preserving the orthonormality constraints. In Appendix A, derivations of closedform evolution equations for U ( x, t ), Σ( t ) and Y ( t ) are provided from the first-order optimality conditionsof the functional (12). The closed form evolution equations of U ( x, t ), Σ( t ) and Y ( t ) are as follows: ∂U∂t = ( M (Φ) Y − U (cid:10) U, M (Φ) Y (cid:11) )Σ − + U ϕ, (13) d Σ dt = (cid:10) U, M (Φ) Y (cid:11) − ϕ Σ − Σ θ, (14) dYdt = ( I − Y Y T ) (cid:10) M (Φ) , U (cid:11) Σ − T + Y θ. (15)As we show in Appendix B, any skew-symmetric choice for matrices θ and ϕ leads to an equivalent DBOdecomposition. In this work, we consider the simplest choice of ϕ = θ = 0. This choice lends itself to a simpleinterpretation: the updates of the low-rank subspaces ( ˙ U and ˙ Y ) are orthogonal to the current subspace, i.e. , (cid:10) U, ˙ U (cid:11) = 0 and Y T ˙ Y = 0. This choice is referred to as dynamically orthogonal condition and it wasalso used for the evolution of time dependent basis for the application of computing sensitivities [41, 47] orstochastic reduced order modeling [38, 40, 48]. We now replace M (Φ) with the right hand side of the speciestransport equation (Eq. (1)) and use the DBO low-rank decomposition for species, i.e. Φ ∼ = U Σ Y T : M (Φ) = − ( v · ∇ ) U Σ Y T + ∇ · ( ∇ U Σ Y T α ) + S ( U Σ Y T , ρ, T ) . (16)The projection of M (Φ) onto the space spanned by the species modes is obtained by multiplying the aboveequation from right by matrix Y . This results in: M (Φ) Y = − ( v · ∇ ) U Σ + ∇ · ( ∇ U Σ α Y ) + S ( U Σ Y T , ρ, T ) Y, (17)where M (Φ) Y ∈ R ∞× r , and α Y = Y T α Y ∈ R r × r is the low rank diffusion matrix and for cases where α is afunction of x and t , α Y is computed at each physical point x and time t . Similarly, the projection of M (Φ)onto the spatial modes is obtained by taking the inner product of M (Φ) with U from right: (cid:10) M (Φ) , U (cid:11) = − Y Σ T (cid:10) ( v · ∇ ) U, U (cid:11) + (cid:10) ∇ · ( ∇ U Σ Y T α ) , U (cid:11) + (cid:10) S ( U Σ Y T , ρ, T ) , U (cid:11) , (18)5here we have used (cid:10) ( v · ∇ ) U Σ Y T , U (cid:11) = Y Σ T (cid:10) ( v · ∇ ) U, U (cid:11) . Substituting Eqs. (17) and (18) into Eqs. (13-15)yields:Evolution of spatial modes (PDE): ∂U∂t = (cid:89) ⊥ U (cid:2) − ( v · ∇ ) U + ∇ · ( ∇ U Σ α Y )Σ − + SY Σ − (cid:3) , (19)Reduced order model (ODE): d Σ dt = − (cid:10) U, ( v · ∇ ) U (cid:11) Σ + (cid:10) U, ∇ · ( ∇ U Σ α Y ) (cid:11) + (cid:10) U, SY (cid:11) , (20)Evolution of species modes (ODE): dYdt = (cid:89) ⊥ Y (cid:2)(cid:10) ∇ · ( ∇ U Σ Y T α , U (cid:11) + (cid:10) S, U (cid:11)(cid:3) Σ − T , (21)where for f ( x ) ∈ R ∞× and z ∈ R n s × : (cid:89) ⊥ U f = f − U (cid:10) U, f (cid:11) and (cid:89) ⊥ Y z = z − Y Y T z, are the orthogonal projections onto the complement of U and Y , respectively. In the above equations S = S ( U Σ Y T , ρ, T ). We make the following observations about Eqs. (19-21):1. The DBO decomposition has an on-the-fly built-in adaptivity to “chase” the low-rank subspace thatthe species belong to. This can be realized by noting that the right hand side of the evolution equationof U , for example, is equal to the orthogonal projection of R U = M (Φ) Y Σ − onto the complement of U . This means that if R U is in the span of U , then R U can be exactly expressed as a linear combinationof U : R U = U C U for some C U ∈ R r × r . In that case, ˙ U = U C U − U (cid:10) U, U C U (cid:11) = U C U − U (cid:10) U, U (cid:11) C U = 0,where we have used (cid:10) U, U (cid:11) = I . However, when R U is not in the span of U , the basis optimally evolvesto follow the right hand side. An analogous mechanism exists for the evolution of Y .2. The contribution of the convective terms to the evolution of Y is zero. This is because the convectiveterm (cid:10) ( v · ∇ ) U Σ Y T , U (cid:11) = Y C Y for C Y = Σ T (cid:10) ( v · ∇ ) U, U (cid:11) ∈ R r × r is in the span of Y . Therefore, theprojection of Y C Y onto the complement of Y is zero. Similarly, for the special case of space independentand equal diffusivity coefficient, i.e. α = α = · · · = α n s , the contribution of diffusive term to ˙ Y isalso zero.3. Eq. (20) is a reduced order model (ROM) obtained by the projection of the full-dimensional speciestransport equation onto spatial and species modes.The approximation error is computed as the difference between the DBO reconstruction versus the full-rank solution of the species. We use a relative error to measure the accuracy of DBO decomposition. Therelative error is defined as: E ( t ) = (cid:13)(cid:13) Φ( x, t ) − U ( x, t )Σ( t ) Y ( t ) T (cid:13)(cid:13) F (cid:14)(cid:13)(cid:13) Φ( x, t ) (cid:13)(cid:13) F , (22)where Φ( x, t ) represent the solution of the full-dimensional species transport equation.There are two main distinctions between the stochastic DBO, which was recently introduced in [40], andEq. (5) that we highlight here: (i) Eq. (5) approximates deterministic fields (i.e., multiple species), whereasthe stochastic DBO decomposes a random field, in which the y i components are infinite-dimensional randomprocesses. (ii) In the stochastic DBO a Reynolds decomposition of the stochastic field is considered, i.e. the stochastic field is decomposed to the mean and fluctuations and the stochastic DBO decomposition isconsidered only for the approximation of the fluctuations. On the other hand, Eq. (5) does not consider aseparate mean field. This simplifies the DBO evolution equations significantly. Moreover, in this work, itis shown that the evolution equations for the DBO components minimize the variational principle given byEq. (12). The variational principle can facilitate further adjustments to the DBO decomposition in a rigorousand unambiguous manner. 6 .4 Computational Cost The main computational advantage of using DBO is that the transport PDE is solved only for r spatialmodes (Eq. (18)) as opposed to n s species in the full transport equation. The computational cost of evolvingΣ and Y is negligible as they are governed by low-rank ordinary differential equations (ODEs). Moreover,in the DBO decomposition, the species are stored in the compressed form , i.e. , matrices U , Σ and Y arekept in the memory as opposed to their multiplication U Σ Y T , i.e., the decompressed form . The memorystorage requirement is dominated by U as Σ and Y are low-rank matrices and their storage cost is negligible.Therefore, in comparison to the full species transport equation, this results in the memory compression ratioof n s /r . As we show below, if care is taken, it is possible to evolve the DBO components and not storethe ∞ × n s species quasimatrix in the decompressed at any point. Evidently, the same compression ratiois gained when storing the time-resolved DBO solution to the disk, in which U , Σ and Y are written tothe disk and any species can be reconstructed from the DBO decomposition. The memory and I/O savingsare very important for the future high performance computing architecture, where power restrictions imposestringent limitations on memory and I/O usage [17]. In this section we discuss how each term in the righthand side of Eqs. (18-21) can be computed. For brevity, we drop the dependence of S on ρ and T . • ∇ · ( ∇ U Σ α Y ): This term can be written as ∇ · ( ∇ U Σ α Y ) = ∇ · ( ∇ ˜ U ) ∈ R ∞× r , where ˜ U = U D and thelow-rank matrix D = Σ α Y ∈ R r × r can be computed point wise. This shows that computing this termrequires computing Laplacian of r field variables. This term is then utilized in Eqs. (19) and (20). • S ( U Σ Y T ): This term can be computed in a loop over all grid points and when the species concen-tration at point x ∗ is required it can be computed as: S (Φ( x ∗ , t )) = S ( U ( x ∗ , t )Σ( t ) Y ( t ) T ), where U ( x ∗ , t )Σ( t ) Y ( t ) T ∈ R × n s is the vector of species concentration that is utilized in detailed kinetics.Again, the matrix S ( U Σ Y T ) ∈ R ∞× n s never have to be stored as only the projection of S onto Y and U , i.e. , SY ∈ R ∞× r and (cid:10) S, U (cid:11) ∈ R n s × r are needed in the DBO equations. Both (cid:10) S, U (cid:11) and SY needto be stored and they can be computed as a running summation and they can be stored in the sameloop where S is calculated. • ∇ · ( ∇ U Σ Y T α ): Similar to the computation of S , this term can be computed for each point x ∗ . Thisrequires the computation of the Laplacian of n s field variables. However, only the projection of thisterm onto U ( (cid:10) ∇ · ( ∇ U Σ Y T α ) , U (cid:11) ) appears in the DBO equations, and this term can be computed asa running summation in the same loop over grid points. The DBO spatial and species modes can be ranked based on “energy” in the second norm sense. The rankingcan be achieved by performing a singular value decomposition (SVD) of Σ( t ):Σ( t ) = R U ( t ) ˜Σ( t ) R Y ( t ) , (23)where ˜Σ( t ) is a diagonal matrix that contains the ranked singular values: ˜ σ ( t ) > ˜ σ ( t ) > ... > ˜ σ r ( t ) and R U ( t ) and R Y ( t ) are orthonormal matrices that can be used to rotate U and Y as follows:˜ U ( x, t ) = U ( x, t ) R U ( t ) , (24a)˜ Y ( t ) = Y ( t ) R Y ( t ) . (24b)The components { ˜ U ( x, t ) , ˜Σ( t ) , ˜ Y ( t ) } represent the DBO decomposition in the canonical form. We notethat the DBO in the canonical form and the form that is computed are equivalent: U Σ Y T = ˜ U ˜Σ ˜ Y T . SeeAppendix B for a more details on equivalent decompositions. In any demonstration figures in this paper, thecomponents of the DBO decomposition are shown in the canonical form.7 .6 Static vs. Time-Dependent Manifolds In this section, we draw contrast between time-dependent manifolds extracted by DBO and the static man-ifolds extracted from PCA. To this end, we write the DBO and PCA decompositions in the matrix form:PCA decomposition: Φ( x, t ) (cid:39) U P CA ( x, t ) Y TP CA , (25a)DBO decomposition: Φ( x, t ) (cid:39) U ( x, t )Σ( t ) Y ( t ) T . (25b)In the PCA decomposition, Y P CA ∈ R n s × r represents the static manifold. In PCA, the full-dimensionalspatiotemporal species data is required in the form of: A = [Φ T ( x ∗ , t ∗ )] ∈ R n s × n , where x ∗ and t ∗ areselected points and times, respectively and n is the total number of space-time points. The data on speciesΦ( x, t ) is obtained by either a high-fidelity simulation, e.g. DNS or from experiment. The matrix Y P CA consists of the first r eigenvetors of the correlation matrix: C = AA T ∈ R n s × n s associated with the r largest eigenvalues of C . An evolution equation for the principal component U P CA ( x, t ) is obtained byplugging the decomposition, given in Eq. (25a), into Eq. (1) and perform the Galerkin projection onto Y P CA . There are several key differences between DBO and PCA decompositions: (i) the PCA decompositionrelies on high-fidelity data, i.e. one must commit to preforming a high-fidelity simulation of a canonicalproblem with all species, or conducting an experiment and storing the high-dimensional data. That maybe a prohibitively expensive undertaking for problems with a very large number of species both in termsof computation/experiment and storage requirements of the space-time resolved matrix A . On the otherhand, the DBO decomposition does not require data: it extracts the low-dimensional structure directly fromthe species transport equation. In that sense, DBO eliminates the potentially expensive and offline step ofhigh-fidelity data generation that is required in the PCA workflow. We note that the PCA workflow is similarto many other common data-driven reduction techniques such as proper orthogonal decomposition (POD)[49, 50] and dynamic mode decomposition (DMD) [51–54]. (ii) Reliance of the PCA approximation to high-fidelity training data means that Y P CA may not be a good low-rank representative of the full compositionspace when used for different operating conditions, e.g. boundary conditions, geometry, Mach and Reynoldsnumbers. This potential limitation does not exist for the DBO decomposition, because DBO is solved for theproblem at hand. (iii) In the DBO decomposition, the low-rank matrix Y ( t ) is a time-dependent manifold,whereas in PCA, Y P CA is static. This allows the DBO decomposition to instantaneously adapt to changes inΦ( x, t ). More specifically, Y P CA is a low-dimensional manifold in a time-averaged sense. This can be realizedby inspecting the continuous analogue of the eigen-decomposition of the correlation matrix C that is formedin PCA: C ij = T (cid:82) T (cid:82) D φ i ( x, t ) φ j ( x, t ) dxdt . As we demonstrate in our results, the DBO decomposition closelyapproximates the eigen-decomposition of the instantaneous correlation matrix: ˆ C ij ( t ) = (cid:82) D φ i ( x, t ) φ j ( x, t ) dx .We compare the performance of DBO with the reduction based on instantaneous principal componentanalysis (I-PCA). The I-PCA components can be computed in a data-driven manner by computing SVDof the instantaneous matrix of full species: Φ( x, t ) = ˆ U ( x, t ) ˆΣ( t ) ˆ Y T ( t ), where ˆ( ∼ ) denotes the componentsof I-PCA, ˆ U ( x, t ) = { ˆ u ( x, t ) , ˆ u ( x, t ) , . . . , ˆ u n s ( x, t ) } is the quasimatrix of left singular functions, ˆΣ( t ) =diag(ˆ σ ( t ) , ˆ σ ( t ) , . . . , ˆ σ n s ( t )) is the diagonal matrix of singular values and ˆ Y ( t ) = { ˆ y ( t ) , ˆ y ( t ) , . . . , ˆ y n s ( t ) } isthe matrix of right singular vectors. It is also straightforward to show that the I-PCA spatial and speciesmodes are the eigenfunctions and eigenvectors of the two-point correlation operator (Eq. (8)) and instanta-neous correlation matrix, respectively, as shown below: (cid:90) D C ( x, x (cid:48) , t )ˆ u i ( x (cid:48) , t ) dx (cid:48) = ˆ σ i ( t )ˆ u i ( x, t ) and ˆ C ( t )ˆ y i ( t ) = ˆ σ i ( t )ˆ y i ( t ) , i = 1 , , . . . , n s . (26)A rank- r truncated I-PCA represents the best approximation that any rank r decomposition in the form ofEq. (5) can achieve at any given time t and therefore, comparison of DBO components in the canonical formwith I-PCA shows how closely DBO is approximating the optimal low-rank decomposition.8 Demonstration Cases
In this section, we demonstrate the application of DBO in solving passive transport equation with manyspecies. The passive species are governed by Eq. (1) for the case of S (Φ , ρ, T ) = 0. We consider n s = 1000species with different diffusivities and initial conditions. The solution of Eq. (1) is considered with v ( x, t )obtained from the solution of Burger’s equation with the presence of shocks. We compare the solutionobtained by the DBO reduction against the same-rank I-PCA reduction.We consider the Burgers’ equation governed by: ∂v∂t + v ∂v∂x = ν ∂ v∂x , x ∈ [0 , π ] , and t ∈ [0 , t f ] , (27)where v is the velocity of the flow and ν is the viscosity of the fluid which is assumed to be 0 .
01. The Burgers’equation is solved with the initial condition of v ( x,
0) = 0 . x ) − .
5) sin( x + 2 π × .
37) and the initialcondition of the passive species transport equation is assumed to be: φ i ( x,
0) = (cid:80) n s n =1 ζ ( n ) i n b sin( nπxL ), where ζ ( n ) i are chosen from an independent normal distribution and b is the rate of decay of the spectrum and isconsidered to be 2. The diffusivity of the i th species is considered to be α i = 0 . / √ i, i = 1 , . . . , n s and thelength of the physical domain is L = 2 π . The reason that each species has a different concentration is becauseof different initial condition and different diffusivity coefficients. The DBO spatial modes, governed by Eq.(19), are solved using spectral Fourier method with N = 512 Fourier modes. The fourth-order Runge-Kuttascheme is used for the time integration of the DBO Eqs. (19-21) with ∆ t = 1 /
256 and up to the final time of t f = 4 time units. In Figure 2(a), solutions of two different passive species ( φ and φ ) using reduction sizesof r = 2 , t = 4 are shown. The presence of the traveling shock in v manifests itself with a sharpchange in the passive species. As it can be seen, even for r = 2, i.e. , a drastic reduction, the location of theshock is predicted correctly and as reduction rank increases from r = 2 to r = 8 the passive species obtainedby DBO converge to the true profile obtained by directly solving the full-dimensional passive species and thelocation and amplitude of the shock are correctly captured. The accuracy of DBO is also compared againstthe I-PCA. As Fig. 2(b) represents, the first two I-PCA spatial modes (ˆ u ( x, t ) and ˆ u ( x, t )), which are thefirst two dominant eigenfunctions of C ( x, x (cid:48) , t ), match relatively well with the DBO spatial modes (˜ u ( x, t )and ˜ u ( x, t )) for r = 2, and as the reduction size increases to r = 4 and r = 8 the agreement between thefirst two DBO modes and those of the I-PCA improves.In Fig. 3(a) the evolution of the error using different reduction sizes are shown. The error improves aboutone order of magnitude by increasing the reduction size from r = 2 to r = 12. In addition, the evolution ofsingular values extracted from I-PCA and DBO solutions are presented in Fig. 3(b). It is observed that DBOcan capture the largest singular values with a better accuracy. However, the accuracy of DBO degrades forcapturing lower singular values. The difference between singular values of I-PCA and those of DBO is theresult of the lost interactions of DBO modes with the unresolved modes. The unresolved error affects themodes with smaller singular values more intensely. Ultimately, the DBO can be augmented with an adaptivestrategy to add/remove modes based on a criterion. Similar strategies have been adopted in the past. Seefor example [55]. For a discussion of the error in ROM based on time-dependent basis, see [47]. We investigate the performance of DBO for an incompressible reactive flow, in which velocity is not affectedby reactions and the reactive source term is a function of Φ only, i.e. , S (Φ). To this end, we use DBO to solvea biochemically reactive flow, in which the coagulation cascade process in a Newtonian fluid is modeled by9 (a) Passive species (b) Spatial modes Figure 2: DBO of passive species transported in 1D Burger’s equation. (a) Comparison of solutions of φ and φ obtained from DBOusing reduction sizes of r = 2 , t = 4. (b) Comparison ofthe first two dominant spatial modes using reduction sizes of r = 2 , t = 4 with the optimal spatial modes obtained from I-PCA. -2 -1 (a) Error (b) Singular values Figure 3: DBO of passive species transported in 1D Burger’s equation: (a) Evolution of solution error using reduction sizes of r = 2 , r = 2 , Anand et al. [56]. In this model, the coagulation cascade process is simulated by solving for n s =23 species.The velocity field is obtained by solving the incompressible Navier-Stokes equations: ∂v∂t + ( v · ∇ ) v = −∇ p + 1 Re ∇ v, ∇ · v = 0 , (28)where v = ( v x , v y ) is the velocity vector field, p is the pressure field, and Re is the Reynolds number. Theinvolved species, their corresponding Schmidt numbers, as well as source terms of each reactant can be foundin [57]. 10 x i at inlet outflow L H i = 0 Figure 4: 2D incompressible turbulent reactive flow: Schematic of the computational domain and boundary conditions.Figure 5: 2D incompressible turbulent reactive flow. The four most dominant DBO spatial modes in different times.
The schematic of the problem is shown in Figure 4. Height and length of the jet are H = 2 and L = 10, correspondingly and the Reynolds number based on reference length of half the height ( H/ ν is Re = wH/ ν = 1000. The species boundary condition at the inlet is φ i (0 , x , t ) = 1 / (cid:0) tanh ( x + H/ /δ − tanh ( x − H/ /δ (cid:1) for all species, where δ = 0 .
1. For the spatialdiscretization, spectral/hp element method is used with N e = 4008 quadrilateral elements and polynomialorder of 5. For more details of the spectral/hp element method see [58–60]. The fourth-order Runge-Kuttascheme is used for the time integration with ∆ t = 5 × − and the system is solved till t f = 5 units oftime. The velocity field at each time step is then interpolated on a Cartesian spectral element mesh in therectangular domain shown by dashed line in Figure 4. The DBO equations are solved in the rectangulardomain with 251 elements in the horizontal direction and 76 elements in the vertical direction. The spectralpolynomial of order 5 is used in each direction.We solve the problem using r = 5 reduction size. Figure 5 shows the first 4 dominant spatial modes at 5different instants. The first mode, associated with the largest singular value, is positive in all x and t . Thismode is the most energetic mode and it roughly approximates the orthonormalized mean of the species, i.e. ,˜ u ( x, t ) (cid:39) φ ( x, t ) / (cid:107) φ ( x, t ) (cid:107) , where φ ( x, t ) = 1 /n s (cid:80) n s i =1 φ i ( x, t ). This occurs when the mean of the species isin fact the most energetic mode. The higher modes change sign in the domain and capture finer structures.Also it is clear that the modes are advecting with the flow from left to right.11 igure 6: Time-dependent low-dimensional manifold of 2D incompressible turbulent reactive flow: dominant species modes in different timesin the form of color-coded matrices.Figure 7: 2D incompressible turbulent reactive flow: species concentration of first and tenth species in different times obtained from DNS(full-dimensional) and DBO using r = 5. Figure 6 displays the species matrix ˜ Y ( t ) in different times. Similar to the spatial modes, the first modeassociated with the most energetic direction is always positive. This mode changes slowly with time. On theother hand the higher modes have positive and negative contribution of each species and they change fasterwith time. Each component of vector ˜ y i ( t ) should be interpreted as the instantaneous contribution of thecorresponding species. The species labels are shown in Figure 6. Unlike conventional reduction schemes such12 a) Singular values -10 -8 -6 -4 -2 (b) Solution error Figure 8: 2D incompressible turbulent reactive flow. (a) Evolution of eigenvalues obtained from DBO solutions using r = 5 reduction sizeand I-PCA solutions. (b) Evolution of solution error using reduction sizes of r = 3 , as skeletal reductions, DBO does not eliminate any species or reactions. Instead, a time-dependent weightedcontributions of all species are considered in the low-rank decomposition. For example, as it is clear fromFigure 6, φ and φ have a small component in ˜ y ( t ), which vanishes with time. On the other hand, φ and φ have a larger footprint in all ˜ y i modes.Figure 7 shows several snapshots of first and tenth species concentration obtained from direct numericalsolution (DNS) of the full-dimensional species transport and DBO. We can observe that solutions obtainedfrom DBO using only r = 5 are very similar to DNS data at all times. For further examining the accuracyof DBO, singular values extracted from DBO decomposition are compared with those of the I-PCA decom-position. As presented in Figure 8(a), the evolution of singular values extracted from DBO solutions closelymatch with I-PCA, and similar to the previous demonstration case the most dominant modes show betteragreement with I-PCA. In Figure 8(b), the approximation errors of DBO for r = 3 , r . The applicability of using DBO decomposition for simulating 2D compressible reactive Navier-Stokes isdemonstrated in this section. For the mathematical description of compressible flows involving n s species,the primary transport variables are the density ρ ( x, t ), velocity vector v ( x, t ), pressure p ( x, t ), total energy E ( x, t ), temperature T ( x, t ), and species mass fractions Φ( x, t ). The reactive compressible Navier-Stokeequations are given by: 13 ρ∂t + ∂ρv j ∂x j = 0 , (29a) ∂ρv i ∂t + ∂ρv i v j ∂x j = − ∂p∂x i + ∂τ ij ∂x j , (29b) ∂ρE∂t + ∂ ( ρEv j ) ∂x j = − ∂pv j ∂x j + ∂ ( τ ij v i ) ∂x j − ∂q j ∂x j + W (Φ , ρ, T ) , (29c) ∂ Φ k ∂t + v j ∂ Φ k ∂x j = − ρ ∂J kj ∂x j + S k (Φ , ρ, T ) . (29d)Here, the viscosity flux τ , heat flux q , mass flux J , and total energy are represented by: τ ij = 1 Re (cid:16) ∂v i ∂x j + ∂v j ∂x i − ∂v k ∂x k δ ij (cid:17) , q j = − Ec.P e ∂T∂x j , J kj = − α k ∂φ k ∂x j , in dimensionless format where e is the internal energy, E = e + v i v i is the total energy, Ec = ( γ − M a , P e = Re.P r , and
M a are Eckert, Peclet, and Mach numbers, respectively. Moreover, W = (cid:80) n s k =1 ρS k ∆ h f,k is the heat release where S k is the dimensionless species source term and ∆ h f,k is the formation enthalpyof species k . We assume a perfect gas with the specific heat ratio γ = c p /c v and R = c p − c v , where R isthe gas constant, c p and c v denote the specific heats at constant pressure and constant volume, respectively,and are assumed to be constants. Equal and constant diffusion coefficients and unity Schmidt ( Sc = µ/D )and Prandtl ( P r = c p µ/λ ) numbers for all species are assumed, where D is the mass diffusivity and λ isthe heat conductivity. The viscosity and molecular diffusion coefficients can, in general, be temperaturedependent but in this study, they are assumed to be constants. As the ability of DBO decomposition inhandling different species with different diffusivities was demonstrated in previous sections, here we use asimple diffusion model. More complex simulations with realistic conditions will be the subject of futurestudies.Simulations are conducted of a 2D temporally evolving jet transport of 15 species ( n s = 15) associatedwith hydrogen-air burning based on the model in Ref. [61]. As is shown in Fig. 9(a) our temporal layerconsists of three parallel streams. The middle stream travels in opposite direction but with the same speedas the bottom and top streams. The transport variables are normalized with respect to L r , V r , ρ r , and T r where L r is the size of domain in each direction and V r = ∆ V is the velocity difference across the layer. ρ r and T r are defined for a mixture of species in which φ mH = 2 φ mO = 2 φ mN = 2 φ mAR = 2 φ mCO = 2 φ mHE = 2 /
7, under T = 300 K and p = 1 atm where Φ mk is the mole fraction of k th species. In our simulations, Re = ρ r V r L r /µ =10 and M a = V r / √ γRT r = 0 .
5. Periodic boundary condition are imposed on all four boundaries (Figure9(a)) and initial temperature and pressure on the entire domain are set to T = 4 . p = 2 . t = 3.For synchronizing hydrodynamics with combustion by equating ignition delay and vortex roll-up time, thespecies source terms are multiplied by 0 .
01. Figure 9(c) portrays the evolution of dimensionless temperatureand mass fractions of H and O in a constant pressure ( p = 2 . T = 4 and initial mass fractions shown in the fuel side of Figure 9(a). It is clear in Figure9(c) that the temperature inflection point happens at t ≈ .
5. The physical domain is discretized using theFourier spectral method with 256 Fourier modes, and we used ode113 [62] of MATLAB for time integrationof DBO components.DBO and PCA-ROM simulations are initialized at t = 1 from the DNS data and evolved with r modestill t = 3, i.e. within the orange rectangle in Fig. 9(c). PCA-ROM modes containing a linear combination of H, H , O, O , OH, H O, N , HO , H O , AR, CO, CO , OH ∗ , HE, HCO Time of temperature inflection point dT /dt | max in the fuel side of the domain igure 9: Compressible reactive flow. (a) Computational domain and boundary conditions for compressible reactive Navier-Stokes simulationwith 15 species for hydrogen-air burning based on the model in Ref. [61], (b) contour of temperature at t = 3, (c) dimensionless temperatureas well as mass fractions of H and O in a zero-dimensional constant pressure reactor initialized with T = 4 and mass fractions from thereference state. In 2D simulation the species and energy source terms are multiplies by a constant (0.01 here) to experience temperatureinflection point in the middle of the domain at the time of vortex roll-ups. n s species mass fractions were extracted from the DNS data of all the spatial points of the domain and at200 equally spaced time steps from t = 1 to t = 3. Here, the results of DBO and PCA-ROM are shown andcompared against DNS within ∆ t = 2 after their initialization ( t ∈ [0 , u i ( x, t )).Figure 11(a) shows the instantaneous singular values of DBO decomposition for r = 6 and 7 as well asthe I-PCA singular values. It is clear that the singular values of DBO very closely match with the mostdominant singular values of I-PCA. For the compressible reactive flow, the DBO approximation error affects ρ , ρv i and E due to two-way nonlinear coupling of species and these variables, since the heat source W (Φ , ρ, T )utilizes Φ (cid:39) U Σ Y T . Figure 11(b),(c) and (d) demonstrate the instantaneous errors in estimating Φ , T, and ρ , respectively. In general the error decreases with increasing r .Figure 12 demonstrates the time-plot of temperature at two fixed spatial positions in the domain Thesetwo points, labeled as P1 and P2, are shown in Figure 9(a). P1 is the center of the domain and initiallyhas the fuel composition, and P2 is very close to the shear layer on fuel side. Here, the variations in DNStemperatures after t = 0 . r <
6, whileDBO reasonably estimates the DNS temperature trends with r = 4. Figure 13 shows the differences between˜ Y modes in PCA-ROM and DBO as described in § y ( t ) barely changesin time, and it is very similar to the first PCA-ROM species mode ˜ y . This mode is associated with the15 igure 10: Compressible reactive flow. First four orthonormal DBO spatial modes in different time. products of combustion, which occupy a large portion of the domain and because of this the singular valueof the first DBO mode ˆ σ is an order of magnitude larger than the next mode, as shown in Fig. 11(a). Y DBO and Y P CA have similar patterns in t = 0 . r ≤
5. This means that the properspecies composition at each spatial point would not remain in space of the first 5 modes of PCA-ROM duringignition. The contribution of CO , for example, in the second and third species modes reveals the differencebetween PCA-ROM and DBO. In DBO, CO component in ˜ y ( t ) continuously increases from before ignition( t = 0 .
5) to after ignition ( t = 2 . CO in the secondPCA mode. Similar observations can be made for other species e.g. H O and CO . In this paper, we presented a variational principle for the determination of a low-rank decomposition (DBO)of the passive and reactive species transport equation. The optimality conditions of the variational principlelead to closed-form evolution equations for the components of the decomposition. The DBO decompositionconsists of a set of time-dependent orthonormal spatial modes, a set of orthonormal species modes and alow-rank factorization of the correlation matrix. The novelty of DBO is that all three components are timedependent – enabling the spatial and species subspaces adapt to the changes of the dynamics on the fly. TheDBO decomposition does not require the offline step of generating high-fidelity data to extract the low rankcomposition space. This step is needed in data-driven reduction techniques such as PCA. Instead, the low-rank composition space is extracted from the species transport equation. Therefore, the DBO decompositionis not fine-tuned for a target problem and it can extract the low-rank structure for the problem at hand.16 -10 -8 -6 -4 -2 (a) Singular values -8 -6 -4 -2 (b) Species error -8 -6 -4 -2 (c) Temperature error -6 -4 -2 (d) Density error Figure 11: Compressible reactive flow. (a) Evolution of eigenvalues extracted from DNS and DBO using r = 6 and 7, and (b) errors inspecies mass fractions, (c) temperature, and (d) density using r = 3 , , , We demonstrated the numerical performance of DBO for passive species transport as well as reactiveincompressible and compressible flow. In all demonstration cases, we show that the DBO decompositionclosely approximates the best instantaneous low-rank decomposition obtained by performing instantaneousPCA of the full-rank species.
Acknowledgements
This work has been supported by the National Science Foundation under Grant No. 2042918 and by NASAunder Grant No. 80NSSC18M0150. The authors thank Prof. Peyman Givi, Prof. Harsha Chelliah and Dr.Jackie Chen for useful comments that led to a number of improvements. The authors also thank PrernaPatil for her technical input. 17 (a) Temperature at P1 (b) Temperature at P2
Figure 12: Compressible reactive flow. Comparing the performance of PCA and DBO using r = 4 , , and 6 against DNS based on predictedtemperature at (a) point 1 and (b) point 2.Figure 13: Low-dimensional time-dependent manifolds for compressible reactive flow. Comparison of species modes obtained from PCA(static), I-PCA and DBO (time-dependent) in different times. A Variational Principle for DBO
For the sake of simplicity in notation, we denote u i ( x, t ) as u i , y i ( t ) as y i and Σ ij ( t ) as Σ ij . In the followingwe show that the first-order optimality condition of the variational principle leads to closed-form evolutionequations for U , Y and Σ. Throughout this derivation, we use index notation, where the repeated index18mplies summation over that index. We begin with the functional form of the variational principle: G ( ˙ U , ˙Σ , ˙ Y , λ, θ ) = (cid:10) ˙ u i , ˙ u k (cid:11) (Σ ij Σ kl )( y Tj y l )+ (cid:10) u i , u k (cid:11) ( ˙Σ ij ˙Σ kl )( y Tj y l )+ (cid:10) u i , u k (cid:11) (Σ ij Σ kl )( ˙ y Tj ˙ y l )+ 2 (cid:10) ˙ u i , u k (cid:11) (Σ ij ˙Σ kl )( y Tj y l )+ 2 (cid:10) u i , u k (cid:11) (Σ ij ˙Σ kl )( y Tj ˙ y l )+ 2 (cid:10) ˙ u i , u k (cid:11) (Σ ij Σ kl )( y Tj ˙ y l ) − (cid:10) ˙ u i , M (Φ)Σ ij y j (cid:11) − (cid:10) u i , M (Φ) ˙Σ ij y j (cid:11) − (cid:10) u i , M (Φ)Σ ij ˙ y j (cid:11) + (cid:13)(cid:13) M (Φ) (cid:13)(cid:13) F + λ ij ( (cid:10) u i , ˙ u j (cid:11) − ϕ ij ) + γ ij ( y Ti ˙ y j − θ ij ) . (30)The first-order optimality condition requires a vanishing gradient of the functional G with respect to allcontrol variables. Starting with the gradient of G with respect to ˙Σ mn . The following is obtained: ∂ G ∂ ˙Σ mn = (cid:10) u m , u k (cid:11) ˙Σ kl ( y Tn y l ) + (cid:10) u i , u m (cid:11) ˙Σ ij ( y Tj y n )+ 2 (cid:10) ˙ u i , u m (cid:11) Σ ij ( y Tj y n )+ 2 (cid:10) u i , u m (cid:11) Σ ij ( y Tj ˙ y n ) − (cid:10) u m , M (Φ) y n (cid:11) = 0 . (31)Using the orthonormality conditions and definitions of ϕ and θ , we have: (cid:10) u m , u k (cid:11) = δ mk , (cid:10) u i , u m (cid:11) = δ im and y Tn y l = δ nl , y Tj y n = δ jn , (cid:10) ˙ u i , u m (cid:11) = ϕ mi and y Tj ˙ y n = θ jn and after simplifications, the evolution equationfor Σ is obtained: d Σ mn dt = (cid:104) u m , M (Φ) y n (cid:105) − ϕ mi Σ in − Σ mi θ in . (32)Next, the gradient of G with respect to ˙ y m is set to zero. This results in: ∂ G ∂ ˙ y m = 2 (cid:10) u i , u k (cid:11) (Σ ij Σ km ) ˙ y Tj + 2 (cid:10) u i , u k (cid:11) (Σ ij ˙Σ km ) y Tj + 2 (cid:10) ˙ u i , u k (cid:11) (Σ ij Σ km ) y Tj − (cid:10) u i , Σ im M (Φ) (cid:11) + γ im y Ti = 0 . (33)To eliminate γ im , we take the inner product of the above equation and y n by multiplying y n from right: ∂ G ∂ ˙ y m y n = 2 δ ik (Σ ij Σ km ) θ nj + 2 δ ik (Σ ij ˙Σ km ) δ jn + 2 ϕ ki (Σ ij Σ km ) δ jn − (cid:10) u i , Σ im M (Φ) y n (cid:11) + γ im δ in = 0 . (34)19hen γ nm is obtained from: γ nm = 2 (cid:10) u i , Σ im M (Φ) y n (cid:11) − ij Σ im θ nj − in ˙Σ im − ϕ ki Σ in Σ km . (35)By Substituting γ nm from Eq. (35) into Eq. (33) and simplifying, We obtain:Σ ij ˙ y Tj = (cid:10) u i , M (Φ) (cid:11) ( I − y n y Tn ) + Σ ij θ nj y Tn , (36)where I ∈ R n s × n s is the identity matrix. Taking the transpose of the Eq. (36) results in:Σ ji ˙ y j = ( I − y n y Tn ) (cid:10) M (Φ) , u i (cid:11) + Σ ji θ nj y n , (37)where we have used: (cid:10) u i , M (Φ) (cid:11) T = (cid:10) M (Φ) , u i (cid:11) and θ Tnj = − θ jn . Dividing both sides of Eq. (37) by theinverse of Σ T yields the evolution equation for the species modes as in the following:˙ y j = ( I − y n y Tn ) (cid:10) M (Φ) , u i (cid:11) Σ − ji + y n θ nj . (38)Next, the gradient of G with respect to ˙ u m is considered. Since u is a function, i.e. infinite-dimensional, weuse the Fr´echet differential: G (cid:48) | ˙ U (cid:44) lim (cid:15) → G ( ˙ U + (cid:15) ˙ U (cid:48) , ˙Σ , ˙ Y , λ, γ ) − G ( ˙ U , ˙Σ , ˙ Y , λ, γ ) (cid:15) . For optimality condition the Fr´echet differential of G with respect to every column of ˙ U must vanish. Thisresults in: G (cid:48) | ˙ u m = 2 (cid:10) ˙ u (cid:48) , ˙ u k (cid:11) (Σ mj Σ kl )( y Tj y l )+ 2 (cid:10) ˙ u (cid:48) , u k (cid:11) (Σ mj ˙Σ kl )( y Tj y l )+ 2 (cid:10) ˙ u (cid:48) , u k (cid:11) (Σ mj Σ kl )( y Tj ˙ y l ) − (cid:10) ˙ u (cid:48) , Σ mj M (Φ) y j (cid:11) + λ im (cid:10) ˙ u (cid:48) , u i (cid:11) = 0 . (39)The above equation can be written as (cid:10) ˙ u (cid:48) , ∇ ˙ u m G (cid:11) and since ˙ u (cid:48) is an arbitrary perturbation, ∇ ˙ u m G mustvanish. This results in: ∇ ˙ u m G = 2 ˙ u k (Σ mj Σ kl )( y Tj y l )+ 2 u k (Σ mj ˙Σ kl )( y Tj y l )+ 2 u k (Σ mj Σ kl )( y Tj ˙ y l ) − mj M (Φ) y j + λ im u i = 0 . (40)Similar to the procedure of deriving the evolution equation for ˙ y m , to eliminate λ im an inner product of Eq.(40) and u n is taken. This results in: (cid:10) u n , ∇ ˙ u m G (cid:11) = 2 ϕ nk Σ mj Σ kl δ jl + 2 δ nk (Σ mj ˙Σ kl ) δ jl + 2 δ nk (Σ mj Σ kl ) θ jl − (cid:10) u n , Σ mj M (Φ) y j (cid:11) + λ im δ in = 0 . (41)20olving for λ nm using Eq. (41) results in: λ nm = 2 (cid:10) u n , Σ mj M (Φ) y j (cid:11) − ϕ nk Σ mj Σ kj − mj ˙Σ nj − mj Σ nl θ jl . (42)Replacing λ mn from Eq. (42) into Eq. (40) and simplifying the results yields:˙ u i Σ ij = ( M (Φ) y j − u n (cid:10) u n , M (Φ) y j (cid:11) ) + u n ϕ ni Σ ij . (43)Multiplying both sides of Eq. (43) by the inverse of Σ from right results in:˙ u i = ( M (Φ) y j − u n (cid:10) u n , M (Φ) y j (cid:11) )Σ − ij + u n ϕ ni , (44)where Σ − ij is the ij element of matrix Σ − . Therefore, evolution equations of DBO components are as in thefollowing: ∂U∂t = ( M (Φ) Y − U (cid:10) U, M (Φ) Y (cid:11) )Σ − + U ϕ, (45) d Σ dt = (cid:104) U, M (Φ) Y (cid:105) − ϕ Σ − Σ θ, (46) dYdt = ( I − Y Y T ) (cid:10) M (Φ) , U (cid:11) Σ − T + Y θ. (47)
B Equivalence of DBO Decomposition
Any choice of skew-symmetric matrices for ϕ and θ lead to equivalent decompositions. Two DBO decompo-sitions are equivalent if they represent the same low-rank subspace instantaneously. Therefore, the two DBOdecompositions { U, Σ , Y } and { ˆ U , ˆΣ , ˆ Y } are equivalent if and only if: U Σ Y T = ˆ U ˆΣ ˆ Y T . As a result, ifˆ U = U R U , ˆ Y = Y R Y and ˆΣ = R TU Σ R Y , (48)for any orthonormal matrices R U ∈ R r × r and R Y ∈ R r × r , then it is straightforward to show that the twodecompositions are equivalent using R TU R U = I and R Y R TY = I , where I ∈ R r × r is the identity matrix. Thematrices R U and R Y are in-subspace rotations . i.e. , both U and ˆ U span the same subspace. The same istrue for Y .Now, let { U, Σ , Y } and { ˆ U , ˆΣ , ˆ Y } be the solution of Eqs. (13-15) with the skew-symmetric matrices { ϕ, θ } and { ˆ ϕ, ˆ θ } , respectively. It can be shown that the two decompositions are equivalent for all t > t = 0 with the corresponding rotation matrices R U and R Y , and (ii) the rotationmatrices R U and R Y evolve according to:˙ R U = R U ˆ ϕ − ϕR U and ˙ R Y = R Y ˆ θ − θR Y , (49)with the initial conditions of R U and R Y , respectively.It was shown in Lemma 2.1 in Ref. [41] that matrices R U and R Y remain orthonormal for all t > { ˆ U , ˆΣ , ˆ Y } is an in-subspace rotation of { U, Σ , Y } according to Eq. (48) and therotation matrices are governed by Eq. (49). To this end, we start from the evolution equations for { ˆ U , ˆΣ , ˆ Y } { U, Σ , Y } .The DBO evolution equations for { ˆ U , ˆΣ , ˆ Y } are given by:˙ˆ U = (cid:89) ⊥ ˆ U M (Φ) ˆ Y ˆΣ − + ˆ U ˆ ϕ. First we note that: (cid:81) ⊥ ˆ U = (cid:81) ⊥ U , since (cid:89) f ⊥ ˆ U = f − ˆ U (cid:10) ˆ U , f (cid:11) = f − U R U (cid:10) U R U , f (cid:11) = f − U R U R TU (cid:10) U, f (cid:11) = f − U (cid:10) U, f (cid:11) , for any f ∈ R ∞× . Therefore,˙ U R U + U ˙ R U = (cid:89) ⊥ U M (Φ) Y R Y R − Y Σ − R − TU + U R U ˆ ϕ = (cid:89) ⊥ U M (Φ) Y Σ − R − TU + U R U ˆ ϕ. Rearranging the above equation and multiplying it by R TU from right results in:˙ U = (cid:89) ⊥ U M (Φ) Y Σ − + U R U ˆ ϕR TU − U ˙ R U R TU , Now replacing ˙ R U = R U ˆ ϕ − ϕR U in the above equation results in:˙ U = (cid:89) ⊥ U M (Φ) Y Σ − + U R U ˆ ϕR TU − U ( R U ˆ ϕ − ϕR U ) R TU = (cid:89) ⊥ U M (Φ) Y Σ − + U ϕ, which is the governing equation for U . The analogous procedure can be repeated for the evolution equationof Y . Replacing the equivalence relationship for the evolution of ˆΣ results in:˙ R TU Σ R Y + R TU ˙Σ R Y + R TU Σ ˙ R Y = (cid:10) U R U , M (Φ) Y R Y (cid:11) − ˆ ϕR TU Σ R Y − R TU Σ R Y ˆ θ = R TU (cid:10) U, M (Φ) Y (cid:11) R Y − ˆ ϕR TU Σ R Y − R TU Σ R Y ˆ θ. Multiplying the above equation from left by R U and from right by R TY and rearranging to obtain the evolutionequation for ˙Σ results in:˙Σ = (cid:10) U, M (Φ) (cid:11) − R U ˆ ϕR TU Σ − Σ R Y ˆ θR TY − R U ˙ R TU Σ − Σ ˙ R Y R TY . The term ˙ R TU can be obtained from Eq. (49): ˙ R TU = ϕ T R TU − R TU ϕ T = − ϕR TU + R TU ϕ , where we have used ϕ T = − ϕ . Using this relation for ˙ R TU and replacing ˙ R Y from Eq. (49) into the above equation and aftersimplification yields the evolution equation for ˙Σ, i.e. , Eq. (14). References [1] T. Lu, C. Law, Toward Accommodating Realistic Fuel Chemistry in Large-Scale Computations, Prog.Energy Combust. Sci. 35 (2) (2009) 192–215.[2] G. T. Johnson, L. J. Hunter, A numerical study of dispersion of passive scalars in city canyons, Boundary-Layer Meteorology 75 (3) (1995) 235–262. 223] Z. Warhaft, Passive scalars in turbulent flows, Annual Review of Fluid Mechanics 32 (1) (2000) 203–240.[4] M. Anand, K. Rajagopal, K. Rajagopal, A model incorporating some of the mechanical and biochemicalfactors underlying clot formation and dissolution in flowing blood, Journal of Theoretical Medicine5 (3-4) (2003) 183–218.[5] R. Antonia, P. Orlandi, Effect of schmidt number on small-scale passive scalar turbulence, Appl. Mech.Rev. 56 (6) (2003) 615–632.[6] C. Safta, C. Madnia, Autoignition and Structure of Nonpremixed CH4/H2 Flames: Detailed and Re-duced Kinetic Models, Combust. Flame 144 (1-2) (2006) 64–73.[7] X.-X. Li, C.-H. Liu, D. Y. Leung, Numerical investigation of pollutant transport characteristics insidedeep urban street canyons, Atmospheric Environment 43 (15) (2009) 2410–2418.[8] H. Wang, M. Sun, Y. Yang, N. Qin, A passive scalar-based method for numerical combustion, Interna-tional Journal of Hydrogen Energy 40 (33) (2015) 10658–10661.[9] M. Cerrolaza, S. Shefelbine, D. Garz´on-Alvarado, Numerical methods and advanced simulation in biome-chanics and biological processes, Academic Press, 2017.[10] S. B. Pope, Small scales, many species and the manifold challenges of turbulent combustion, Proceedingsof the Combustion Institute 34 (1) (2013) 1–31. doi:https://doi.org/10.1016/j.proci.2012.09.009 .URL [11] D. A. Sheen, X. You, H. Wang, T. Løv˚as, Spectral Uncertainty Quantification, Propagation and Op-timization of a Detailed Kinetic Model for Ethylene Combustion, Proc. Combust. Inst. 32 (1) (2009)535–542.[12] C. K. Law, Combustion at a crossroads: Status and prospects, Proceedings of the Combustion Institute31 (1) (2007) 1–29. doi:https://doi.org/10.1016/j.proci.2006.08.124 .URL [13] M. E. Coltrin, P. Glarborg, Chemically reacting flow: theory and practice, Wiley-Interscience, 2003.[14] A. Nouri, P. Givi, D. Livescu, Modeling and Simulation of Turbulent Nuclear Flames in Type Ia Super-novae, Prog. Aerosp. Sci. 108 (2019) 156–179.[15] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, D. Mavriplis, CFD Vision2030 Study: A Path to Revolutionary Computational Aerosciences.[16] J. C. Bennett, H. Abbasi, P. Bremer, R. Grout, A. Gyulassy, T. Jin, S. Klasky, H. Kolla, M. Parashar,V. Pascucci, P. Pebay, D. Thompson, H. Yu, F. Zhang, J. Chen, Combining in-situ and in-transitprocessing to enable extreme-scale scientific analysis, in: SC ’12: Proceedings of the InternationalConference on High Performance Computing, Networking, Storage and Analysis, 2012, pp. 1–9. doi:10.1109/SC.2012.31 .[17] Applied mathematics research for exascale computing, U.S. Department of Energy, Office of Science,Advanced Scientific Computing Research Program.[18] S. Vajda, T. Tur´anyi, Principal Component Analysis for Reducing the Edelson-Field-Noyes Model ofthe Belousov-Zhabotinskii Reaction, J. Phys. Chem. 90 (8) (1986) 1664–1670.2319] G. Esposito, H. Chelliah, Skeletal Reaction Models based on Principal Component Analysis: Applicationto Ethylene–Air Ignition, Propagation, and Extinction Phenomena, Combust. Flame 158 (3) (2011) 477–489.[20] A. Stagni, A. Frassoldati, A. Cuoci, T. Faravelli, E. Ranzi, Skeletal Mechanism Reduction ThroughSpecies-Targeted Sensitivity Analysis, Combust. Flame 163 (2016) 382–393.[21] T. Lu, C. K. Law, A directed relation graph method for mechanism reduction, Proceedings of theCombustion Institute 30 (1) (2005) 1333–1341.[22] P. Pepiot-Desjardins, H. Pitsch, An efficient error-propagation-based reduction method for large chemicalkinetic mechanisms, Combust. Flame 154 (1-2) (2008) 67–81.[23] K. E. Niemeyer, C.-J. Sung, M. P. Raju, Skeletal mechanism generation for surrogate fuels using directedrelation graph with error propagation and sensitivity analysis, Combustion and flame 157 (9) (2010)1760–1770.[24] L. Elliott, D. Ingham, A. Kyne, N. Mera, M. Pourkashanian, C. Wilson, Genetic Algorithms for Optimi-sation of Chemical Kinetics Reaction Mechanisms, Prog. Energy Combust. Sci. 30 (3) (2004) 297–328.[25] N. Sikalo, O. Hasemann, C. Schulz, A. Kempf, I. Wlokas, A Genetic Algorithm-Based Method for theAutomatic Reduction of Reaction Mechanisms, Int. J. Chem. Kinet. 46 (1) (2014) 41–59.[26] R. Djouad, B. Sportisse, N. Audiffren, Reduction of Multiphase Atmospheric Chemistry, J. Atmos.Chem. 46 (2) (2003) 131–157.[27] Y. Gao, R. Shan, S. Lyra, C. Li, H. Wang, J. Chen, T. Lu, On Lumped-Reduced Reaction Model forCombustion of Liquid Fuels, Combust. Flame 163 (2016) 437–446.[28] M. Stiefenhofer, Quasi-Steady-State Approximation for Chemical Reaction Networks, J. Math. Biol.36 (6) (1998) 593–609.[29] M. Rein, The Partial-Equilibrium Approximation in Reacting Flows, Phys. Fluids A: Fluid Dyn. 4 (5)(1992) 873–886.[30] J. Keck, Rate-Controlled Constrained-Equilibrium Theory of Chemical Reactions in Complex Systems,Prog. Energy Combust. Sci. 16 (2) (1990) 125–154.[31] S. Lam, D. Goussis, Understanding Complex Chemical Kinetics with Computational Singular Pertur-bation, in: Symp. Combust. Proc., Vol. 22, Elsevier, 1989, pp. 931–941.[32] S. Gupta, H. Im, M. Valorani, Analysis of n-Heptane Auto-Ignition Characteristics using ComputationalSingular Perturbation, Proc. Combust. Inst. 34 (1) (2013) 1125–1133.[33] U. Maas, S. Pope, Simplifying Chemical Kinetics: Intrinsic Low-Dimensional Manifolds in CompositionSpace, Combust. Flame 88 (3-4) (1992) 239–264.[34] J. C. Sutherland, A. Parente, Combustion modeling using principal component analysis, Proceedingsof the Combustion Institute 32 (1) (2009) 1563–1570. doi:https://doi.org/10.1016/j.proci.2008.06.147 .URL doi:https://doi.org/10.1016/j.combustflame.2014.11.027 .URL [36] O. Owoyele, T. Echekki, Toward Computationally Efficient Combustion DNS with Complex Fuels viaPrincipal Component Transport, Combust. Theory Model 21 (4) (2017) 770–798.[37] M. R. Malik, P. Obando Vega, A. Coussement, A. Parente, Combustion modeling using principal com-ponent analysis: A posteriori validation on sandia flames d, e and f, Proceedings of the CombustionInstitute doi:https://doi.org/10.1016/j.proci.2020.07.014 .URL [38] T. Sapsis, P. Lermusiaux, Dynamically orthogonal field equations for continuous stochastic dynamicalsystems, Physica D: Nonlinear Phenomena 238 (23-24) (2009) 2347–2360.[39] M. Cheng, T. Y. Hou, Z. Zhang, A dynamically bi-orthogonal method for time-dependent stochasticpartial differential equations i: Derivation and algorithms, Journal of Computational Physics 242 (0)(2013) 843 – 868. doi:http://dx.doi.org/10.1016/j.jcp.2013.02.033 .URL [40] P. Patil, H. Babaee, Real-time reduced-order modeling of stochastic partial differential equations viatime-dependent subspaces, Journal of Computational Physics 415 (2020) 109511. doi:https://doi.org/10.1016/j.jcp.2020.109511 .URL [41] H. Babaee, T. P. Sapsis, A minimization principle for the description of modes associated with finite-timeinstabilities, Proceedings of the Royal Society of London A: Mathematical, Physical and EngineeringSciences 472 (2186).URL http://dx.doi.org/10.1098/rspa.2015.0779 [42] M. Choi, T. P. Sapsis, G. E. Karniadakis, On the equivalence of dynamically orthogonal and bi-orthogonal methods: Theory and numerical simulations, Journal of Computational Physics 270 (2014)1 – 20. doi:http://dx.doi.org/10.1016/j.jcp.2014.03.050 .URL [43] O. Koch, C. Lubich, Dynamical low-rank approximation, SIAM Journal on Matrix Analysis and Appli-cations 29 (2) (2007) 434–454. doi:10.1137/050639703 .URL http://dx.doi.org/10.1137/050639703 [44] E. Musharbash, F. Nobile, T. Zhou, Error analysis of the dynamically orthogonal approximation oftime dependent random pdes, SIAM Journal on Scientific Computing 37 (2) (2015) A776–A810. doi:10.1137/140967787 .URL https://doi.org/10.1137/140967787 [45] H. Babaee, M. Farazmand, G. Haller, T. P. Sapsis, Reduced-order description of transient instabilitiesand computation of finite-time Lyapunov exponents, Chaos: An Interdisciplinary Journal of NonlinearScience 27 (6) (2017) 063103. doi:10.1063/1.4984627 .URL http://dx.doi.org/10.1063/1.4984627 [46] Z. Battles, L. N. Trefethen, An extension of matlab to continuous functions and operators, SIAM Journalon Scientific Computing 25 (5) (2004) 1743–1770.2547] M. Donello, M. Carpenter, H. Babaee, Computing sensitivities in evolutionary systems: A real-timereduced order modeling strategy (2020). arXiv:2012.14028 .[48] H. Babaee, An observation-driven time-dependent basis for a reduced description of transient stochas-tic systems, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences475 (2231) (2019) 20190506. doi:10.1098/rspa.2019.0506 .URL https://doi.org/10.1098/rspa.2019.0506 [49] N. Aubry, On the hidden beauty of the proper orthogonal decomposition 2 (5-6) (1991) 339–352. doi:10.1007/BF00271473 .URL [50] L. Alvergue, H. Babaee, G. Gu, S. Acharya, Feedback stabilization of a reduced-order model of a jet incrossflow, AIAA Journal 53 (9) (2015) 2472–2481. doi:10.2514/1.J053295 .URL http://dx.doi.org/10.2514/1.J053295 [51] P. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Me-chanics 656 (2010) 5–28.[52] J. Kutz, S. Brunton, B. Brunton, J. Proctor, Dynamic Mode Decomposition, Society for Industrial andApplied Mathematics, 2016. doi:doi:10.1137/1.9781611974508 .URL https://doi.org/10.1137/1.9781611974508 [53] K. Laksari, M. Kurt, H. Babaee, S. Kleiven, D. Camarillo, Mechanistic insights into human brainimpact dynamics through modal analysis, Physical Review Letters 120 (13) (2018) 138101–. doi:10.1103/PhysRevLett.120.138101 .URL https://link.aps.org/doi/10.1103/PhysRevLett.120.138101 [54] M. Velegar, N. B. Erichson, C. A. Keller, J. N. Kutz, Scalable diagnostics for global atmosphericchemistry using ristretto library (version 1.0), Geoscientific Model Development 12 (4) (2019) 1525–1539. doi:10.5194/gmd-12-1525-2019 .URL [55] H. Babaee, M. Choi, T. P. Sapsis, G. E. Karniadakis, A robust bi-orthogonal/dynamically-orthogonalmethod using the covariance pseudo-inverse with application to stochastic flow problems, Journal ofComputational Physics 344 (2017) 303–319. doi:https://doi.org/10.1016/j.jcp.2017.04.057 .URL [56] M. Anand, K. Rajagopal, K. Rajagopal, A model for the formation and lysis of blood clots, Pathophys-iology of haemostasis and thrombosis 34 (2-3) (2005) 109–120.[57] Z. Li, A. Yazdani, A. Tartakovsky, G. E. Karniadakis, Transport dissipative particle dynamics modelfor mesoscopic advection-diffusion-reaction problems, The Journal of chemical physics 143 (1) (2015)014101.[58] G. E. Karniadakis, S. J. Sherwin, Spectral/hp element methods for computational fluid dynamics, OxfordUniversity Press, USA, 2005.[59] H. Babaee, S. Acharya, X. Wan, Optimization of forcing parameters of film cooling effectiveness, Journalof Turbomachinery 136 (6) (2013) 061016–061016.URL http://dx.doi.org/10.1115/1.4025732 http://dx.doi.org/10.1115/1.4025562http://dx.doi.org/10.1115/1.4025562