An ensemble-Proper OrthogonalDecomposition method for theNonstationary Navier-Stokes EquationsSupported by the US Air Force Office of Scientific Research grant FA9550-15-1-0001 and US Department of Energy Office of Science grants DE-SC0009324 and DE-SC0010678.

An ensemble-Proper Orthogonal
Decomposition method for the
Nonstationary Navier-Stokes Equationsthanks: Supported by the US Air Force Office of Scientific Research grant FA9550-15-1-0001 and US Department of Energy Office of Science grants DE-SC0009324 and DE-SC0010678.

Max Gunzburger, Nan Jiang, and Michael Schneier Department of Scientific Computing, Florida State University, Tallahassee, FL 32306-4120 (mgunzburger@fsu.edu, njiang@fsu.edu, mschneier89@gmail.com).
Abstract

The definition of partial differential equation (PDE) models usually involves a set of parameters whose values may vary over a wide range. The solution of even a single set of parameter values may be quite expensive. In many cases, e.g., optimization, control, uncertainty quantification, and other settings, solutions are needed for many sets of parameter values. We consider the case of the time-dependent Navier-Stokes equations for which a recently developed ensemble-based method allows for the efficient determination of the multiple solutions corresponding to many parameter sets. The method uses the average of the multiple solutions at any time step to define a linear set of equations that determines the solutions at the next time step. To significantly further reduce the costs of determining multiple solutions of the Navier-Stokes equations, we incorporate a proper orthogonal decomposition (POD) reduced-order model into the ensemble-based method. The stability and convergence results for the ensemble-based method are extended to the ensemble-POD approach. Numerical experiments are provided that illustrate the accuracy and efficiency of computations determined using the new approach.

Key words. Ensemble methods, proper orthogonal decomposition, reduced-order models, Navier-Stokes equations.

1 Introduction

Computing an ensemble of solutions of fluid flow equations for a set of parameters or initial/boundary conditions for, e.g., quantifying uncertainty or sensitivity analyses or to make predictions, is a common procedure in many engineering and geophysical applications. One common problem faced in these calculations is the excessive cost in terms of both storage and computing time. Thanks to recent rapid advances in parallel computing as well as intensive research in ensemble-based data assimilation, it is now possible, in certain settings, to obtain reliable ensemble predictions using only a small set of realizations. Successful methods that are currently used to generate perturbations in initial conditions include the Bred-vector method, [33], the singular vector method, [8], and the ensemble transform Kalman filter, [7]. Despite all these efforts, the current level of available computing power is still insufficient to perform high-accuracy ensemble computations for applications that deal with large spatial scales such as numerical weather prediction. In such applications, spatial resolution is often sacrificed to reduce the total computational time. For these reasons the development of efficient methods that allow for fast calculation of flow ensembles at a sufficiently fine spatial resolution is of great practical interest and significance.

Only recently, a first step was taken in [25, 26] where a new algorithm was proposed for computing an ensemble of solutions of the time-dependent Navier-Stokes equations (NSE) with different initial condition and/or body forces. At each time step, the new method employs the same coefficient matrix for all ensemble members. This reduces the problem of solving multiple linear systems to solving one linear system with multiple right-hand sides. There have been many studies devoted to this type of linear algebra problem and efficient iterative methods have been developed to significantly save both storage and computing time, e.g., block CG [13], block QMR [14], and block GMRES [15]. Even for some direct methods, such as the simple LU factorization, one can save considerable computing cost.

Because the main goal of the ensemble algorithm is computational efficiency, it is natural to consider using reduced-order modeling (ROM) techniques to further reduce the computational cost. Specifically, we consider the proper orthogonal decomposition (POD) method which has been extensively used in the engineering community since it was introduced in [30] to extract energetically coherent structures from turbulent velocity fields. POD provides an optimally ordered, orthonormal basis in the least-squares sense, for given sets of experimental or computational data. The reduced order model is then obtained by truncating the optimal basis.

Research on POD and its application to the unsteady NSE has been and remains a highly active field. Recent works improving upon POD have dealt with the combination of Galerkin strategies with POD [6, 11], stabilization techniques [2, 10, 32], and regularized/large eddy simulation POD models for turbulent flows [35, 36].

In this paper, we study a Galerkin proper orthogonal decomposition (POD-G-ROM) based ensemble algorithm for approximating solutions of the NSE. Accordingly, our aim in this paper is to develop and demonstrate a procedure for the rapid solution of multiple solutions of the NSE, requiring only the solution of one reduced linear system with multiple right-hand sides at each time step.

1.1 Previous works on ensemble algorithms

The ensemble method given in [25] is first-order accurate in time and requires a CFL-like time step condition to ensure stability and convergence. Two ensemble eddy viscosity numerical regularizations are studied in [26] to relax the time step restriction. These two methods utilized the available ensemble data to parametrize the eddy viscosity based on a direct calculation of the kinetic energy in fluctuations without further modeling. They both give the same parametrization for each ensemble member and thus preserve the efficiency of the ensemble algorithm. The extension of the ensemble method to higher-order accurate ensemble time discretization is nontrivial. For instance, the method is not extensible to the most commonly used Crank-Nicolson scheme. Making use of a special combination of a second-order in time backward difference formula and an explicit second-order Adams-Bashforth treatment of the nonlinear term, a second-order accurate in time ensemble method was developed in [22].Another second-order ensemble method with improved accuracy is presented in [23]. The ensemble algorithm was further used in [24] to model turbulence. By analyzing the evolution of the model variance, it was proved that the proposed ensemble based turbulence model converges to statistical equilibrium, which is a desired property of turbulence models.

2 Notation and preliminaries

Let , , denote an open regular domain with boundary and let denote a time interval. Consider Navier-Stokes equations on a bounded domain, each subject to the no-slip boundary condition, and driven by different initial conditions and body force densities , i.e., for , we have

\hb@xt@.01(2.1)

where denotes the given constant kinematic viscosity of the fluid and and respectively denote the velocity and pressure of the fluid flow.

We denote by and the norm and inner product, and denote by and the norms and the Sobolev norms respectively. The space is the Sobolev space , equipped with norm . The space denotes the dual space of bounded linear functionals defined on ; this space is endowed with the norm

The solutions spaces for the velocity and for the pressure are respectively defined as

Let denote the usual norm. For a function that is well defined on we define the norms

The subspace of consisting of weakly divergence free functions is defined as

A weak formulation of (LABEL:eq:NSE) is given as follows: for , find and that, for almost all , satisfy

\hb@xt@.01(2.2)

Conforming velocity and pressure finite element spaces based on a regular triangulation of having maximum triangle diameter are respectively denoted by

We assume that the pair of spaces satisfy the discrete inf-sup (or ) condition required for stability of finite element approximation; we also assume that the finite element spaces satisfy the approximation properties

for a constant having value independent of . The total number of finite element degrees of freedom is given by . A concrete example for which the stability condition approximation estimates are known to hold is the family of Taylor-Hood -, , element pairs [16, 17]. For the most commonly used Taylor-Hood element pair based on a tetrahedral grid, is roughly equal to three times the number of vertices plus twice the number of edges.

Further, in this paper we will need to solve the NSE (LABEL:eq:NSE) using a second order time stepping scheme (e.g., Crank Nicolson). We will assume the FE approximations satisfy the following error estimates:

\hb@xt@.01(2.3)
\hb@xt@.01(2.4)

The subspace of consisting of discretely divergence free functions is defined as

Note that in most cases, and for the Taylor-Hood element pair in particular, , i.e., discretly divergence free functions are not divergence free.

As is common to do, we define the explicitly skew-symmetric trilinear form introduced by Temam given by

This form satisfies the bounds, [28]

\hb@xt@.01(2.5)
\hb@xt@.01(2.6)

Moreover, we have that for all so that we may replace the nonlinear term in (LABEL:wfwf) by . The advantage garnered through the use of compared to is that for all whereas only if .

Definition 2.1

Let , , where , denote a partition of the interval . For and , let . Then, the ensemble mean is defined, for , by

For , let denote approximations, e.g., interpolants or projections, of the initial conditions . Then, the full space-time discretization of (LABEL:eq:NSE), or more precisely of (LABEL:wfwf), we consider is given as follows: given, for , and for almost every , find, for and for , and satisfying

\hb@xt@.01(2.7)

We refer to this discretization as En-full-FE indicating that we are referring to an ensemble-based discretization of (LABEL:wfwf) using a high-dimensional finite element space. This ensemble-based discretization of the NSE is noteworthy because the system (LABEL:En-full-FE) is not only linear in the unknown functions and , but because of the use of ensembles, we also have that the coefficient matrix associated with (LABEL:En-full-FE) is independent of , i.e., at each time step, all members of the ensemble can be determined from linear algebraic systems all of which have the same coefficient matrix. On the other hand, the linear system can be very large because in practice can be very large. This observation, in fact, motivates interest in building reduced-order discretizations of the NSE.

Because and are assumed to satisfy the condition, (LABEL:En-full-FE) can be more compactly expressed as follows: given, for , and for almost every , find, for and for , satisfying

\hb@xt@.01(2.8)

Note that in general it is a difficult matter to construct a basis for the space so that in practice, one still works with (LABEL:En-full-FE). We introduce the reduced system (LABEL:eq:_conv) so as to facilitate the analyses given in later sections.

3 Proper orthogonal decomposition (POD) reduced-order modeling

The POD model reduction scheme can be split into two main stages: an offline portion and an online portion. In the offline portion, one collects into what is known as a snapshot set the solution of a partial differential equation (PDE), or more precisely, of a discrete approximation to that solution, for a number of different input functions and/or evaluated at several time instants. The snapshot set is hopefully generated in such a way that it is representative of the behavior of the exact solution. The snapshot set is then used to generate a POD basis, hopefully of much smaller cardinality compared to that of the full finite element space, that provides a good approximation to the data present in the snapshot set itself. In the online stage, the POD basis is used to generate approximate solutions of the PDE for other input functions; ideally these will be accurate approximations achieved much more cheaply compared to the use of a standard method such as a standard finite element method.

In the rest of this section, we delve into further detail about the generation of the snapshot set, the construction of the POD basis in a finite element setting, and how the POD basis can be used to construct a reduced-order model for the NSE in the ensemble framework. This section will focus on the framework specific to this paper; for more detailed presentations about POD, see, e.g., [12, 18, 19, 34].

3.1 Snapshot set generation

The offline portion of the algorithm begins with the construction of the snapshot set which consists of the solution of the PDE for a number of different input functions and/or evaluated at several different time instants. Given a positive integer , let denote a uniform partition of the time interval . Note that this partition is usually much coarser than the partition of into intervals, introduced in Definition LABEL:def21, which is used to discretize the PDE, i.e., we have . We first define the set of snapshots corresponding to exact solutions of the weak form of the NSE (LABEL:wfwf). For , we select different initial conditions and denote by the exact velocity field satisfying (LABEL:wfwf), evaluated at , , which corresponds to the initial condition . Then, the space spanned by the so obtained snapshots is defined as

\hb@xt@.01(3.1)

In the same manner, we can construct a set of snapshots , , , of finite element approximations of the velocity solution determined from a standard finite element discretization of (LABEL:wfwf). Note that one could also determine, at lesser cost but with some loss of accuracy, the snapshots from the ensemble-based discretization (LABEL:En-full-FE). We can then also define the space spanned by the discrete snapshots as

\hb@xt@.01(3.2)

Note that . The snapshots are finite element solutions so the span of the snapshots is a subset of the finite element space . Additionally, it is important to note that by construction, the snapshots satisfy the discrete continuity equation so that the span of the snapshots is indeed a subspace of the discretly divergence free subspace .

If we denote by the vector of coefficients corresponding to the finite element function . With , we may also define the snapshot matrix as

i.e., the columns of are the finite element coefficient vectors of the discrete snapshots.

To construct a reduced basis that results in accurate approximations, the snapshot set must contain sufficient information about the dynamics of the solution of the PDE. In our context, this requires one to not only take a sufficient number of snapshots with respect to time, but also to select a set of initial conditions that generate a set of solutions that is representative of the possible dynamics one may encounter when using other initial conditions. In the POD framework for the NSE, the literature on selecting this set is limited. One of the few algorithms which has been explored in the ensemble framework is the previously mentioned Bred-vectors algorithm given in [33]. Further exploration of this and other approaches for the selection of initial conditions is a subject for future research.

3.2 Construction of the POD basis

Using the set of discrete snapshots, we next construct the POD basis . We define the POD function space as

There are a number of equivalent ways in which one may characterize the problem of determining ; for a full discussion see [9, Section 2]. For example, the POD basis construction problem can be defined as follows: determine an orthonormal basis for such that for all , solves the following constrained minimization problem

\hb@xt@.01(3.3)

where is the Kronecker delta and the minimization is with respect to all orthonormal bases for . We note that by defining our basis in this manner we elect to view the snapshots as finite element functions as opposed to finite element coefficient vectors.

Define the correlation matrix , where denotes the Gram matrix corresponding to full finite element space. Then, the problem (LABEL:Min) is equivalent to determine the dominant eigenpairs satifying

\hb@xt@.01(3.4)

where denotes the Euclidean norm of a vector. The finite element coefficient vectors corresponding to the POD basis functions are then given by

\hb@xt@.01(3.5)

Alternatively, we can let , and define so that and then determine the singular value decomposition of the modified snapshot matrix ; the vectors , are then given as the first left singular vectors of which correspond to the first singular values .

3.3 POD reduced-order modeling

We next illustrate how a POD basis is used to construct a reduced-order model for the NSE within the ensemble framework. The discretized system that defines the POD approximation mimics that for the full finite element approximation, except that now we seek an approximation in the POD space having the basis . Specifically, for , we define the POD approximate initial conditions as and then pose the following problem: given , for and for , find satisfying

\hb@xt@.01(3.6)

We refer to this discretization as En-POD indicating that we are referring to an ensemble-based discretization of (2.2) using a low-dimensional POD space. Note that because , i.e., the POD approximation is by construction discretely divergence free, the pressure term in the POD-discretized NSE (LABEL:En-POD-Weak) drops out and we are left with a system involving only the POD approximation to the velocity. One further point of emphasis is that the initial conditions used in (LABEL:En-POD-Weak) are different from the initial conditions used to construct the snapshot set, i.e., we use initial conditions to solve the full finite element system (LABEL:En-full-FE) to determine the snapshots, and now solve additional approximations of the NSE by solving the much smaller POD system (LABEL:En-POD-Weak).

As was the case for (LABEL:En-full-FE), the POD system (LABEL:En-POD-Weak) is linear in the unknown and the associated coefficient matrix does not depend on , i.e., it is the same for all realizations of the initial condition. On the other hand, (LABEL:En-POD-Weak) is a system of equations in unknowns whereas (LABEL:En-full-FE) involves equations in the same number of unknowns, where and denote the total number of POD and finite element degrees of freedom, respectively. Thus, if , solving (LABEL:En-POD-Weak) requires much less cost compared to solving (LABEL:En-full-FE). In this way the offline cost of constructing the POD basis can be amortized over many online solves using the much smaller POD system. We address the assembly costs related to (LABEL:En-POD-Weak) in Section LABEL:numexp.

4 Stability analysis of En-POD

We prove the conditional, nonlinear, long-time stability of solutions of (LABEL:En-POD-Weak).

The projection operator : is defined by

\hb@xt@.01(4.1)

Denote by the spectral norm for symmetric matrices and let denote the POD mass matrix with entries and denote the matrix with entries , . It is shown in [27] that

\hb@xt@.01(4.2)

As , we have the following lemma, see of [25, page 276] for proof.

Lemma 4.1

For any ,

Theorem 4.2

[Stability of En-POD] For and , let satisfy (LABEL:En-POD-Weak). Suppose the time-step condition

\hb@xt@.01(4.3)

holds. Then, for ,

\hb@xt@.01(4.4)

Proof. The proof is provided in Appendix LABEL:app1.     

Remark 4.3

In the time-step condition, the constant C is dependent on the shape of the domain and the mesh as a result of the use of inverse inequality in the proof. For a fixed mesh on a fixed domain, C is a generic constant that is independent of the time step , the solution and viscosity .

5 Error analysis of En-POD

We next provide an error analysis for En-POD solutions.

Lemma 5.1

[ norm of the error between snapshots and their projections onto the POD space] We have

\hb@xt@.01(5.1)

and thus for ,

\hb@xt@.01(5.2)

Proof. The proof of (LABEL:errL2_1) follows exactly the proof of [34, Theorem 3]; (LABEL:errL2_2) is then a direct consequence of (LABEL:errL2_1).     

Lemma 5.2

[ norm of the error between snapshots and their projections in the POD space] We have

\hb@xt@.01(5.3)

and thus for ,

\hb@xt@.01(5.4)

Proof. The proof of (LABEL:errH1_1) follows exactly the proof of [21, Lemma 3.2]; (LABEL:errH1_2) is then direct consequence of (LABEL:errH1_1).     

Lemma 5.3

[Error in the projection onto the POD space] Consider the partition used in Section LABEL:snapshot. For any , let . Then, the error in the projection onto the POD space satisfies the estimates

\hb@xt@.01(5.5)
\hb@xt@.01(5.6)

Proof. The proof is provided in Appendix LABEL:app2.     

To bound the error between the POD based approximations and the true solutions, we assume the following regularity for the true solutions and body forces:

We assume the following estimate is also valid as done in [21].

Assumption 5.4

Consider the partition used in Section LABEL:snapshot. For any , let . Then, the error in the projection onto the POD space satisfies the estimates

\hb@xt@.01(5.7)

Let be the error between the true solution and the POD approximation, then we have the following error estimates.

Theorem 5.5 (Error analysis of En-POD)

Consider the method (LABEL:eq:_conv) and the partition used in Section 3.1. Suppose that for any , the following conditions hold

\hb@xt@.01(5.8)

Then, for any , there is a positive constant such that

\hb@xt@.01(5.9)

Proof. The proof is provided in Appendix LABEL:app3.     

6 Numerical simulations

We investigate the efficacy of our algorithm via the numerical simulation of a flow between two offset circles [25]. Before we discuss the examples and the numerical results, we briefly discuss the computational costs associated with the En-POD algorithm and how they compare to those of the En-full-FE algorithm.

6.1 Computational costs

As stated in Section LABEL:PODsec, we can split the computational cost of our algorithm into offline and online portions. In the offline portion, we generate the snapshot matrix by solving the Navier-Stokes equations for perturbations. Using , we then generate a reduced basis to be used in our online calculations. It is fair to assume that the cost of creating the snapshot matrix will dominate the cost of generating the reduced basis associated with the eigenvalue problem (LABEL:eigProb), especially when we consider that there exist very efficient techniques [20] for determining the partial SVD of matrices.

Turning to the cost of solving the Navier-Stokes equation, the discrete systems that arise from a FEM discretization have been studied at great length. Whereas it is possible to use a nonlinear solver such as Newton’s method or a nonlinear multigrid iteration, these methods often suffer from a lack of robustness. Instead, it is more popular to linearize the system and then to use the Schur complement approach. This allows for the use of a linear multigrid solver or Krylov method such as GMRES to solve the problem. For full details, see, e.g., [31]. Unfortunately, there are a number of factors such as the mesh size, the value of the Reynolds number, and the choice of pre-conditioner which make it very difficult to precisely estimate how quickly these methods converge.

Estimating the online cost of the En-POD method, however, is much easier. Because the POD discrete system is small and dense and the ensemble method has right-hand sides, the most efficient way to solve this problem is, at each time step, to do a single LU factorization and a backsolve for each right-hand side. Denoting again by the cardinality of the reduced basis, the online cost of the En-POD method is

\hb@xt@.01(6.1)

We note that this process is highly parallelizable. For example, if we have access to total processors, then we can remove the factor in the second term.

It is important to note that the assembly of the low-dimensional reduced basis system requires manipulations involving the reduced basis which, as we have seen, are finite element functions so that, in general, that assembly involves computational costs that depend on the dimension of the finite element space. Thus, naive implementations of a reduced basis method involve assembly costs that are substantially greater than solving costs and which, given the availability of very efficient solvers, do not result in significant savings compared to that incurred by the full finite element discretization. For linear problems the stiffness matrix is independent of the solution so that one can assemble the small reduced basis stiffness matrix during the offline stage. For nonlinear problems, the discrete system changes at each time step (and generally at each interrogation of a nonlinear solver) so that, in general, it is not an easy matter to avoid the high assembly costs. However, because the nonlinearity in the Navier-Stokes system is quadratic, the assembly costs can again be shifted to the offline stage during which one assembles a low-dimensional third-order tensor that can be reused throughout the calculations.

Turning to the computational cost for the FEM ensemble method, as mentioned previously, the most efficient way to solve the resulting systems is a block solver (e.g., block GMRES). In trying to estimate the computational cost, we run into the same problem as we do for estimating the cost of solving the standard FEM discretization of the Navier-Stokes problem; specifically, it is very difficult to precisely determine how quickly any block solver converges.

Due to the difficulties outlined above in a priori estimation of the computational costs for both our algorithms we omit any CPU time comparison in the numerical experiments. Instead, we focus on the accuracy of our En-POD method, demonstrating that it is possible to achieve similar results as those given by the En-full-FE method. A more rigorous and thorough analysis comparing the computational cost of the En-POD and En-full-FE method is a subject of future research.

6.2 Flow between two offset circles

For the numerical experiment we examine the two-dimensional flow between two offset circles with viscosity coefficient . Specifically, the domain is a disk with a smaller offset disc inside. Let , , , and ; then the domain is given by

No-slip, no-penetration boundary conditions are imposed on both circles. All computations are done using the FEniCS software suite [29]. The deterministic flow driven by the counterclockwise rotational body force

displays interesting structures interacting with the inner circle. A Kármán vortex street is formed which then re-interacts with the inner circle and with itself, generating complex flow patterns.

For our test problems, we generate perturbed initial conditions by solving a steady Stokes problem with perturbed body forces given by

with different perturbations defined by varying . We discretize in space via the - Taylor-Hood element pair. Meshes were generated using the FEniCS built-in mshr package with varying refinement levels. An example mesh is given in Figure LABEL:meshEx1.

Figure 6.1: Mesh for flow between offset circles resulting in 16,457 total degrees of freedom for the Taylor-Hood element pair.

In order to generate the POD basis, we use two perturbations of the initial conditions corresponding to and . Using a mesh that results in 16,457 total degrees of freedom and a fixed time step , we run a standard full finite element code111We also generated snapshots using the En-full-FE method. We found that we obtained exactly the same results as those reported here if instead we use a standard finite element method. for each perturbation from to . For the time discretization we use the Crank-Nicolson method and take snapshots every seconds. In Figure LABEL:eigvaldecay, we illustrate the decay of the singular values generated the snapshot matrix.

Figure 6.2: The 40 largest singular values of the snapshot matrix.

6.3 Example 1

The purpose of this example is to illustrate our theoretical error estimates and to show the efficacy of our method in a “data mining” setting, i.e., to show that we can accurately represent the information contained in the En-full-FE approximation which requires the specification of 16,457 coefficients by the En-POD approximation that requires the specification of a much smaller number of coefficients, in fact, merely 10 will do. Thus, we determine the En-POD approximation using the same perturbations, mesh, and time step as were used in the generation of the POD basis. We verify at each time step that condition (LABEL:ineq:CFL-h) is satisfied. In order to illustrate the accuracy of our approach, we provide, in Figure LABEL:noRomEx1, plots of the velocity field of the ensemble average at the final time for both the En-full-FE and En-POD approximations. We also provide in Figure LABEL:errorDiff (left) the difference between the two ensemble averages at the final time . In addition, in Figure LABEL:Energy_ex1, we plot, for and for both methods, the energy and the enstrophy .

Figure 6.3: For Example 1, the ensemble average of the velocity field at the final time of the En-full-FE (left) approximation and the En-POD approximation with reduced basis vectors (right).
Figure 6.4: For Example 1, The difference between the ensemble average of the velocity field at the final time of the En-full-FE (approximation and the En-POD approximation with reduced basis vectors.

We need POD basis functions to reproduce the flow with a reasonable level of accuracy. This is seen in Table LABEL:tabEx1(a) which shows a small discrete error corresponding to 10 basis vectors and, as the number of basis vectors increases beyond that, the error appears to decreases monotonically. Visual confirmation is given by comparing the two plots in Figure LABEL:noRomEx1 as well as Figure LABEL:errorDiff (left); at time the En-POD method appears to produce a flow which is very similar to that for the En-full method. Additionally, in Figure LABEL:Energy_ex1, we plot the energy and enstrophy of En-POD with varying cardinalities for the POD basis and for the En-full-FE method. It can be seen that as the number of POD basis vectors increases our approximation improves with the En-POD energy and enstrophy becoming indistinguishable from that for the En-full-FE for or more POD basis functions.

(a) Example 1