A Second-Order Time-Stepping Scheme for Simulating Ensembles of Parameterized Flow Problems

1 Introduction

Many computational fluid dynamics applications require multiple simulations of a flow under different input conditions. For example, the ensemble Kalman filter approach used in data assimilation first simulates a forward model a large number of times by perturbing either the initial condition data, boundary condition data, or uncertain parameters, then corrects the model based on the model forecasts and observational data. A second example is the construction of low-dimensional surrogates for partial differential equation (PDE) solutions such as sparse-grid interpolants or proper orthogonal decomposition approximations, for which one has to first obtain expensive approximations of solutions corresponding to several parameter samples. Another example is sensitivity analyses of solutions for which one often has to determine approximate solutions for a number of perturbed inputs such as the values of certain physical parameters. In this paper, we consider such applications and develop a second-order time-stepping scheme for efficiently simulating an ensemble of flows. In particular, we consider the setting in which one wishes to determine the PDE solutions for several different choices of initial condition and boundary condition data, forcing functions, and physical parameters appearing in the PDE model.

The ensemble algorithm was first developed in [18] to find a set of J solutions of the Navier–Stokes equations (NSE) subject to different initial condition and forcing functions. The main idea is that, based on the introduction of an ensemble average and a special semi-implicit time discretization, the discrete systems for the multiple flow simulations share a common coefficient matrix. Thus, instead of solving J linear systems with J different matrices and right-hand sides, one only need solve a single linear system with J right-hand side vectors. This leads to very significant computational savings in linear solver costs when LU factorization (for small-scale systems) or a block iterative algorithm (for large-scale systems) are used. High-order ensemble algorithms were designed in [15, 16]. For high Reynolds number flows, ensemble regularization methods and a turbulence model based on ensemble averaging have been developed in [15, 19, 17, 24]. The method has also been extended to simulate MHD flows in [23] and to develop ensemble-based reduced-order modeling techniques in [9, 8]. In [10], the authors proposed a first-order ensemble algorithm that deals with a number of flow simulations subject to not only different initial condition, boundary conditions, and/or body force data, but also distinct viscosity coefficients appearing in the NSE model. In this paper, we follow the same direction and develop an ensemble scheme having higher accuracy.

To begin, consider an ensemble of incompressible flow simulations on a bounded domain subject to Dirichlet boundary conditions. The j-th member of the ensemble is a simulation associated with the positive viscosity coefficient νj, initial condition data uj0, boundary condition data gj, and body force fj. All of these data may vary from one simulation to another. Then, for j=1,…,J, we need to solve

(1.1)

There is a long list of work in developing time discretization methods for the NSE including explicit, implicit, and semi-implicit schemes, for example, [11, 13, 12, 14, 20, 22, 25]. In general, explicit schemes are easier to implement, but they suffer the severely restricted time step size from stability requirement. The fully implicit and semi-implicit schemes have better stability conditions, but the discretization would lead to a varying coefficient matrix of the system. As a result, a different linear system has to be solved for each member at every time step, thus totally J linear systems need to be solved per time step. To overcome this issue, we propose a new, second-order accurate in time, ensemble-based scheme that improves the computational efficiency. The scheme is semi-implicit that permits the use of a known quantity (the ensemble mean defined below), which is independent of the ensemble index j, in the advection term and, therefore, leads to a single coefficient matrix for all the ensemble members.

For keeping the exposition simple, we consider a uniform time step Δ⁢t and let tn=n⁢Δ⁢t for n=0,1,…. We then consider the ensemble of semi-discrete in time systems

where ujn, pjn and fjn denote approximations of uj⁢(⋅,tn), pj⁢(⋅,tn) and fj⁢(⋅,tn) of (1.1), respectively. In (1.2), u¯n and ν¯ denote the ensemble means of the velocity field and viscosity coefficient, respectively, defined by

and uj′⁣n represents the fluctuation defined by

It is easy to see that the coefficient matrix in the spatial discretization of (1.2) does not depend on j. Thus, all the members in the ensemble do share a common coefficient matrix. To advance one time step, one only need solve a single linear system with J right-hand side vectors, which is more efficient than solving J individual simulations.

We assume homogeneous flow boundary conditions (gj=0) in the following derivation and analysis of the proposed ensemble algorithm. But flows with inhomogeneous essential boundary conditions are considered in our first numerical experiment presented in Section 5, where, in the implementation, the data gj at each time step is first replaced by its interpolant on the Lagrangian finite element space and then is enforced on the boundary nodes. The extension of our analysis to the inhomogeneous cases will follow the idea presented in [3], which will be discussed elsewhere.

In what follows, we present a rigorous theoretical analysis of the proposed scheme. In Section 2, we provide some notations and preliminaries; in Section 3, the stability conditions of the scheme are obtained; and in Section 4, an error estimate is derived. Then, several numerical experiments are presented in Section 5.

2 Notation and Preliminaries

Let Ω be an open, regular domain in ℝd(d=2⁢ or ⁢3). The space L2⁢(Ω) is equipped with the norm ∥⋅∥ and inner product (⋅,⋅). Denote by ∥⋅∥Lp and ∥⋅∥Wpk, respectively, the norms for Lp⁢(Ω) and the Sobolev space Wpk⁢(Ω). Let Hk⁢(Ω) be the Sobolev space W2k⁢(Ω) equipped with the norm ∥⋅∥k. For functions v⁢(x,t) defined on (0,T), we define (1≤m<∞)

Given a time step Δ⁢t, let vn=v⁢(tn) and define the discrete norms

Denote by H-1⁢(Ω) the dual space of bounded linear functionals defined on

A norm for H-1⁢(Ω) is given by

We choose the velocity space X and pressure space Q to be

The space of weakly divergence free functions is then

A weak formulation of (1.1) reads: find uj:[0,T]→X and pj:[0,T]→Q for a.e. t∈(0,T] satisfying, for j=1,…,J,

with uj⁢(x,0)=uj0⁢(x).

For the spatial discretization, we use a finite element (FE) method. However, the results can be extended to many other variational methods without much difficulty. Denote by Xh⊂X and Qh⊂Q the conforming velocity and pressure FE spaces on an edge to edge triangulation of Ω with h denoting the maximum diameter of the triangles. Assume that the pair of spaces (Xh,Qh) satisfy the discrete inf-sup (or LBBh) condition, that is required to guarantee the stability of FE approximations. We also assume that the FE spaces satisfy the following approximation properties [21]:

where the generic constant C>0 is independent of the mesh size h. One example for which the LBBh stability condition is satisfied is the family of Taylor–Hood Ps+1-Ps element pairs (i.e., k=s+1 in the definition of Xh), for s≥1 [7]. The discrete divergence free subspace of Xh is

We assume the mesh and FE spaces satisfy the following standard inverse inequality (typical for locally quasi-uniform meshes and standard FE spaces, see, e.g., [2]): for all vh∈Xh,

Define the explicitly skew-symmetric trilinear form

that satisfies the bounds [21]

where C⁢(Ω) is a constant depending on the domain. Denote the exact solution to (1.1) and the FE approximate solution to (2.7) at t=tn by ujn and uj,hn, respectively.

The fully discrete finite element discretization of (1.2) at tn+1 is as follows: given uj,hn, find uj,hn+1∈Xh and pj,hn+1∈Qh satisfying

This is a two-step method that requires uj,h0 and uj,h1 to start the time stepping; uj,h0 is determined by the initial condition and uj,h1 can be computed by the first-order ensemble algorithm developed by the authors in [10] (which is locally second-order accurate) or by using a standard, non-ensemble time-stepping method to compute each individual simulation at the very first time step. Compared to the second-order ensemble scheme developed in [15] for the NSEs without variations in the viscosity coefficient, the scheme (2.7) for parametrized flows introduces an additional average of the viscosity coefficients for the parameterized flow ensemble. As is shown in the next section, the deviation of the flow viscosity coefficients from the ensemble average will play an important role in the stability analysis of the scheme.

3 Stability Analysis

We begin by proving the conditional, nonlinear, long time stability of (2.7) under conditions on the time step and viscosity coefficient deviation: for any j=1,…,J, there exist 0≤μ<1 and 0<ϵ≤2-2⁢μ such that

where C denotes a generic constant depending on the domain and the minimum angle of the mesh.

4 Error Analysis

In this section we derive the numerical error estimate of the proposed ensemble scheme (2.7). We first give a lemma on the estimate of the consistency error of the backward differentiation formula that will be used in the error analysis for the fully discrete ensemble scheme.

By assuming that Xh and Qh satisfy the LBBh condition, the ensemble scheme (2.7) is equivalent to: forn=1,…,N-1, finduj,hn+1∈Vhsuch that

(4.2)

To analyze the rate of convergence of the approximation, we assume the regularity assumptions on the NSE given by

Let ejn=ujn-uj,hn be the error between the true solution of (1.1) and the approximate solution determined from (4.2). We then have the following error estimates.

It is well known that the Taylor–Hood Ps+1-Ps element pairs for which the LBBh stability condition and the approximation properties (2.1)–(2.3) are all satisfied [2, 7, 21]. In particular, when the popular P2-P1 Taylor–Hood FE is used (i.e., k=2 and s=1 in the definitions of Xh and Qh, respectively), we have the following optimal convergence results.

5 Numerical Experiments

The goal of the numerical experiments is two-fold:

1. to numerically illustrate the convergence rate of the ensemble algorithm (2.7), that is, illustrate the second-order accuracy in time.

2. to explore the stability of the algorithm.

In particular, the numerical results strongly indicate that the stability condition (3.2) is sharp.

5.2 Stability Tests

Two conditions, (3.1) and (3.2), guarantee the stability of the proposed scheme. Condition (3.2), in many applications, relates to the probability distribution of the uncertain physical parameters. This requirement on the parameter deviation ratio can be easily checked. If it is not fulfilled, one could divide the parameter sample set into smaller subsets so that it holds on each subset. Condition (3.1) depends on the nature of nonlinear problems. Its severity depends on the governing equations, domain, model parameters, initial/boundary conditions, forcing terms, etc. In practice, once condition (3.2) holds, condition (3.1) can be satisfied by making Δ⁢t sufficiently small and/or by dividing the ensemble into smaller ensembles. Of course, when the ensemble consists of high Reynolds number flows, this condition could easily fail due to the requirement of having an extremely small time-step size leading to a prohibitive computational cost. Condition (3.1) has been investigated extensively in [18, 19, 15]. Hence, in the following, we are mostly interested on the optimality of condition (3.2). However, we do consider the conditional stability due to (3.1) because we want to determine values of the parameter for which that condition is satisfied; this is not directly computable from (3.1) because of the generic constant appearing in that condition. Note that (3.2) contains no such constant so that we can directly study the sharpness of the condition.

We check the stability of our algorithm by using the problem of a flow between two offset circles (see, e.g., [18, 15, 19, 17]). The domain is a disk with a smaller off-center obstacle inside. By letting r1=1, r2=0.1, and c=(c1,c2)=(12,0), the domain is given by

Ω={(x,y):x2+y2≤r12⁢ and ⁢(x-c1)2+(y-c2)2≥r22}.

The flow is driven by a counterclockwise rotational body force

(5.1)f⁢(x,y,t)=[-6⁢y⁢(1-x2-y2),6⁢x⁢(1-x2-y2)]⊤

with no-slip boundary conditions imposed on both circles. A von Kármán vortex street forms behind the inner circle and then re-interacts with that circle and with itself, generating complex flow patterns. We consider multiple numerical simulations of the flow with different viscosity coefficients using the ensemble-based algorithm (2.7). For spatial discretization, we apply the P2-P1 Taylor–Hood element pair on a triangular mesh that is generated by a Delaunay triangulation with 80 mesh points on the outer circle and 60 mesh points on the inner circle and with refinement near the inner circle, resulting in 18,638 degrees of freedom; see Figure 1.

Figure 1

Mesh for the flow between two offset cylinders.

In order to illustrate the stability analysis based on conditions (3.1) and (3.2), we design two numerical tests involving two different sets of viscosity coefficients within an ensemble of three members, keeping the rest of computational setting, including the initial and boundary conditions and body force, the same for all the members. In particular, the initial condition is generated by solving the steady Stokes problem with viscosity ν=0.03 and the same body force f⁢(x,y,t) given by (5.1). We have two test cases:

1. Case 1: ν1=0.021, ν2=0.030, ν3=0.039.

2. Case 2: ν1=0.019, ν2=0.030, ν3=0.041.

Note that the viscosity coefficients ν1 and ν3 for Case 2 are obtained by making small perturbations from those for Case 1 with the average of the viscosity coefficients ν¯=0.03 being the same for in both cases. However, the stability condition (3.2) holds in the first case but breaks down in the second case. In fact, the parameter deviation ratios are given by

1. Case 1: |ν1-ν¯|ν¯=310<13, |ν2-ν¯|ν¯=0<13, |ν3-ν¯|ν¯=310<13.

2. Case 2: |ν1-ν¯|ν¯=1130>13, |ν2-ν¯|ν¯=0<13, |ν3-ν¯|ν¯=1130>13.

For the stability test, we use the kinetic energy as a criterion and compare the ensemble simulation results to the independent simulations using the same mesh and time-step size.

Case 1.

Condition (3.2) is satisfied so that this case illustrates the conditional stability due to (3.1). As mentioned above, we also use this test to determine a value for Δ⁢t for which (3.1) is satisfied so that, in Case 2, we can study the sharpness of condition (3.2). We first test the ensemble-based algorithm at a large time step size Δ⁢t=0.5. The corresponding evolutions of the energy of all the three members are plotted in Figure 2. It is seen that for Δ⁢t=0.5, the algorithm is unstable because the energy of the third member increases dramatically after t=4 and that of the first member after t=4.5. Although not shown in this figure, the energy of the second member also blows up but not until after t=20. This implies that the stability condition (3.1) does not hold. Therefore, we next decrease the time step size to Δ⁢t=0.05 and re-run the ensemble simulations. The associated evolutions of the energies are shown in Figure 2, indicating that the algorithm is now stable over the same time interval. Indeed, additional numerical experiments show that, for any time step size smaller than 0.05, the algorithm is always stable in Case 1 over a much longer time interval, for instance, [0,50]. When an even smaller time step size Δ⁢t=0.01 is selected, the comparison of the energy evolutions of ensemble-based simulations with the corresponding independent simulations over the time interval [0,5] is given in Figures 3. The ensemble simulation is obviously stable and the output energy approximations are very close to those of the independent simulations.

Figure 2

For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble simulations for Δ⁢t=0.5and Δ⁢t=0.05.

Figure 3

For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independentsimulations (Ind.) for Δ⁢t=0.01.

Case 2.

We run ensemble simulations using the small time step size Δ⁢t=0.01 over the same time interval as that for Case 1. As we mentioned above, the viscosity coefficients in Case 2 are obtained by slightly perturbing those in Case 1; this is the only difference between the two computational settings. Since Δ⁢t is chosen small, we believe condition (3.1) still holds for Case 2. But condition (3.2) no longer holds. Therefore, we expect the ensemble simulation to be unstable even when using the small time step size Δ⁢t=0.01. The plots of energy evolutions in Figure 4 matches our expectation as they clearly indicate that the ensemble simulation is unstable in this case. In fact, the energy of the third member blows up after t=1.95 and then affects the other two members and results in their energy dramatically increasing after t=2.45.

Figure 4

For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.) for Δ⁢t=0.01.

6 Conclusions

In this paper, we develop a second-order time-stepping ensemble scheme to compute a set of Navier–Stokes equations in which every member is subject to an independent computational setting including a distinct viscosity coefficient, initial condition data, boundary condition data, and/or body force. By using the ensemble algorithm, all ensemble members share a common coefficient matrix after discretization, although with different right-hand side vectors. Therefore, many efficient block iterative solvers such as the block CG and block GMRES can be applied to solve such a single linear system with multiple right-hand side vectors, leading to great savings in both storage and simulation time. A rigorous analysis shows the proposed algorithm is conditionally, nonlinearly and long-time stable provided two explicit conditions hold and is second-order accurate in time. Two numerical experiments are presented that illustrate our theoretical analysis. In particular, the first is a test problem having an analytic solution that serves to illustrate that the rate of convergence with respect to the time-step size is indeed second order whereas the second example is for a flow between two offset cylinders and shows that the stability condition is sharp. For future work, we plan to investigate the performance of the ensemble algorithm in data assimilation applications.