Error estimates for a POD method for solving viscous G-equations in incompressible cellular flows
EError estimates for a POD method for solving viscous G-equations inincompressible cellular flows
Haotian Gu a , Jack Xin b , Zhiwen Zhang c, ∗ a Department of Mathematics, University of California at Berkeley, Berkeley, CA 94720, USA. b Department of Mathematics, University of California at Irvine, Irvine, CA 92697, USA. c Department of Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China.
Abstract
The G-equation is a well-known model for studying front propagation in turbulent combustion.In this paper, we develop an efficient model reduction method for computing regular solutionsof viscous G-equations in incompressible steady and time-periodic cellular flows. Our methodis based on the Galerkin proper orthogonal decomposition (POD) method. To facilitate the al-gorithm design and convergence analysis, we decompose the solution of the viscous G-equationinto a mean-free part and a mean part, where their evolution equations can be derived ac-cordingly. We construct the POD basis from the solution snapshots of the mean-free part.With the POD basis, we can efficiently solve the evolution equation for the mean-free part ofthe solution to the viscous G-equation. After we get the mean-free part of the solution, themean of the solution can be recovered. We also provide rigorous convergence analysis for ourmethod. Numerical results for viscous G-equations and curvature G-equations are presentedto demonstrate the accuracy and efficiency of the proposed method. In addition, we study theturbulent flame speeds of the viscous G-equations in incompressible cellular flows.
AMS subject classification:
Keywords:
Viscous G-equation; Hamilton-Jacobi type equation; front speed computation;cellular flows; proper orthogonal decomposition (POD) method; convergence analysis.
1. Introduction
Front propagation in turbulent combustion is a nonlinear and complicated dynamical process.The G-equation has been a very popular field model in combustion and physics literature forstudying premixed turbulent flame propagation [30, 31, 13, 36, 33, 45]. The G-equation modelis a sound phenomenological approach to study turbulent combustion, which uses the level-setformulation to study the flame front motion laws with the front width ignored. The simplestmotion law is that the normal velocity of the front is equal to a constant S l (the laminar speed)plus the projection of fluid velocity V ( x , t ) along the normal. This gives the inviscid G-equation G t + V · ∇ G + S l |∇ G | = 0 , (1) ∗ Corresponding author
Email addresses: [email protected] (Haotian Gu), [email protected] (Jack Xin), [email protected] (Zhiwen Zhang) a r X i v : . [ m a t h . NA ] N ov here the set { ( x , t ) : G ( x , t ) = 0 } corresponds to the location of the flame front at time t . Thederivation of the inviscid G-equation will be elaborated in detail in Section 2.1. Developingnumerical methods for (1) has been an active research topic for decades; see e.g. [34, 21, 12,26, 15] and references therein.As the fluid turbulence is known to cause stretching and corrugation of flames, additionalmodeling terms need to be incorporated into the basic G-equation. If the curvature term isadded into the basic equation to model the curvature effects and the curvature term is furtherlinearized, we will arrive at the viscous G-equation G t + V · ∇ G + S l |∇ G | = dS l ∆ G. (2)In order to compute numerical solutions of Eq.(2), the authors of [27] first approximated the G-equations by a monotone discrete system, then applied high resolution numerical methods suchas WENO (weighted essentially non-oscillatory) finite difference methods [20] with a combina-tion of explicit and semi-implicit time stepping strategies, depending on the size and propertyof dissipation in the equations. However, these existing numerical methods become expensivewhen one needs to discretize the Eq.(2) with a fine mesh or one needs to solve Eq.(2) manytimes with different parameters. This motivates us to exploit the low-dimensional structures ofthe regular solutions of viscous G-equation (2) and develop model reduction methods to solvethem efficiently.One of the successful model reduction ideas in the study of turbulent flows is the properorthogonal decomposition (POD) method [42, 7]. The POD method uses the data from anaccurate numerical simulation and extracts the most energetic modes in the system by usingthe singular value decomposition. This approach generates low-dimensional structures thatplay an important role in the dynamics of the flow. The Galerkin POD method has beenused to solve many types of partial differential equations, including linear parabolic equations,Burgers equations, Navier-Stokes equations and Hamilton–Jacobi–Bellman (HJB) equations;see [24, 25, 8, 44, 2, 22, 6, 23, 4] and references therein for details. The interested reader isreferred to [37, 6, 19] for a comprehensive introduction of the model reduction methods.In this paper, we study the POD method to solve the viscous G-equation (2), which is aHamilton-Jacobi type equation. To deal with the periodic boundary condition of the problemand facilitate the convergence analysis of the POD method, we decompose the solution of theviscous G-equation into a mean-free part and a mean part, where their evolution equationscan be derived accordingly; see Eq.(15). We construct the POD basis from the snapshots ofthe mean-free part of the solutions, which can be used to solve the evolution equation for themean-free part. Then, the mean part of the solution can be computed from the mean-freepart; see Eq.(14). Notice that the bilinear form of the evolution equation for the mean-freepart satisfies the coercive condition, which follows from the Poincar´e-Wirtinger inequality . Weprovide rigorous convergence analysis and show that the accuracy of our method is guaranteed.We remark that the idea of decomposing the data into a mean-free part and a mean part plays anessential role in the convergence analysis of the POD method for solving the viscous G-equation.Finally, we conduct numerical experiments to demonstrate the accuracy and efficiency of theproposed method.We find that the POD basis can be used to compute long-time solution of the viscousG-equation or the viscous G-equations with different parameters. Moreover, we study the tur-bulent flame speeds of viscous G-equations in incompressible steady and time periodic cellular2ows, which help our understanding of turbulent combustion. To further reduce the numericalerror, we develop an adaptive strategy to dynamically enrich the POD basis and demonstrateits effectiveness through a viscous G-equation with time-periodic fluid velocity. We remarkthat our POD method can easily be extended to solve other types of G-equations [27]. Thoughwe are unable to provide rigorous convergence analysis, we find that the POD method is stillefficient in solving curvature G-equation when its solution is smooth. To the best of our knowl-edge, our result is the first one in the literature that develops POD based model reductionmethod to solve G-equations and compute their front speeds.The rest of this paper will be organized as follows. In Section 2, we shall give a briefderivation of G-equation models. In Section 3, we show the detailed derivations of the modelreduction method for G-equations. In Section 4, we provide the convergence analysis of theproposed method. Our proof is based on the backward Euler-Galerkin-POD approximationscheme. The proofs for other discretization schemes, such as Crank-Nicolson scheme, can beobtained in a similar manner. In Section 5, we shall perform numerical experiments to test theperformance and accuracy of the proposed method. We find that POD method can provideconsiderable savings over finite difference methods in solving the G-equation while its numericalerror is relatively small. Finally, the concluding remarks will be given in Section 6. In theappendix Section, we provide the derivation of a finite difference scheme to solve G-equationsproposed in [27] and the procedure of constructing POD basis.
2. Turbulent combustion and G-equations
In this section, we briefly introduce the derivation of the G-equation in turbulent combus-tion. In a thin reaction zone regime and the corrugated flamelet regime of premixed tur-bulent combustion (Chapter 2 of [36]), the flame front is modeled by a level set function: { ( x , t ) : G ( x , t ) = 0 , x ∈ R n } , which is the interface between the burned area, denoted by { G < } and the unburned area, denoted by { G > } , respectively. Therefore, one can studythe propagation of the flame front by solving the dynamic equation for the level set function,namely the G-equation. Burned Fuel Unburned Fuel Flame Front V L s n ( , ) 0 G x t = G < G > Figure 1: Illustration of local interface velocities in the G-equation and a flame front in a 2D space. S l and the projection of fluid velocity V ( x , t ) along thenormal direction (see Fig.1). Hence, the trajectory of a particle x ( t ) on the interface satisfies,d x d t = V ( x , t ) + S l n , (3)where S l is the laminar flame speed and n is the normal vector. In terms of the level setfunction, the motion law (3) gives the inviscid G-equation, G t + V · ∇ G + S l |∇ G | = 0 , (4)where ∇ denotes the gradient operator. Thus, the normal vector in (3) can be computed using n = ∇ G/ |∇ G | . Notice that the set { G ( x , t ) = 0 } corresponds to the location of the flame frontat time t .To take into account the effect of flame stretching and corrugation, additional terms areadded into the inviscid G-equation (4) and one can obtain extended G-equation models involvingcurvature effects. The curvature G-equation is G t + V · ∇ G + S l |∇ G | = dS l |∇ G | (cid:0) ∇ · ∇ G |∇ G | (cid:1) , (5)where d is the Markstein number. The mean curvature term leads to a nonlinear degeneratesecond-order PDE that has been widely studied in the literature from both the analytical andnumerical point of views. The curvature dependent motion is well-known; see [14, 34, 33, 9]and references therein. However, we are not aware of a POD approximation of curvaturemotion in the general setting. We remark that the classical boundary conditions to deal withfirst-order level set equations and second-order mean curvature motion equations are Neumanntype boundary conditions. We are interested in computing the turbulent flame speeds of theG-equations in incompressible flows. Thus, we choose periodic boundary conditions in theG-equation (5), which will be further clarified in the numerical experiment; see Section 5.7.If the curvature term is linearized, we obtain the viscous G-equation as follows, which is theresearch focus in this paper. G t + V · ∇ G + S l |∇ G | = dS l ∆ G. (6)The linearized version of the curvature G-equation, i.e. the viscous G-equation (6) allows us tocarry out a convergent Galerkin approximation and POD error analysis, which is unavailablefor (5). Our POD method remains efficient for smooth solutions of the curvature G-equation(5), as shown numerically later. Now given a unit vector P ∈ R n , where n is the dimension of the physical space, we shallconsider the viscous G-equation (6) with planar initial condition (cid:40) G t + V · ∇ G + S l |∇ G | = dS l ∆ G in R n × (0 , ∞ ) ,G ( x ,
0) = P · x on R n × { t = 0 } . (7)4ere, we assume the flame front propagates in direction P with the initial front being { P · x = 0 } and x = ( x , x , ..., x n ). In addition, we assume V ( x , t ) is spatially periodic with C in x and C in t . Moreover, V ( x , t ) is divergence-free, i.e., ∇ x · V ( x , t ) = 0, ∀ t and uniformly bounded,i.e., || V ( x , t ) || ∞ ≤ A .If we write G ( x , t ) = P · x + u ( x , t ), then u ( x , t ) is also spatially periodic and satisfies thefollowing periodic initial value problem (cid:40) u t + V · ( P + ∇ u ) + S l | P + ∇ u | = dS l ∆ u in T n × (0 , ∞ ) ,u ( x ,
0) = 0 on T n × { t = 0 } , (8)where T n = [0 , n . Hence we can reduce the problem (7) defined on the whole domain R n tothe problem (8) defined on [0 , n by imposing the affine periodic condition (cid:40) G t + V · ∇ G + S l |∇ G | = dS l ∆ G in T n × (0 , ∞ ) ,G ( x ,
0) = P · x on T n × { t = 0 } , (9)with the assumption G ( x + z , t ) = G ( x , t ) + P · z .To illustrate the main idea of our method, we choose P = e = (1 , , ...,
0) and have G ( x , t ) = x + u ( x , t ). Let A ( t ) denote the volume of the burned area that has invaded duringtime interval (0 , t ). Then, the turbulent flame speed S T is defined as the linear growth rate of A ( t ). Notice that G ( x ,
0) = x and G ( x + e , t ) = G ( x , t ) + 1. Then, the S T can be evaluatedby G ( x , t ) or u ( x , t ) in [0 , n , S T = lim t →∞ − t (cid:90) [0 , n G ( x , t ) d x = lim t →∞ − t (cid:90) [0 , n ( x + u ( x , t )) d x . (10)In this work, we are interested in computing the turbulent flame speeds of the viscous G-equations in incompressible cellular flows, which is an effective (homogenized) quantity andinvolves a long-time computation. To this end, we impose periodic boundary conditions for (9).We will study the POD method for solving other nonlinear PDEs with other types of boundaryconditions in our future research.There are many numerical methods to solve the viscous G-equation (9) or Eq.(8) on abounded physical domain T n . For instance, we can apply a higher order Hamilton-Jacobiweighted essentially non-oscillatory (HJ-WENO) scheme and total variation diminishing Runge-Kutta (TVD-RK) scheme in the spatial and time discretization, respectively. See [41, 20, 33, 27]for details of the schemes. To make our paper self-contained, we show the numerical schemesproposed by [27]; see the Appendix Section for more details.
3. Model reduction method for G-equations
In practice, we are interested in studying the dependence of the turbulent flame speed S T ondifferent parameters in the viscous G-equation, such as the Markstein number d . As such, wehave to solve the viscous G-equation many times, which is an expensive task. Therefore, weneed to design numerical methods that allow us to efficiently and accurately solve the viscousG-equation. We shall develop an efficient model reduction method to achieve this goal. Wepoint out that the model reduction methods also bring computational savings when one needsto solve high-dimensional problems, e.g. 3D viscous G-equations.5 .1. A decomposition strategy According to the definition in (10), we can either solve Eq.(8) to obtain u ( x, t ) or solve Eq.(9)to obtain G ( x, t ), in order to compute the turbulent flame speed S T . We shall consider to solveEq.(8) since it is easier to deal with the boundary condition. Let us decompose the solution u of Eq.(8) into ˆ u + ¯ u , where ˆ u is the mean-free part and ¯ u is the mean of u . This decompositionmeans that (cid:90) T n ˆ u ( x , t ) d x = 0 ∀ t and ¯ u ( t ) = (cid:90) T n u ( x , t ) d x . (11)Then, substituting u = ˆ u + ¯ u into (8), we obtainˆ u t + ¯ u t + V · (cid:0) P + ∇ ˆ u (cid:1) + S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) − dS l ∆ˆ u = 0 . (12)Integrating the above equation over the domain T n gives (cid:90) T n (cid:2) ˆ u t + ¯ u t + V · ( P + ∇ ˆ u ) − dS l ∆ˆ u (cid:3) d x = − (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x . (13)The definitions of ¯ u and ˆ u in (11) imply that (cid:82) T n ˆ u t = 0 and (cid:82) T n ¯ u t = ¯ u t . Furthermore, theperiodic conditions of V and ˆ u imply (cid:82) T n V · P = 0, (cid:82) T n ∆ˆ u = 0, and (cid:82) T n V · ∇ ˆ u = 0, where wehave used the facts that V is divergence-free and ˆ u is the mean-free. Combining these results,we can simplify the Eq.(13) as, ¯ u t = − (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x . (14)Finally, we find that the Eq.(8) is equivalent to ˆ u t + V · ∇ ˆ u − dS l ∆ˆ u + S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) − (cid:82) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x + V · P = 0 in T n × (0 , ∞ )¯ u t = − (cid:82) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , t ) (cid:12)(cid:12) d x on t ∈ (0 , ∞ )ˆ u ( x ,
0) = 0 on T n × { t = 0 } ¯ u (0) = 0 (15)The strategy of decomposing the solution or data into a mean-free part and a mean part isoften used in scientific computing and data science. However, it has not been studied in theG-equation setting. We shall demonstrate later that this decomposition plays a crucial role inthe algorithm design and facilitates the convergence analysis of our proposed method. In this section, we shall present our model reduction method to solve Eq.(15). Since theevolution equation for ¯ u depends on ˆ u , we first consider to solve the equation for ˆ u , i.e.,ˆ u t + V · ∇ ˆ u − dS l ∆ˆ u + S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) − (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x + V · P = 0 , in T n × (0 , ∞ ) . (16)Let H per ( T n ) denote the Sobolev space on the domain T n with a periodic boundary conditionand let (cid:104)· , ·(cid:105) denote the standard inner product on L ( T n ). Let H ⊂ H per ( T n ) be the subspace6onsisting of all mean-free functions. Since H is a closed subspace of H per ( T n ), H itself is aHilbert space. Let (cid:104)· , ·(cid:105) H denote the standard inner product on H . Define the bilinear form a ( · , · ) : H × H → R to be a ( u, v ) = (cid:90) T n (cid:2) ( V · ∇ u ) v + dS l ( ∇ u · ∇ v ) (cid:3) d x . (17)Also, define a nonlinear map from H to L ( T n ) to be F ( u ) = S l (cid:12)(cid:12) P + ∇ u (cid:12)(cid:12) − (cid:90) T n S l (cid:12)(cid:12) P + ∇ u (cid:12)(cid:12) d x . (18)The weak formulation associated with the Eq.(16) is (cid:104) ˆ u t , ψ (cid:105) + a (ˆ u, ψ ) + (cid:104) F (ˆ u ) , ψ (cid:105) = (cid:104)− P · V , ψ (cid:105) , ∀ ψ ∈ H, (19)which can be solved by using numerical methods.The POD technique for discretization of Eq.(16) requires snapshots, which can be obtainedby an independent numerical method or by the appropriate technological means related to aspecific application, for instance the data from an experiment. We apply a higher-order WENOscheme and TVD-RK scheme to solve the problem (8), and extract the mean part to obtaina set of numerical solutions { ˆ u ( · , t k ) } to Eq.(16), where t k = k ∆ t , ∆ t = T /m , k = 0 , ..., m .Then, we use 2 m + 1 snapshots { ˆ u ( · , t ) , . . . , ˆ u ( · , t m ) , ¯ ∂ ˆ u ( · , t ) , . . . , ¯ ∂ ˆ u ( · , t m ) } , where ¯ ∂ ˆ u ( t i ) =(ˆ u ( t i ) − ˆ u ( t i − )) / ∆ t , to generate the POD basis functions. Let S r = { ψ ( · ) , ψ ( · ) , · · · , ψ r ( · ) } denote the r -dimensional orthonormal POD basis functions obtained from the 2 m +1 snapshots,which minimize the following error12 m + 1 m (cid:88) i =0 (cid:13)(cid:13)(cid:13) ˆ u ( t i ) − r (cid:88) j =1 (cid:104) ˆ u ( t i ) , ψ j (cid:105) H ψ j (cid:13)(cid:13)(cid:13) H + 12 m + 1 m (cid:88) i =1 (cid:13)(cid:13)(cid:13) ¯ ∂ ˆ u ( t i ) − r (cid:88) j =1 (cid:104) ¯ ∂ ˆ u ( t i ) , ψ j (cid:105) H ψ j (cid:13)(cid:13)(cid:13) H . (20)See the appendix 7.3 for the details of the POD method. By adding the terms ¯ ∂ ˆ u ( t i ), i = 1 , ..., m into the snapshots, we obtain more accurate POD basis functions and avoid an extra (∆ t ) − factor in the convergence analysis; see Section 4 for more details.In practice, we can use experimental data or reference numerical solutions to generate so-lution snapshots. The choice of snapshots is critical to the quality of the POD method. Thesnapshots should span a linear space that approximates the solution space of the original PDEswell. However, the POD method itself gives no guidance on how to select the snapshots (page503 of [6]). One may need some a posteriori error estimate for the solution obtained by thePOD method. For time-dependent parametric PDEs, it is difficult to obtain a posteriori errorestimate for the PDEs. We propose an adaptive strategy to enrich the POD basis; see Section5.6 for more details. We refer the interested reader to [18, 39, 1, 35, 38, 28, 17, 29] and referencetherein for the adaptive techniques in the model reduction methods. Remark . The construction of the POD basis can be costly. One can apply randomizedprojection methods [3] to reduce the cost in the offline stage. Once the construction is done,the POD basis can be used to solve viscous G-equations with different parameters, which bringssignificant computational savings in the online stage.7 .3. A backward Euler and POD-based Galerkin method
The POD basis functions provide an efficient approach to approximate the solution in thephysical space. If we choose the backward Euler scheme to discretize the time space, weobtain a backward Euler and POD-based Galerkin method to solve Eq.(16). Specifically, letˆ U k ≡ (cid:80) ri =1 a i ( t k ) ψ i denote the numerical solution at t = t k , where ψ i ’s are the POD basisfunctions. We want to find solutions { ˆ U k } mk =0 ⊂ S r satisfying (cid:40) (cid:104) ¯ ∂ ˆ U k , ψ (cid:105) + a ( ˆ U k , ψ ) + (cid:104) F ( ˆ U k ) , ψ (cid:105) = (cid:104)− P · V , ψ (cid:105) , ∀ ψ ∈ S r , ˆ U = 0 , (21)where ¯ ∂ ˆ U k = ( ˆ U k − ˆ U k − ) / ∆ t . By choosing the test function ψ to be ψ i , i = 1 , ..., r in (21) andletting a k = (cid:0) a ( t k ) , · · · , a r ( t k ) (cid:1) T denote the coefficient vector, we obtain a nonlinear equationsystem for a k as M a k = M a k − + c + f k , (22)where M , M ∈ R r × r with ( M ) ij = (cid:104) ψ i , ψ j (cid:105) + ∆ t · a ( ψ i , ψ j ) and ( M ) ij = (cid:104) ψ i , ψ j (cid:105) , c ∈ R r with c i = − ∆ t (cid:104) V · P , ψ i (cid:105) , and F k ∈ R r with ( f k ) i = − ∆ t (cid:10) F ( (cid:80) ri a i ( t k ) ψ i ) , ψ i (cid:11) . The matrices M , M and c can be pre-computed and saved.The Eq.(22) can be efficiently computed using the Newton’s method, where the solutionat time t k − can be chosen as the initial guess for a k . The proposed method is very efficientsince the number of the POD basis is small in general. For more challenging problems, such as3D problems or convection-dominated problems (leading to an increase of the number of thePOD basis), the effect of the nonlinear term becomes significant. One can choose the discreteempirical interpolation method (DEIM) for nonlinear model reduction [10].We solve Eq.(22) to get the POD solution ˆ U k , which is the numerical solution of the mean-free part ˆ u ( x , t ) of the G-equation (16). As indicated by the Eq.(15), u ( x , t ) is recovered byadd the mean solution back u ( x , t ) = ˆ u ( x , t ) − (cid:90) t (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , s ) (cid:12)(cid:12) dxds . (23)Since the solution ˆ u ( x , t ) is smooth, we use a high-order difference scheme to compute thespatial derivatives of ∇ ˆ u ( x , · ) and use extrapolations at boundaries to approximate the spatialderivatives to maintain the second-order accuracy. Then, we apply a composite trapezoidalrule to compute the spatial integration and temporal integration in (23). Finally, we obtain thenumerical solution to u ( x , t ) as follows U k = ˆ U k −
12 ∆ t k (cid:88) i =1 (cid:20) I (cid:0) S l | P + ∇ ˆ U i − ( x ) | (cid:1) + I (cid:0) S l | P + ∇ ˆ U i ( x ) | (cid:1)(cid:21) , (24)where I ( · ) denotes the numerical integrator by the composite trapezoidal rule. Moreover, weassume the mesh size of the composite trapezoidal rule is h , which is determined by the numeri-cal method in computing the snapshots. One can also choose a wider finite difference scheme tocompute the spatial derivatives and high-order numerical methods to compute the integrationin (23). 8 . Convergence analysis In this section, we shall present some convergence analysis to show that the accuracy of the nu-merical solution is guaranteed. Our convergence analysis follows the framework of the Galerkinfinite element methods for parabolic problems [43]. To deal with the nonlinearity of the viscousG-equation, the following lemmas are useful.
Lemma 4.1.
There exists a constant γ > such that for any u, v ∈ H , || F ( u ) − F ( v ) || L ≤ γ || u − v || H . (25) Proof.
According to the definition of F in Eq.(18), the proof of this lemma can be easilyobtained by using the triangle inequality. Lemma 4.2.
The bilinear form a ( · , · ) defined in (17) is continuous and coercive, which meansthat there exist constants β > and κ > such that for any ψ, φ ∈ Ha ( ψ, φ ) ≤ β || ψ || H || φ || H , a ( ψ, ψ ) ≥ κ || ψ || H . (26) Proof.
Let || V || ∞ denote the maximum amplitude of the vector field V . One has the estimate a ( ψ, φ ) ≤ (cid:90) T n || V || ∞ |∇ ψ | · | φ | + dS l |∇ ψ | · |∇ φ | , ≤ || V || ∞ || ψ || H || φ || H + dS l || ψ || H || φ || H . The last inequality follows from the Cauchy-Schwarz inequality. Moreover, since V isdivergence-free, V · ∇ is skew-symmetric, which means a ( ψ, ψ ) = dS l (cid:90) T n |∇ ψ | . Since ψ is mean-free, the Poincar´e-Wirtinger inequality implies that there exists a constant C ,depending only on the domain T n , such that (cid:82) T n | ψ | ≤ C (cid:82) T n |∇ ψ | . Therefore, we know that a ( · , · ) is coercive, i.e., a ( ψ, ψ ) ≥ κ || ψ || H , where κ = dS l / ( C + 1) > P r : H → S r , u (cid:55)→ P r u such that a ( P r u, ψ ) = a ( u, ψ ) , ∀ ψ ∈ S r . (27)Facts from functional analysis guarantee that P r is well-defined and bounded because a ( · , · )is continuous and coercive. More specifically, 9 | P r u || H ≤ βκ || u || H , ∀ u ∈ H. (28)Using the same argument in Lemma 3 and Lemma 4 in [24], we can prove that P r hasthe following approximation property. More details of the approximation property of the PODbasis can be found in the appendix 7.3. Lemma 4.3. m m (cid:88) k =0 || ˆ u ( t k ) − P r ˆ u ( t k ) || H ≤ βκ m (cid:88) l = r +1 λ l , (29) and m m (cid:88) k =1 || ¯ ∂ ˆ u ( t k ) − P r ¯ ∂ ˆ u ( t k ) || H ≤ βκ m (cid:88) l = r +1 λ l , (30) where λ l is the l -th largest eigenvalues of the correlation matrix K associated with the solutionsnapshots. Theorem 4.4.
Let ˆ u ( · ) and { ˆ U k } mk =0 be the solutions to Eq. (19) and its backward Euler andPOD-based Galerkin approximation, respectively. Then for sufficiently small ∆ t , there exists aconstant C > depending on ˆ u , d , S l , V , P and T but independent of r , ∆ t , and m such that m m (cid:88) k =1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ U k − ˆ u ( t k ) (cid:12)(cid:12)(cid:12)(cid:12) L ≤ C (cid:0) (∆ t ) + m (cid:88) l = r +1 λ l (cid:1) . (31) Proof.
For k = 0 , , . . . , m , define ϑ k = ˆ U k − P r ˆ u ( t k ) and (cid:37) k = P r ˆ u ( t k ) − ˆ u ( t k ). Then1 m m (cid:88) k =1 || ˆ U k − ˆ u ( t k ) || L ≤ m m (cid:88) k =1 || ϑ k || L + 2 m m (cid:88) k =1 || (cid:37) k || L , (32)2 m m (cid:88) k =1 || (cid:37) k || L ≤ m m (cid:88) k =1 || (cid:37) k || H ≤ βκ m (cid:88) l = r +1 λ l . (33)Define ¯ ∂ϑ k = ( ϑ k − ϑ k − ) / ∆ t . For all ψ ∈ S r , (cid:104) ¯ ∂ϑ k , ψ (cid:105) + a ( ϑ k , ψ ) = (cid:104) υ k , ψ (cid:105) + (cid:104) F (ˆ u ( t k )) − F ( ˆ U k ) , ψ (cid:105) , (34)where υ k = ˆ u t ( t k ) − P r ¯ ∂ ˆ u ( t k ) = ˆ u t ( t k ) − ¯ ∂ ˆ u ( t k ) + ¯ ∂ ˆ u ( t k ) − P r ¯ ∂ ˆ u ( t k ). Define ω k = ˆ u t ( t k ) − ¯ ∂ ˆ u ( t k )and η k = ¯ ∂ ˆ u ( t k ) − P r ¯ ∂ ˆ u ( t k ). Taking ψ = ϑ k ∈ S r in the previous equality, we obtain1∆ t (cid:0) || ϑ k || L − (cid:104) ϑ k , ϑ k − (cid:105) (cid:1) + κ || ϑ k || L ≤ || F (ˆ u ( t k )) − F ( ˆ U k ) || L || ϑ k || L + || ϑ k || L || υ k || L . (35)By using the Lemma 4.1, we get || F (ˆ u ( t k )) − F ( ˆ U k ) || L || ϑ k || L ≤ γ || ˆ u ( t k ) − ˆ U k || H || ϑ k || L , ≤ γ (cid:0) || (cid:37) k || H + || ϑ k || H (cid:1) || ϑ k || L . ϑ k ∈ S r , which is a finite dimensional space, the norms defined on S r are equivalent. Thismeans that there exists some constant C > || F (ˆ u ( t k )) − F ( ˆ U k ) || L || ϑ k || L ≤ C (cid:0) || (cid:37) k || H + || ϑ k || L (cid:1) || ϑ k || L , ≤ C || (cid:37) k || H + 3 C || ϑ k || L . Combining these inequalities, we have || ϑ k || L − (cid:104) ϑ k , ϑ k − (cid:105) ≤ ∆ t (cid:16) || ϑ k || L || υ k || L + C || (cid:37) k || H + 3 C || ϑ k || L (cid:17) . By using the inequality of arithmetic and geometric means, we obtain (cid:0) − (1 + C )∆ t (cid:1) || ϑ k || L ≤ || ϑ k − || L + ∆ t (cid:0) || υ k || L + 3 C || (cid:37) k || H (cid:1) . For sufficiently small ∆ t , there exists a constant C > || ϑ k || L ≤ (1 + C ∆ t ) (cid:18) || ϑ k − || L + ∆ t (cid:0) || υ k || L + 3 C || (cid:37) k || H (cid:1)(cid:19) . (36)By iteratively using the inequality (36), we have || ϑ k || L ≤ e C T (cid:18) || ϑ || L + ∆ t k (cid:88) j =1 (cid:0) || υ j || L + 3 C || (cid:37) j || H (cid:1)(cid:19) . (37)Note that ϑ = 0 in our case. Therefore, we sum the inequality (37) from k = 1 to m andarrive at m (cid:88) k =1 || ϑ k || L ≤ e C T ∆ t m (cid:88) k =1 k (cid:88) j =1 (cid:0) || υ j || L + 3 C || (cid:37) j || H (cid:1) , ≤ e C T T m (cid:88) j =1 (cid:0) || υ j || L + 3 C || (cid:37) j || H (cid:1) , ≤ e C T T m (cid:88) j =1 (cid:0) || ω j || L + 2 || η j || L + 3 C || (cid:37) j || H (cid:1) . Therefore,1 m m (cid:88) k =1 || ˆ U k − ˆ u ( t k ) || L ≤ m m (cid:88) k =1 || ϑ k || L + 6 βκ m (cid:88) k = r +1 λ k , ≤ e C T Tm m (cid:88) j =1 (cid:0) || ω j || L + 2 || η j || L + 3 C || (cid:37) j || H (cid:1) + 6 βκ m (cid:88) l = r +1 λ l . ω j = ˆ u t ( t j ) − ¯ ∂ ˆ u ( t j ) = t (cid:82) t j t j − ( t − t j − )ˆ u tt ( t ) dt , we obtain that m (cid:88) j =1 || ω j || L ≤ m (cid:88) j =1 t ) (cid:90) t j t j − ( t − t j − ) dt (cid:90) t j t j − || ˆ u tt ( t ) || L dt ≤ ∆ t (cid:90) T || ˆ u tt ( t ) || L . (38)Together with Lemma 4.3, we obtain1 m m (cid:88) k =1 || ˆ U k − ˆ u ( t k ) || L ≤ (cid:18) e C T (cid:90) T || ˆ u tt ( t ) || L (cid:19) (∆ t ) + C m (cid:88) l = r +1 λ l , where C = βκ (cid:0) e C T T (2 + 3 C ) + 2 (cid:1) and the fact ∆ t = T /m is used. This completes the prooffor Theorem 4.4.Theorem 4.4 provides an error estimate for the mean-free part. The mean part depends onthe mean-free part and the original solution can be recovered afterward. Thus, we obtain theerror estimate for the original solution as follows.
Theorem 4.5.
Let u ( · ) and { U k } mk =0 be the solutions to Eq. (8) and its numerical approximationEq. (24) based on the backward Euler and POD-based Galerkin scheme, respectively. Then forsufficiently small ∆ t , there exists a constant C (cid:48) ≥ depending on β , κ , ˆ u , d , S l , V , P and T but independent of r , ∆ t , h , and m , such that m m (cid:88) k =1 || U k − u ( t k ) || L ≤ C (cid:48) (cid:0) (∆ t ) + m (cid:88) l = r +1 λ l + h (cid:1) . (39) Proof.
Let C be the constant appearing in Theorem 4.4; see Eq.(31). Recall that u ( t k ) = ˆ u ( t k ) + ¯ u ( t k ) = ˆ u ( t k ) − (cid:90) t k (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , t ) (cid:12)(cid:12) d x dt (40)and U k = ˆ U k + ¯ U k = ˆ U k −
12 ∆ t k (cid:88) i =1 (cid:20) I (cid:0) S l | P + ∇ ˆ U i − ( x ) | (cid:1) + I (cid:0) S l | P + ∇ ˆ U i ( x ) | (cid:1)(cid:21) , (41)where I ( · ) is the numerical integrator by the composite trapezoidal rule. We also denote˜ U i = − I (cid:0) S l | P + ∇ ˆ U i ( x ) | (cid:1) , i = 0 , , ..., k . Then, we obtain1 m m (cid:88) k =1 (cid:12)(cid:12)(cid:12)(cid:12) U k − u ( t k ) (cid:12)(cid:12)(cid:12)(cid:12) L ≤ m m (cid:88) k =1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ U k − ˆ u ( t k ) (cid:12)(cid:12)(cid:12)(cid:12) L + 2 m m (cid:88) k =1 (cid:12)(cid:12) ¯ u ( t k ) − ¯ U k (cid:12)(cid:12) . (42)From Theorem 4.4, we get an estimate for the first part on the RHS of the inequality (42).In the sequel, we shall estimate the second part of the RHS of the inequality (42). We considerthe following decomposition 12 (cid:88) k =1 (cid:12)(cid:12) ¯ u ( t k ) − ¯ U k (cid:12)(cid:12) ≤ m (cid:88) k =1 (cid:12)(cid:12)(cid:12) ¯ u ( t k ) − k (cid:88) j =1
12 ∆ t (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1)(cid:12)(cid:12)(cid:12) + 2 m (cid:88) k =1 (cid:12)(cid:12)(cid:12) k (cid:88) j =1
12 ∆ t (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1) − k (cid:88) j =1
12 ∆ t (cid:0) ˜ U j − + ˜ U j (cid:1)(cid:12)(cid:12)(cid:12) . (43)The first term on the RHS of the inequality (43) is bounded by2 m (cid:88) k =1 (cid:12)(cid:12)(cid:12) ¯ u ( t k ) − k (cid:88) j =1
12 ∆ t (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1)(cid:12)(cid:12)(cid:12) , ≤ t ) m (cid:88) k =1 (cid:12)(cid:12)(cid:12) k (cid:88) j =1 (cid:16) ¯ u ( t j ) − ¯ u ( t j − )∆ t − (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1)(cid:17)(cid:12)(cid:12)(cid:12) , ≤ t ) m (cid:88) k =1 k k (cid:88) j =1 (cid:12)(cid:12)(cid:12) ¯ u ( t j ) − ¯ u ( t j − )∆ t − (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1)(cid:12)(cid:12)(cid:12) , ≤ t ) m (cid:88) k =1 k (∆ t ) (cid:90) T (cid:12)(cid:12)(cid:12)(cid:12) ¯ u ttt ( t ) (cid:12)(cid:12)(cid:12)(cid:12) L ≤ T
12 (∆ t ) (cid:90) T | ¯ u ttt ( t ) | L . (44)The second term on the RHS of Eq.(43) is bounded by2 m (cid:88) k =1 (cid:12)(cid:12)(cid:12) k (cid:88) j =1
12 ∆ t (cid:0) ¯ u t ( t j − ) + ¯ u t ( t j ) (cid:1) − k (cid:88) j =1
12 ∆ t (cid:0) ˜ U j − + ˜ U j (cid:1)(cid:12)(cid:12)(cid:12) ≤ (∆ t ) m (cid:88) k =1 k k (cid:88) j =1 (cid:16)(cid:12)(cid:12) ¯ u t ( t j − ) − ˜ U j − (cid:12)(cid:12) + (cid:12)(cid:12) ¯ u t ( t j ) − ˜ U j (cid:12)(cid:12) (cid:17) , = (∆ t ) m (cid:88) k =1 k k (cid:88) j =1 (cid:16)(cid:12)(cid:12) (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , t j − ) d x − I (cid:0) S l | P + ∇ ˆ U j − ( x ) | (cid:1)(cid:12)(cid:12) + (cid:12)(cid:12) (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , t j ) d x − I (cid:0) S l | P + ∇ ˆ U j ( x ) | (cid:1)(cid:12)(cid:12) (cid:17) , ≤ t ) m m (cid:88) j =0 (cid:16) S l (cid:90) T n (cid:12)(cid:12) ∇ ˆ u ( x , t j ) − ∇ ˆ U j ( x ) (cid:12)(cid:12) d x + (cid:12)(cid:12)(cid:12) (cid:90) T n S l (cid:12)(cid:12) P + ∇ ˆ U j ( x ) (cid:12)(cid:12) d x − I (cid:0) S l | P + ∇ ˆ U j ( x ) | (cid:1)(cid:12)(cid:12)(cid:12) (cid:17) , ≤ T S l m (cid:88) j =0 (cid:16) || ˆ u ( t j ) − ˆ U j || H + C h (cid:17) ≤ T S l ( 3 βκ m (cid:88) l = r +1 λ l + mC h ) , (45)where C comes from the error estimate of the trapezoidal rule and depends on third-orderderivatives of the solution with respect to physical variables, and the last inequality followsfrom the fact that ˆ u ( t j ) − ˆ U j ∈ S r and norms are equivalent in a finite dimensional space.Substituting the estimate results (31), (44) and (45) into (42), we obtain13 m m (cid:88) k =1 || U k − u ( t k ) || L ≤ C (cid:48) (cid:0) (∆ t ) + m (cid:88) l = r +1 λ l + h (cid:1) , (46)where C (cid:48) ≥ β , κ , ˆ u , d , S l , V , P and T , but is independent of r , ∆ t , h , and m .Theorem 4.4 and Theorem 4.5 show that the errors depend on the sum of eigenvalues of thecorrelation matrix K (associated with the solution snapshots) that are not included in the PODbasis. Thus, the decay of eigenvalues provides a quantitative guidance for choosing the numberof the POD basis. A typical approach is to choose r so that (cid:80) ml = r +1 λ l ≤ min { (∆ t ) , h } .
5. Numerical Results
We shall perform numerical experiments to test the performance and accuracy of the proposedmethod. The numerical experiments are all carried out on a laptop with a 2.5GHz quad-coreIntel Core i7 processor and 8GB RAM. The Python codes are published on GitHub. Weconsider the following viscous G-equation on [0 , with planar initial condition (cid:40) G t + V · ∇ G + S l |∇ G | = dS l ∆ G in [0 , × (0 , T ) ,G ( x ,
0) = P · x on [0 , × { t = 0 } , (47)where x = ( x, y ), P is the flame front propagation direction, and the assumption G ( x + z , t ) = G ( x , t ) + P · z is used. The setting of the fluid velocity V ( x ) is an important issue in turbulentcombustion modeling. We first consider a steady incompressible cellular flow, V ( x ) = ( V , V ) = ∇ ⊥ H = ( −H y , H x ) , H = A π sin(2 πx ) sin(2 πy ) , (48)where A is the amplitude of the velocity filed. In the following numerical experiments on thissteady flow (48), we choose A = 4 . S l = 1, and P = (1 , , into ( N h +1) × ( N h +1)grids with mesh size h = 1 /N h and N h = 80. For the time discretization, we choose time step∆ t = 1 /N t with N t = 1000. As such, solving the viscous G-equation from 0s to T s will generate2 T N t + 1 snapshots.In the POD scheme, the algorithm in the Section 3.2 is first applied to construct a setof POD basis, denoted by S = { ψ , · · · , ψ r } , r ≤ min { ( N h + 1) × ( N h + 1) , T N t + 1 } .In practice, we choose the number of the POD basis r to be the smallest integer such that (cid:80) N snap l = r +1 λ l / (cid:80) N snap l =1 λ l ≤ e P OD , where λ l ’s are the eigenvalues of the covariance matrix of thesolution snapshots, N snap is the number of solution snapshots, and the error threshold e P OD is https://github.com/Harry970804/POD G-Equation. H norm incomputing the POD basis. Usually, r is assumed to be much smaller than the degree of freedomin the physical space discretization. In each of the following experiments, r will be listed indetail later. But we want to point out in advance that, under most of the parameter settingswe have experimented, a small r (around 10) will be good enough to capture true solutions wellwith a high accuracy. With the POD basis, we further apply the scheme in Section 3.3 to getthe POD solution. The time discretization for the POD scheme is ∆ t P OD = 1 / In the first numerical experiment, we choose d = 0 . t = 0 to 1s using the reference numerical scheme and obtain mean-free solution snapshots.Then, the POD basis functions are constructed from the entire snapshots. In this case, thenumber of basis functions is r = 4. Finally, we compute the (mean-free) POD solution on thesame time period.We compare the mean-free parts obtained by two difference methods at T = 1 .
0. Fig.2shows that the POD solution agrees well with the reference solution. Moreover, based on themean-free POD solution, we recover the numerical solution of the G-equation according toEq.(24). In Fig.3, we show the recovered solution using our POD method and the referencesolution. We can see the good performance of the POD method.We continue our experiment by computing Eq.(47) with different choices of d from 0 .
01 to0 .
1. We compare the errors between the POD method and the reference method in computingthe mean-free part of the solution, the mean of the solution and the recovered solution, respec-tively. Table 1 shows the relative error in L norm of these results with different choices of d ’s, where “Rela. error, Mean-free”, “Rela. error, Mean”, and “Rela. error, Recovered” meanrelative error of the mean-free part of the solution, relative error of the mean of the solutionand relative error of the recovered solution, respectively. We can see that: (1) the relative errorof the mean-free component is small, which shows the good performance of the POD method;(2) the relative error of the mean part is much smaller than that of the mean-free part, whichwe think is due to the cancellation of the error in computing the mean; and (3) the relativeerror of the recovered solution is smaller than that of the mean-free part since the recoveredsolution usually has a larger magnitude than the mean-free part.For the dimension of the subspace spanned by the POD basis, we observe that under a fixedthreshold of e P OD , in general, it will increase as d decreases. But it still remains small and itfurther implies the efficiency of the POD method. In all the experiments we have reported, thedimensions vary in the range from 4 to 10.In Fig.4 and Fig.5, we show the comparison between the recovered solutions obtained usingthe POD method and the reference solutions for d = 0 .
05 and d = 0 .
01, respectively. Thoughwe observe some sharp layers appearing in the solution as d decreases, the POD basis cancapture these layered structures and give accurate numerical results.These results show that the POD method can capture the low-dimensional structures in theregular solutions of the viscous G-equations and provide an efficient model reduction methodto approximate the solutions. We remark that when d is extremely small or even d = 0 thesolutions of G-equations are not smooth and cannot be computed by the POD-based Galerkinmethod. We will exploit some features of the discontinuous Galerkin method and develop newPOD methods to address this challenge in our future work.15 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−0.2−0.10.00.10.2 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.00050.00100.00150.00200.00250.00300.00350.0040 (c)
Absolute difference.
Figure 2: Mean-free component of the solution at T = 1 with d = 0 . x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−1.2−1.0−0.8−0.6−0.4 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−1.2−1.0−0.8−0.6−0.4 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0020.0030.0040.0050.0060.0070.008 (c)
Absolute difference.
Figure 3: Solution of the viscous G-equation at T = 1 with d = 0 . In this numerical experiment, we first solve the viscous G-equation (47) from T = 0 to T = 1using the finite difference scheme to construct POD basis. Then, we adopt the POD methodto solve the equation on a longer time from T = 0 to T = 2 using that POD basis.We conduct the experiments for d = 0 . , . , . . . , .
1. The number of POD basis functionsis fixed to be 6 across all different d .In Fig.6 and Fig.7, the results for d = 0 . d .In addition, for d = 0 .
05 and d = 0 .
01, we plot their recovered solutions obtained fromboth methods in Fig.8 and Fig.9 at T = 2, respectively. Again, we find that the POD methodperforms well for those 2 settings. Moreover, the profiles in Figs.8 and 9 can also be viewed16 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−1.6−1.4−1.2−1.0−0.8−0.6 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0020.0040.0060.008 (c)
Absolute difference.
Figure 4: Solution of the viscous G-equation at T = 1 with d = 0 . x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−2.2−2.0−1.8−1.6−1.4−1.2 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−2.2−2.0−1.8−1.6−1.4−1.2 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0200.0250.0300.0350.0400.045 (c)
Absolute difference.
Figure 5: Solution of the viscous G-equation at T = 1 with d = 0 . Rela. error, Mean
Rela. error, Recovered d 0.06 0.07 0.08 0.09 0.1Rela. error, Mean-free
Rela. error, Mean
Rela. error, Recovered
Table 1: The relative errors between POD solution and reference solution. approximately as downward shifts of the profiles in Figs.4 and 5, respectively. The smaller the d is, the bigger shift will be observed. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−0.2−0.10.00.10.2 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−0.2−0.10.00.10.2 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0010.0020.0030.0040.0050.0060.007 (c)
Absolute difference.
Figure 6: Mean-free component of the solution at T = 2 with d = 0 . d 0.01 0.02 0.03 0.04 0.05Rela. error, Mean-free Rela. error, Recovered d 0.06 0.07 0.08 0.09 0.1Rela. error, Mean-free
Rela. error, Recovered
Table 2: The relative errors between POD solution and reference solution at T = 2. In this subsection, we investigate the robustness of the POD basis across different diffusionparameter d ’s. Specifically, we build the POD basis from the solution snapshots of the viscousG-equation with the diffusivity d . Then, we use that pre-computed POD basis to computesolutions of the viscous G-equation with other d ’s.In our numerical experiment, we solve the viscous G-equation from T = 0 to T = 1 with d = 0 .
05 to build the POD basis. In Table 3, we show the relative errors between the POD18 (a)
Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−2.6−2.4−2.2−2.0−1.8 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0160.0180.0200.0220.0240.0260.028 (c)
Absolute difference.
Figure 7: Solution of the viscous G-equation at T = 2 with d = 0 . x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−3.0−2.8−2.6−2.4−2.2 (a) Reference solution x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−3.2−3.0−2.8−2.6−2.4−2.2 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.0600.0650.0700.0750.080 (c)
Absolute difference.
Figure 8: Solution of the viscous G-equation at T = 2 with d = 0 . x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−4.4−4.2−4.0−3.8−3.6−3.4 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−4.4−4.2−4.0−3.8−3.6−3.4−3.2 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.040.060.080.100.12 (c)
Absolute difference.
Figure 9: Solution of the viscous G-equation at T = 2 with d = 0 . d . One can findthat POD basis provides good approximations to the solution of the viscous G-equation when d is chosen from d = 0 .
01 to d = 0 .
1. The closer d is to d , the smaller the error will be. The factthat the pre-computed POD basis can be used to compute solutions associated with differentparameters with high accuracy is of great importance in practical computations. d 0.01 0.02 0.03 0.04 0.05Rela. error, Recovered d 0.06 0.07 0.08 0.09 0.1Rela. error, Recovered Table 3: The relative errors between POD solutions and reference solutions. S T Since the POD method is very efficient to approximately solve the viscous G-equation, we applyit to compute the turbulent flame speed S T , which is defined in the Eq.(10). In this experiment,we compare the numerical results of the turbulent flame speeds obtained by two methods. Firstof all, for each parameter d , we compute the solution snapshots from T = 0 to T = 1 to extractthe first 6 POD basis functions. Then, with the POD basis, we solve the viscous G-equationfor a much longer time, from T = 0 to T = 8, to obtain the approximated burned area A ( t ).Finally, we show the growth rate of the A ( t ) over time t in Fig.10. We find that the resultsobtained by our method agree well with the reference method. These results show that thePOD basis can be used to track the long-time evolution of the turbulent flame speed. Werefer the interested reader to the theoretical results on the long-time behavior of solutions toHamilton-Jacobi equations or quasi-linear parabolic equations; see e.g. [32, 5, 40]. Our resultsprovide an interesting numerical illustration of these theoretical results. A ( t ) speed by PODspeed by finite difference (a) d = 0 . A ( t ) speed by PODspeed by finite difference (b) d = 0 . A ( t ) speed by PODspeed by finite difference (c) d = 0 . Figure 10: Numerical results of the turbulent flame speed obtained by different methods.
In Table 4 we compare the CPU time of two methods when solving the viscous G-equation fromtime T = 0 to T = 1 with different d . The unit of the computational time is second. One cansee that computational cost of the POD method is far less than the finite difference method,20hich shows that model reduction method is very useful in many applications. Besides, onecan find that the CPU time of the POD method slightly increases when d decreases. This isdue to the fact that we are using more POD basis for small d in order to achieve the same PODtruncation error threshold e P OD = 0 . f k in Eq.(22)), which may deteriorate the performance of the POD method formore challenging problems, such as 3D problems or convection-dominated problems. However,one can apply the discrete empirical interpolation method [10] to maintain the efficiency of thePOD method. d 0.01 0.02 0.03 0.04 0.05POD method Reference method d 0.06 0.07 0.08 0.09 0.1POD method
Reference method
Table 4: CPU time of two methods in solving the viscous G-equation from T = 0 to 1 with different d . We consider a one-parameter family of time-dependent periodic cellular flow V ( x , t ) = (cid:0) V , V (cid:1) = A (cid:0) cos(2 πy ) , cos(2 πx ) (cid:1) + Aθ cos(2 πt ) (cid:0) (sin(2 πy ) , sin(2 πx ) (cid:1) , (49)where A is the amplitude of the velocity field. The first component of (49) is a steady cellularflow, while the second term of (49) is a time-periodic perturbation controlled by θ > θ increases. In theexperiments reported in Fig.11, we fix A = 4, θ = 1 . d ’s. For each d , the POD basis is extracted from the reference solution from time T = 0 to T = 1. Here, thenumbers of POD basis functions are r = 11 , ,
27 for d = 0 . , . , .
01, respectively. Thenusing the POD basis, we compute the POD solution from time T = 0 to T = 4. The turbulentflame speeds A ( t ) /t obtained by the POD method are compared with the reference solutions.The relative errors for the burning speed are basically negligible. These results justify thatthe POD method is still efficient in solving viscous G-equation with the time-periodic cellularflow. Compared with the results in Fig.10, the turbulent flame speed in Fig.11 has differentpatterns, which is caused by the mixing and chaotic features of the fluid velocity (49).We should point out that the strategy to choose the snapshots is a crucial point in the PODmethod. More generally, the POD basis obtained by minimizing the projection error over aset of snapshots chosen a priori , e.g. (20) may not be able to capture the whole dynamicsof equations with time-periodic or chaotic flows. Therefore, we propose an adaptive strategyto enrich the POD basis functions, which was also used in [11]. We provide the step-by-stepalgorithm of the adaptive strategy in the appendix 7.4. We consider the viscous G-equationwith the time-dependent velocity field (49) to demonstrate the main idea, where A = 4 . θ = 1 . d = 0 .
1. 21 .0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0t0.81.01.21.41.61.82.0 A ( t ) / t speed by PODspeed by finite difference (a) d = 0 . A ( t ) / t speed by PODspeed by finite difference (b) d = 0 . A ( t ) / t speed by PODspeed by finite difference (c) d = 0 . Figure 11: Turbulent flame speeds obtained by different methods for the time-dependent periodic cellular flow.
We first compute the snapshots from time T = 0 to T = 0 . r = 7 basis functions, which arecalled the initial POD basis functions. Notice that the time period of the velocity field (49) isone. Thus, the initial POD basis functions may not be able to capture the whole time-varyingdynamics of the viscous G-equation.In our adaptive strategy, we will check the approximation ability of the current POD basisafter every ∆ T time period. At each checking time, we use the current POD solution as initialdata and apply the finite difference scheme to compute a short period of time, say N ∆ t , toobtain N new snapshots. Here, ∆ t is the time step in the finite difference scheme and N > N ∆ t (cid:28) ∆ T . Then, we project the N snapshots onto the currentPOD basis and compute the norm of the projection error. If the norm of the projection erroris bigger than a prescribed threshold, we enrich the POD basis from the information in theprojection error.In this experiment, we choose ∆ T = 0 .
5, ∆ t = 0 .
001 and N = 50. We observe thatevery time, one or two basis functions will be added. At the end of the experiment, thereare r = 17 basis functions. We compare the performances of the fixed-basis POD methodwith the adaptive POD method, in terms of estimating the flame speed. Fig.12 clearly showsthat the adaptive strategy indeed improves the performance of the POD method and obtainsa more accurate estimation for the flame speed. More specifically, at the final time T = 4, therelative error between the speed calculated from the adaptive POD method and that from thefinite difference method is 1.0676%, while the relative error for the fixed-basis POD methodis 5.6648%. Meanwhile, the adaptive method still achieves high computational efficiency, asthe adaptive method is running the POD scheme for most time steps and only uses the finitedifference scheme for approximately 10% of the time steps. The CPU time for solving theequation from time t = 0 to t = 4 with the finite difference method is 2598.83 seconds, whilefor the adaptive POD method, it only takes 249.48 seconds.Notice that Fig.12 still shows a gap between the finite difference solution and the adaptivePOD solution. The reason is that we only check and enrich the POD basis every ∆ T time periodand we may lose a small amount of solution information between any two checking times. Howto choose ∆ T is important in the adaptivity of the POD method. To address this issue, oneneeds to obtain some a posteriori error estimate for the solution obtained by the POD method.22e will study this issue in our subsequent research. A ( t ) / t speed by AdaPODspeed by POD without Adaptationspeed by finite difference Figure 12: A comparison of the POD method with a fixed basis strategy and an adaptive strategy.
To further investigate the performance of the POD method, we solve the curvature G-equations(5) using the POD basis. Specifically, we consider the curvature G-equation with the periodicvelocity field (48) and planar initial condition. If we write G ( x , t ) = P · x + u ( x , t ), then u ( x , t )is spatially periodic and satisfies the following periodic initial value problem u t + V · ( P + ∇ u ) + S l | P + ∇ u | = dS l | P + ∇ u | (cid:18) ∇ · P + ∇ u | P + ∇ u | (cid:19) in T n × (0 , ∞ ) ,u ( x ,
0) = 0 on T n × { t = 0 } , (50)where T n = [0 , n . Again, we decompose u into the mean-free part ˆ u and the mean ¯ u . Thederivation is almost same as what we have done in Section 3.1. The evolution equations for ˆ u and ¯ u are as follows ˆ u t + V · ∇ ˆ u − dS l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12)(cid:18) ∇ · P + ∇ ˆ u | P + ∇ ˆ u | (cid:19) + (cid:82) T n dS l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12)(cid:18) ∇ · P + ∇ ˆ u | P + ∇ ˆ u | (cid:19) d x + S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) − (cid:82) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x + V · P = 0 in T n × (0 , ∞ )¯ u t = − (cid:82) T n S l (cid:12)(cid:12) P + ∇ ˆ u ( x , t ) (cid:12)(cid:12) d x + (cid:82) T n dS l | P + ∇ ˆ u | (cid:18) ∇ · P + ∇ ˆ u | P + ∇ ˆ u | (cid:19) d x on t ∈ (0 , ∞ )ˆ u ( x ,
0) = 0 on T n × { t = 0 } ¯ u (0) = 0 (51)The implementation of the POD method is the same as we have done for the viscous G-equation in Section 3. First of all, we use a reference numerical method to solve the curvatureG-equation (see Section 7.2) and construct the POD basis by using the mean-free part of the23olution snapshots. Then, we use the POD-based Galerkin method to solve the evolutionequation for ˆ u in Eq.(51). The corresponding weak formulation is (cid:104) ˆ u t , ψ (cid:105) + a (ˆ u, ψ ) + (cid:104) F (ˆ u ) , ψ (cid:105) = (cid:104)− P · V , ψ (cid:105) , ∀ ψ ∈ H, (52)where the bilinear form a (ˆ u, ψ ) = (cid:82) T n ( V · ∇ ˆ u ) ψd x and the nonlinear part F (ˆ u ) = S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) − (cid:82) T n S l (cid:12)(cid:12) P + ∇ ˆ u (cid:12)(cid:12) d x − dS l | P + ∇ ˆ u | (cid:0) ∇ · P + ∇ ˆ u | P + ∇ ˆ u | (cid:1) + (cid:82) T n dS l | P + ∇ ˆ u | (cid:0) ∇ · P + ∇ ˆ u | P + ∇ ˆ u | (cid:1) d x . Finally, we solvea nonlinear equation system to compute the evolution equation for ˆ u . The evolution equationfor the mean ¯ u can be solved subsequently.In our numerical experiment, we use the POD method to solve (51) with the steady flow(48). We choose A = 4 . d = 0 . S l = 1, and P = (1 , r = 6 POD basisfunctions from the reference solution snapshots from time T = 0 to T = 1. Then, we computethe mean-free POD solution on the same time interval using the POD-based Galerkin method.In the meanwhile, we compute the mean ¯ u solution. Finally, we can recover the solution u .The relative error of the recovered solutions between the POD method and reference method is0.0220. Fig.13 shows the recovered solutions using the POD method and the reference method.We find that the POD method still performs well in solving the curvature G-equations.If we compare Fig.13 with Fig.3, we can see that with the same diffusion parameter d = 0 . d in the curvature G-equation will bring more challenges in the method design and convergenceanalysis, which will be studied in our future research. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−2.0−1.8−1.6−1.4−1.2−1.0−0.8 (a) Reference solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.0−2.0−1.8−1.6−1.4−1.2−1.0−0.8 (b)
POD solution. x0.0 0.2 0.4 0.6 0.8 1.0y 0.00.20.40.60.81.00.010.020.030.040.050.06 (c)
Absolute difference.
Figure 13: Solution of the curvature G-equation at T = 1 with d = 0 .
6. Conclusion
We have proposed an efficient model reduction method to solve viscous G-equations, whichhave been very popular field models in combustion and physics literature for studying turbulentflame propagation. We constructed the POD basis based on learning the solution informationfrom the snapshots. Then, we applied the Galerkin project method to solve the viscous G-equation and smooth solutions of curvature G-equation by using the POD basis. We providedrigorous error analysis for our numerical methods based on a decomposition strategy, where wedecomposed the solution into a mean part and a mean-free part. We showed through numerical24xperiments that our methods can accurately compute the various G-equations with significantcomputational savings. In addition, we found that the POD basis allows us to compute long-time solution of the various G-equations. Thus, we can compute the corresponding turbulentflame speeds in cellular flows. In our future work, we plan to study turbulent flame speedsof G-equations in three dimensional spatially or spatiotemporally periodic vortical flows. Inaddition, we will further investigate the adaptivity for the POD method and using the PODmethod to solve curvature G-equations.
7. Appendix
We first apply the finite difference scheme proposed in [27] to solve the viscous G-equation fromtime 0 to T seconds on the domain D = [0 , × [0 ,
1] to get the snapshots. Specifically, we employa higher-order HJ-WENO scheme and TVD-RK scheme in spatial and time discretization,respectively; see [34, 16, 20, 33] for more details of these schemes.For a small d , the viscous G-equation is convection dominated and it should be treated likea hyperbolic equation. The forward Euler time discretization is given by G n +1 − G n ∆ t + H n ( G − x , G + x , G − y , G + y ) − dS l ∆ G n = 0 , (53)where G − i and G + i denote the left and right discretization of G i in the WENO5 scheme [20]. H is a consistent and monotone numerical Hamiltonian. Then, we have H ( G − x , G + x , G − y , G + y ) = V G velx + V G vely + S l (cid:113) ( G norx ) + ( G nory ) , (54)where the upwinding scheme and the Godunov scheme are applied for the convection term andthe nonlinear term separately [33]. G velx = (cid:40) G − x if V > ,G + x if V < , (55) G vely = (cid:40) G − y if V > ,G + y if V < , (56)( G norx ) = ( G − x ) if V > S l , max(max( G − x , , min( G + x , ) if | V | ≤ S l , ( G + x ) if V < − S l , (57)( G nory ) = ( G − y ) if V > S l , max(max( G − y , , min( G + y , ) if | V | ≤ S l , ( G + y ) if V < − S l . (58)For the diffusion term, we apply the central difference. For the time discretization, we applythe RK3 scheme [16]. The CFL condition in this case is∆ t (cid:18) S l + | V | ∆ x + S l + | V | ∆ y (cid:19) < . (59)25hen d is large, the time step size for the forward Euler scheme is very small, i.e., ∆ t = O (cid:0) (∆ x ) + (∆ y ) ) (cid:1) . To alleviate the stringent time step restriction, we introduce the followingsemi-implicit scheme: G n +1 − G n dt + V · ∇ G n +1 + S l |∇ G n | = dS l ∆ G n +1 , (60)where the convection and diffusion terms are discretized by the central difference, and thenormal direction term is discretized by the Godunov and WENO5 scheme. In this case, theCFL condition is ∆ t (cid:18) S l ∆ x + S l ∆ y (cid:19) < . (61)Numerical results in [27] show that the proposed method is efficient in solving G-equations.The limitation is that the computational cost becomes large when one needs to choose a finemesh to discretize the problem. We also apply the finite difference scheme to solve the curvature G-equation (5), which wasproposed in [27]. The spatial and time discretization are the same as for the viscous G-equationsin Section 7.1. The curvature term is given by |∇ G | (cid:18) ∇ · ∇ G |∇ G | (cid:19) = G y G xx − G x G y G xy + G x G yy G x + G y , (62)where each partial derivative is approximated by the central difference. The computation ofthe numerical Hamiltonian is same as before; see Eqns.(54)-(58). For the time discretization,we apply the RK3 scheme. The time step restriction is∆ t (cid:18) S l + | V | ∆ x + S l + | V | ∆ y + 2 S l d (∆ x ) + 2 S l d (∆ y ) (cid:19) < . (63)When d is large, in order to avoid a small time step size ∆ t = O ((∆ x ) ), the curvature termis decomposed as |∇ G | (cid:18) ∇ · ∇ G |∇ G | (cid:19) = ∆ G − ∆ ∞ G = G xx + G yy − G x G xx + G x G y G xy + G y G yy G x + G y . (64)If we apply the backward Euler scheme on ∆ G and forward Euler scheme on ∆ ∞ G , then wehave the following semi-implicit time discretization scheme G n +1 − G n dt + V · ∇ G n +1 + S l |∇ G n | = dS l (∆ G n +1 − ∆ ∞ G n ) . (65)The convection and diffusion terms are discretized by the central difference, and the normaldirection term is discretized by the Godunov and WENO5 scheme. In this case, the CFLcondition is same as (59). 26 .3. Model reduction using the POD method Let X be a Hilbert space equipped with the inner product (cid:104)· , ·(cid:105) X and u ( · , t ) ∈ X , t ∈ [0 , T ]be the solution of a dynamic system. In practice, we approximate the space X by a linearfinite dimensional space V with dim ( V ) = N dof , where N dof represents the degree of freedomof the solution space and N dof can be extremely large for high-dimensional problem (considerthe finite element method or finite difference method as examples). Given a set of snapshot ofsolutions, denoted as { u ( · , t ) , u ( · , t ) , · · · , u ( · , t m ) } , where t , · · · , t m ∈ [0 , T ] are different timeinstances and m is the number of the solution snapshots.The POD method aims to build a set of orthonormal basis { ψ ( · ) , ψ ( · ) , · · · , ψ r ( · ) } with r ≤ min( m, N dof ) that optimally approximates the solution snapshots. More specifically, thePOD method constructs the basis by solving the following optimization problemmin ψ , ··· ,ψ r ∈ X m m (cid:88) i =1 (cid:13)(cid:13)(cid:13) u ( · , t i ) − r (cid:88) j =1 (cid:104) u ( · , t i ) , ψ j ( · ) (cid:105) X ψ j ( · ) (cid:13)(cid:13)(cid:13) X , s.t. (cid:104) ψ i , ψ j (cid:105) X = δ ij . (66)In order to solve Eq.(66), we consider the eigenvalue problem Kv = λv, where K ∈ R m × m and K ij = m (cid:104) u ( · , t i ) , u ( · , t j ) (cid:105) X is the snapshot correlation matrix. Let v k , k = 1 , · · · , n be the eigenvectors and λ ≥ λ ≥ · · · ≥ λ m > ψ k ( · ) = 1 √ λ k m (cid:88) j =1 ( v k ) j u ( · , t j ) , ≤ k ≤ r. (67)It can also be shown that the following error formula holds1 m m (cid:88) i =1 (cid:107) u ( · , t i ) − r (cid:88) j =1 (cid:104) u ( · , t i ) , ψ j ( · ) (cid:105) X ψ j ( · ) (cid:107) X = m (cid:88) l = r +1 λ l . (68)Finally, let S r denote the r -dimensional space spanned by { ψ ( · ) , ψ ( · ) , · · · , ψ r ( · ) } . In this subsection, we formalize the proposed adaptive strategy for enriching the POD basis inAlgorithm 1, where the notations have been defined before.
Acknowledgement
The research of J. Xin is partially supported by NSF grants DMS-1211179, DMS-1522383 andIIS-1632935. The research of Z. Zhang is supported by Hong Kong RGC grants (Projects27300616, 17300817, and 17300318), National Natural Science Foundation of China (Project11601457), Basic Research Programme (JCYJ20180307151603959) of The Science, Technologyand Innovation Commission of Shenzhen Municipality, Seed Funding Programme for BasicResearch (HKU), and an RAE Improvement Fund from the Faculty of Science (HKU).27 lgorithm 1 Adaptive strategy for enriching the POD basis Input : POD basis S = { ϕ i } , error threshold (cid:15) >
0, computational time T , period ofchecking time ∆ T , time step of POD or finite difference method ∆ t , N T = (cid:100) T / ∆ t (cid:101) , N check = (cid:100) ∆ T / ∆ t (cid:101) , and snapshots number N (where N ∆ t (cid:28) ∆ T ). Parameters in the G-equations.Solution at time t = 0: U . for i = 1 : N T do if i > i mod N check ) = 0 then Set U = ∅ . for j = 1 : N do Apply one step of the finite-difference scheme proposed by [27] (see Appendix 7.1and 7.2) to solve the equation from time ( i + j − t to time ( i + j )∆ t . Denote the finite-difference solution at time ∆ t ( i + j ) by u j . U = U ∪ { u j } . end for Let P S be the projection onto the space spanned by S . Compute P S U = { P S u j } . Compute the POD basis S (cid:48) for P S U , such that the approximation error (68) is lessthan (cid:15) (see Appendix 7.3). S = S ∪ S (cid:48) . end if Apply one step of the backward POD-based scheme (22), with basis S , to solve theequation from time ( i − t to time i ∆ t . Denote the solution at time i ∆ t by U i . end for Output : { U t } References [1]
A. Alla and M. Falcone , An adaptive POD approximation method for the controlof advection-diffusion equations , in Control and Optimization with PDE Constraints,Springer, 2013, pp. 1–17.[2]
A. Alla and M. Falcone , A time-adaptive POD method for optimal control problems ,IFAC Proceedings Volumes, 46 (2013), pp. 245–250.[3]
A. Alla and J. Kutz , Randomized model order reduction , Advances in ComputationalMathematics, 45 (2019), pp. 1251–1271.[4]
A. Alla and L. Saluzzi , A HJB-POD approach for the control of nonlinear PDEs ona tree structure , Applied Numerical Mathematics, 155 (2020), pp. 192–207.[5]
G. Barles and P. Souganidis , Space-time periodic solutions and long-time behavior ofsolutions to quasi-linear parabolic equations , SIAM Journal on Mathematical Analysis, 32(2001), pp. 1311–1323.[6]
P. Benner, S. Gugercin, and K. Willcox , A survey of projection-based model re-duction methods for parametric dynamical systems , SIAM Review, 57 (2015), pp. 483–531.287]
G. Berkooz, P. Holmes, and J. Lumley , The proper orthogonal decomposition in theanalysis of turbulent flows , Annual review of fluid mechanics, 25 (1993), pp. 539–575.[8]
J. Borggaard, T. Iliescu, and Z. Wang , Artificial viscosity proper orthogonal de-composition , Mathematical and Computer Modelling, 53 (2011), pp. 269–279.[9]
K. Brakke , The Motion of a Surface by Its Mean Curvature. , Princeton University Press,2015.[10]
S. Chaturantabut and D. Sorensen , Nonlinear model reduction via discrete empiricalinterpolation , SIAM Journal on Scientific Computing, 32 (2010), pp. 2737–2764.[11]
M. Cheng, T. Hou, and Z. Zhang , A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations II: adaptivity and generalizations , Jour-nal of Computational Physics, 242 (2013), pp. 753–776.[12]
B. Cockburn, G. Karniadakis, and C. Shu , The development of discontinuousGalerkin methods , in Discontinuous Galerkin Methods, Springer, 2000, pp. 3–50.[13]
P. Embid, A. Majda, and P. Souganidis , Comparison of turbulent flame speeds fromcomplete averaging and the G-equation , Physics of Fluids, 7 (1995), pp. 2052–2060.[14]
L. Evans and J. Spruck , Motion of level sets by mean curvature. I , Journal of Differ-ential Geometry, 33 (1991), pp. 635–681.[15]
M. Falcone and R. Ferretti , Semi-Lagrangian approximation schemes for linear andHamilton–Jacobi equations , SIAM, 2013.[16]
S. Gottlieb and C. Shu , Total variation diminishing Runge-Kutta schemes , Mathe-matics of computation of the American Mathematical Society, 67 (1998), pp. 73–85.[17]
C. Gr¨aßle and M. Hinze , POD reduced-order modeling for evolution equations utiliz-ing arbitrary finite element discretizations , Advances in Computational Mathematics, 44(2018), pp. 1941–1978.[18]
S. Gugercin and A. C. Antoulas , A survey of model reduction by balanced truncationand some new results , International Journal of Control, 77 (2004), pp. 748–766.[19]
J. Hesthaven, G. Rozza, and B. Stamm , Certified reduced basis methods forparametrized partial differential equations , Springer, 2016.[20]
G. Jiang and D. Peng , Weighted ENO schemes for Hamilton–Jacobi equations , SIAMJournal on Scientific Computing, 21 (2000), pp. 2126–2143.[21]
S. Jin and Z. Xin , Numerical passage from Systems of Conservation Laws to Hamilton–Jacobi Equations, and Relaxation Schemes , SIAM Journal on Numerical Analysis, 35(1998), pp. 2385–2404. 2922]
D. Kalise and A. Kr¨oner , Reduced-order minimum time control of advection-reaction-diffusion systems via dynamic programming. , in 21st International Symposium on Mathe-matical Theory of Networks and Systems, 2014, pp. 1196–1202.[23]
D. Kalise and K. Kunisch , Polynomial approximation of high-dimensional Hamilton–Jacobi–Bellman equations and applications to feedback control of semilinear parabolicPDEs , SIAM Journal on Scientific Computing, 40 (2018), pp. A629–A652.[24]
K. Kunisch and S. Volkwein , Galerkin proper orthogonal decomposition methods forparabolic problems , Numerische Mathematik, 90 (2001), pp. 117–148.[25]
K. Kunisch, S. Volkwein, and L. Xie , HJB-POD-based feedback design for the optimalcontrol of evolution problems , SIAM Journal on Applied Dynamical Systems, 3 (2004),pp. 701–722.[26]
A. Kurganov and E. Tadmor , New high-resolution semi-discrete central schemes forHamilton–Jacobi equations , Journal of Computational Physics, 160 (2000), pp. 720–742.[27]
Y. Liu, J. Xin, and Y. Yu , A numerical study of turbulent flame speeds of curvatureand strain G-equations in cellular flows , Physica D: Nonlinear Phenomena, 243 (2013),pp. 20–31.[28]
J. Lyu, J. Xin, and Y. Yu , Computing residual diffusivity by adaptive basis learning viaspectral method , Numerical Mathematics: Theory, Methods and Applications, 10 (2017),pp. 351–372.[29] ,
Computing residual diffusivity by adaptive basis learning via super-resolution deepneural networks , in International Conference on Computer Science, Applied Mathematicsand Applications, Springer, 2019, pp. 279–290.[30]
G. Markstein , Nonsteady flame propagation , Pergamon Press, Oxford, 1964.[31]
B. Matkowsky and G. Sivashinsky , An asymptotic derivation of two models in flametheory associated with the constant density approximation , SIAM Journal on Applied Math-ematics, 37 (1979), pp. 686–699.[32]
G. Namah and J. Roquejoffre , Convergence to periodic fronts in a class of semilinearparabolic equations , Nonlinear Differential Equations and Applications NoDEA, 4 (1997),pp. 521–536.[33]
S. Osher and R. Fedkiw , Level set methods and dynamic implicit surfaces , vol. 153,Springer Science & Business Media, 2006.[34]
S. Osher and J. Sethian , Fronts propagating with curvature-dependent speed: al-gorithms based on Hamilton-Jacobi formulations , Journal of computational physics, 79(1988), pp. 12–49.[35]
B. Peherstorfer and K. Willcox , Dynamic data-driven reduced-order models , Com-puter Methods in Applied Mechanics and Engineering, 291 (2015), pp. 21–41.3036]
N. Peters , Turbulent combustion , Cambridge university press, 2000.[37]
A. Quarteroni, A. Manzoni, and F. Negri , Reduced basis methods for partial dif-ferential equations: an introduction , vol. 92, Springer, 2015.[38]
M.-L. Rap´un, F. Terragni, and J. Vega , Adaptive POD-based low-dimensional mod-eling supported by residual estimates , International Journal for Numerical Methods in En-gineering, 104 (2015), pp. 844–868.[39]
M.-L. Rap´un and J. M. Vega , Reduced order models based on local pod plus galerkinprojection , Journal of Computational Physics, 229 (2010), pp. 3046–3063.[40]
J. Roquejoffre , Convergence to steady states or periodic solutions in a class ofHamilton–Jacobi equations , Journal de math´ematiques pures et appliqu´ees, 80 (2001),pp. 85–104.[41]
C. Shu and S. Osher , Efficient implementation of essentially non-oscillatory shock-capturing schemes , Journal of computational physics, 77 (1988), pp. 439–471.[42]
L. Sirovich , Turbulence and the dynamics of coherent structures. I. Coherent structures ,Quarterly of applied mathematics, 45 (1987), pp. 561–571.[43]
V. Thom´ee , Galerkin finite element methods for parabolic problems , vol. 1054, Springer,1984.[44]
S. Volkwein , Proper orthogonal decomposition: Theory and reduced-order modelling ,Lecture Notes, University of Konstanz, 4 (2013).[45]