An efficient algorithm for simulatingensembles of parameterized flow problems

An efficient algorithm for simulating
ensembles of parameterized flow problems

Max Gunzburger Department of Scientific Computing, Florida State University, Tallahassee, FL 32306-4120 Research supported by the U.S. Air Force Office of Scientific Research grant FA9550-15-1-0001 and the U.S. Department of Energy Office of Science grants DE-SC0009324 and DE-SC0016591.    Nan Jiang Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409-0020    Zhu Wang Department of Mathematics, University of South Carolina, Columbia, SC 29208 Research supported by the U.S. National Science Foundation grant DMS-1522672 and the U.S. Department of Energy grant DE-SC0016540.

Many applications of computational fluid dynamics require multiple simulations of a flow under different input conditions. In this paper, a numerical algorithm is developed to efficiently determine a set of such simulations in which the individually independent members of the set are subject to different viscosity coefficients, initial conditions, and/or body forces. The proposed scheme applied to the flow ensemble leads to need to solve a single linear system with multiple right-hand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and long-term stable under certain conditions on the time-step size and a parameter deviation ratio. Rigorous numerical error estimate shows the scheme is of first-order accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results.

Key words. Navier-Stokes equations, ensemble simulations, ensemble method

1 Introduction

Numerical simulations of incompressible viscous flows have important applications in engineering and science. In this paper, we consider settings in which one wishes to obtain solutions for several different values of the physical parameters and several different choices for the forcing functions appearing in the partial differential equation (PDE) model. For example, in building low-dimensional surrogates for the PDE solution such as sparse-grid interpolants or proper orthogonal decomposition approximations, one has to first determine expensive approximation of solutions corresponding to several values of the parameters. Sensitivity analysis of solutions is setting in which one often has to determine approximate solutions for several parameter values and/or forcing functions. An important third example is quantifying the uncertainties of outputs from the model equations. Mathematical models should take into account the uncertainties invariably present in the specification of physical parameters and/or forcing functions appearing in the model equations. For flow problems, because the viscosity of the liquid or gas often depends on the temperature, an inaccurate measurement of the temperature would introduce some uncertainty into the viscosity of the flow. Direct measurements of the viscosity using flow meters and measurements of the state of the system are also prone to uncertainties. Of course, forcing functions, e.g., initial condition data, can and usually are also subject to uncertainty. In such cases, due to the lack of of exact information, stochastic modeling is used to describe flows subject to a random viscosity coefficient and/or random forcing. Subsequently, numerical methods are employed to quantify the uncertainties in system output. It is known that uncertainty quantification (UQ), when a random sampling method such as Monte Carlo method is used, could be computationally expensive for large-scale problems because each individual realization requires a large-scale computation but on the other hand, many realizations may be needed in order to obtain accurate statistical information about the outputs of interest. Therefore, for all the examples discussed and for many others, how to design efficient algorithms for performing multiple numerical simulations becomes a matter of great interest.

The ensemble method which forms the basis for our approach was proposed in [16]; there, a set of solutions of the Navier-Stokes equations (NSE) with distinct initial conditions and forcing terms is considered. All solutions are found, at each time step, by solving a linear system with one shared coefficient matrix and right-hand sides (RHS), reducing both the storage requirements and computational costs of the solution process. The algorithm of [16] is first-order accurate in time; it is extended to higher-order accurate schemes in [14, 15]. Ensemble regularization methods are developed in [14, 17, 26] for high Reynolds number flows, and a turbulence model based on ensemble averaging is developed in [18]. The ensemble algorithm has also been extended to simulate MHD flows in [21]. Ensemble algorithms incorporating reduced-order modeling techniques are studied in [10, 11]. It is worth mentioning that all the ensemble algorithms developed so far can only deal with simulations subject to different initial conditions and/or body forces, but not other model parameters.

In this paper, we develop a numerical scheme for ensemble-based simulations of the NSE in which not only the initial data and body force function, but also the viscosity coefficient, may vary from one ensemble member to another. Specifically, we consider a set of NSE simulations on a bounded domain subject to no-slip boundary conditions for which, for , an individual member solves the system


which corresponds, for each , to a different viscosity coefficient and/or distinct initial data and/or body forces .

Due to the nonlinear convection term, implicit and semi-implicit schemes are invariably used for time integration. For a semi-implicit scheme, the associated discrete linear systems would be different for each individually independent simulation, i.e., for each . As a result, at each time step, linear systems need to be solved to determine the ensemble, resulting in a huge computational effort. For a fully implicit scheme, the situation is even worse because one would have to solve many more linear systems due to the nonlinear solver iteration. To tackle this issue, we propose a novel discretization scheme that results, at each time step, in a common coefficient matrix for all the ensemble members.

1.1 The ensemble-based semi-implicit scheme

For clarity, we temporarily suppress the spatial discretization and only consider the ensemble-based implicit-explicit temporal integration scheme


where and are the ensemble means of the velocity and viscosity coefficient, respectively, defined as

After rearranging the system, we have, at time ,


It is clear that the coefficient matrix of the resulting linear system will be independent of . Thus, for the flow ensemble, to advance all members of the ensemble one time step, we need only solve a single linear system with right-hand sides. Compared with solving individually independent simulations, this approach used with a block solver such as a block generalized CG method [6, 22] is much more efficient and significantly reduces the required storage. When the size of the ensemble becomes huge, it can be subdivided into sub-ensembles so as to balance memory, communication, and computational costs and then (LABEL:First-Order) can be applied to each sub-ensemble.

The rest of this section is devoted to establishing notation and to providing other preliminary information. Then, in §LABEL:sec:stab, we prove a conditional stability result for a fully discrete finite element discretization of (LABEL:First-Order). In §LABEL:sec:err, we derive an error estimate for the fully-discrete approximation. Results of the preliminary numerical simulations that illustrate the theoretical results are given in §LABEL:sec:num and §LABEL:sec:con provides some concluding remarks.

1.2 Notation and preliminaries

Let denote an open, regular domain in for or having boundary denoted by . The norm and inner product are denoted by and , respectively. The norms and the Sobolev norms are denoted by and , respectively. The Sobolev space is simply denoted by and its norm by . For functions defined on , we define, for ,

Given a time step , associated discrete norms are defined as

where and . Denote by the dual space of bounded linear functions on ; a norm on is given by

The velocity space and pressure space are given by

respectively. The space of weakly divergence free functions is

A weak formulation of (LABEL:eq:NSE) reads: for , find and for a.e. satisfying

with .

Our analysis is based on a finite element method (FEM) for spatial discretization. However, the results also extend, without much difficulty, to other variational discretization methods. Let and denote families of conforming velocity and pressure finite element spaces on regular subdivision of into simplicies; the family is parameterized by the maximum diameter of any of the simplicies. Assume that the pair of spaces satisfy the discrete inf-sup (or ) condition required for the stability of the finite element approximation and that the finite element spaces satisfy the approximation properties


where the generic constant is independent of mesh size . An example for which the stability condition and the approximation properties are satisfied is the family of Taylor-Hood , , element pairs. For details concerning finite element methods see [5] and see [7, 8, 9, 20] for finite element methods for the Navier-Stokes equations.

The discretely divergence free subspace of is defined as

Note that, in general, . We assume the mesh and finite element spaces satisfy the standard inverse inequality


that is known to hold for standard finite element spaces with locally quasi-uniform meshes [3]. We also define the standard explicitly skew-symmetric trilinear form

that satisfies the bounds [20]


We also denote the exact and approximate solutions at as and , respectively.

2 Stability analysis

The fully-discrete finite element discretization of (LABEL:First-Order) is given as follows. Given , for , find and satisfying


We begin by proving the conditional, nonlinear, long-time stability of the scheme (LABEL:First-Order-h) under a time-step condition and a parameter deviation condition.

Theorem 2.1 (Stability)

For all , if for some , , and some , , the following time-step condition and parameter deviation condition both hold


then, the scheme (LABEL:First-Order-h) is nonlinearly, long time stable. In particular, for and for any , we have

Proof. The proof is given in Appendix LABEL:proofa1.     

Remark 2.2

It is seen from (LABEL:ineq:CFL-h1) that the upper bound in the time-step condition increases as decreases. As , the bound approaches . Because the upper bound for the relative deviation of viscosity coefficient in (LABEL:ineq:CFL-h2) is bounded by , the two stability conditions are oppositional to each other.

Remark 2.3

Noting that the condition (LABEL:ineq:CFL-h1) only depends on known quantities such as the solution at and that the scheme (LABEL:First-Order-h) is a one-step method, (LABEL:ineq:CFL-h1) can be used to adapt in order to guarantee the stability for the ensemble simulations.

3 Error Analysis

In this section, we give a detailed error analysis of the proposed method under the same type of time-step condition (with possibly different constant on the left hand side of the inequality) and the same parameter deviation condition. Assuming that and satisfy the condition, the scheme (LABEL:First-Order-h) is equivalent to: Given , for , find such that


To analyze the rate of convergence of the approximation, we assume that the following regularity for the exact solutions:

Let denote the approximation error of the -th simulation at the time instance . We then have the following error estimates.

Theorem 3.1 (Convergence of scheme (LABEL:First-Order-h))

For all , if for some , , and some , , the following time-step condition and parameter deviation condition both hold


then, there exists a positive constant independent of the time step such that

Proof. The proof is given in Appendix LABEL:proofa2.     

In particular, when Taylor-Hood elements (, ) are used, i.e., the piecewise-quadratic velocity space and the piecewise-linear pressure space , we have the following estimate.

Corollary 3.2

Assume that and are both accurate or better. Then, if is chosen as the Taylor-Hood element pair, we have

4 Numerical experiments

In this section, we illustrate the effectiveness of our proposed method (LABEL:First-Order) and the associated theoretical analyses in §LABEL:sec:stab and §LABEL:sec:err by considering two examples: a Green-Taylor vortex problem and a flow between two offset cylinders. The first problem has a known exact solution that is used to illustrate the error analysis. The second example does not have an analytic solution but has complex flow structures; it is used to check the stability analysis.

4.1 The Green-Taylor vortex problem

The Green-Taylor vortex flow is commonly used for testing convergence rates, e.g., see [1, 2, 4, 13, 19, 13, 25]. The Green-Taylor vortex solution given by


satisfies the NSE in for and initial condition

The solution consists of an array of oppositely signed vortices that decay as . In the following numerical tests, we take , , , , and . The boundary condition is assumed to be inhomogeneous Dirichlet, that is, the boundary values match that of the exact solution.

We consider an ensemble of two members, and , corresponding to two incompressible NSE simulations with different viscosity coefficients and initial conditions . We investigate the ensemble simulations and compare it with independent simulations. For , we define by the approximation error of the -th member of the ensemble simulation and by the approximation error of the -th independently determined simulation. Here, the superscript “” stands for “ensemble” whereas “” stands for “independent.”

Case 1

We set the viscosity coefficient and initial condition for the first member and and for the second member , where . For this choice of parameters, we have for both and so that the condition (LABEL:ineq:CFL-h2) is satisfied. We first apply the ensemble algorithm; results are shown in Table LABEL:tab:t6ensemble. It is seen that the convergence rate for and is first order.

rate rate rate rate
0.85 0.91 0.93 0.95
0.92 0.95 0.94 0.97
0.96 0.97 0.97 0.98
Table 4.1: For the Green-Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes , errors for ensemble simulations of two members with inputs , and , .

We next compare the ensemble simulations with independent simulations. To this end, we perform each NSE simulation independently using the same discretization setup. The associated approximation errors are listed in Table LABEL:tab:t6u1. Comparing with Table LABEL:tab:t6ensemble, we observe that the ensemble simulation is able to achieve accuracies close to that of the independent, more costly simulations.

rate rate rate rate
0.89 0.93 0.90 0.93
0.94 0.96 0.93 0.96
0.97 0.98 0.97 0.98
Table 4.2: For the Green-Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes , errors in independent simulations of two members with inputs , and , .

Case 2

We now set and while keeping the same initial conditions as for Case 1. With this choice of parameters, for both and , which still satisfies (LABEL:ineq:CFL-h2) but is closer to the upper limit. The ensemble simulation errors are listed in Table LABEL:tab:t2ensemble, which shows the rate of convergence for the second member is nearly 1 and for the first member is approaching 1.

rate rate rate rate
0.65 0.71 1.08 1
0.78 0.83 0.95 0.98
0.87 0.90 0.98 0.98
Table 4.3: For the Green-Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes , errors in ensemble simulations of two members: , and , .

The approximation errors for two independent simulations under using the same discretization setup are listed in Table LABEL:tab:t2u1. Comparing the ensemble simulation results in Table LABEL:tab:t2ensemble with the independent simulations, we find that the accuracy of first member in the ensemble simulation degrades slightly whereas that of the second member in the ensemble simulation improves a bit. Overall, the ensemble simulation is able to achieve the same order of accuracy as the independent simulations.

rate rate rate rate
0.93 0.93 0.86 0.94
0.97 0.97 0.93 0.96
0.98 0.98 0.96 0.98
Table 4.4: For the Green-Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes , errors in independent simulations of two members: , and , .

4.2 Flow between two offset cylinders

Next, we check the stability of our algorithm by considering the problem of a flow between two offset circles [14, 16, 17, 18]. The domain is a disk with a smaller off center obstacle inside. Letting , , and , the domain is given by

Figure 4.1: Mesh for the flow between two offset cylinders example.

The flow is driven by a counterclockwise rotational body force

with no-slip boundary conditions imposed on both circles. The flow between the two circles shows interesting structures interacting with the inner circle. A Von Krmn vortex street is formed 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 (LABEL:First-Order-h). For spatial discretization, we apply the Taylor-Hood element pair on a triangular mesh that is generated by Delaunay triangulation with mesh points on the outer circle and mesh points on the inner circle and with refinement near the inner circle, resulting in degrees of freedom; see Figure LABEL:mesh.

In order to illustrate the stability analysis, we select two different sets of viscosity coefficients for:

  • Case 1:  ,   ,  

  • Case 2:  ,   ,   .

The average of the viscosity coefficients is for both cases. However, the stability condition (LABEL:ineq:CFL-h2) is satisfied in the first case but is not satisfied in the second one, i.e., we have

  • Case 1:  ,   ,  

  • Case 2:  ,   ,   .

The second member of Case 2 has a perturbation ratio greater than 1. Simulations of both cases are subject to the same initial condition and body forces for all ensemble members. In particular, the initial condition is generated by solving the steady Stokes problem with viscosity and the same body force . All the simulations are run over the time interval with a time step size . For the stability test, we use the kinetic energy as a criterion and compare the ensemble simulation results with independent simulations using the same mesh and time-step size.

The comparison of the energy evolution of ensemble-based simulations with the corresponding independent simulations is shown in Figures LABEL:egy_1 and LABEL:egy_2. It is seen that, for Case 1, the ensemble simulation is stable, but for Case 2, it becomes unstable. This phenomena coincides with our stability analysis since the condition (LABEL:ineq:CFL-h2) holds for all members of Case 1, but does not hold for the second member of Case 2. Indeed, it is observed from Figure LABEL:egy_2 that the energy of the second member in Case 2 blows up after , then affecting other two members and results in their energy dramatically increase after .

Figure 4.2: For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.).
Figure 4.3: For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.).

5 Conclusions

In this paper, we consider a set of Navier-Stokes simulations in which each member may be subject to a distinct viscosity coefficient, initial conditions, and/or body forces. An ensemble algorithm is developed for the group by which all the flow ensemble members, after discretization, share a common linear system with different right-hand side vectors. This leads to great saving in both storage requirements and computational costs. The stability and accuracy of the ensemble method are analyzed. Two numerical experiments are presented. The first is for Green-Taylor flow and serves to illustrate the first-order accuracy in time of the ensemble-based scheme. The second is for a flow between two offset cylinders and serves to show that our stability analysis is sharp. As a next step, we will investigate higher-order accurate schemes for the flow ensemble simulations.


  • [1] D. Barbato, L. Berselli, and C. Grisanti, Analytical and numerical results for the rational large eddy simulation model, J. Math. Fluid Mech., 9 (2007), 44-74.
  • [2] L. Berselli, On the large eddy simulation of the Taylor-Green vortex, J. Math. Fluid Mech., 7 (2005), S164-S191.
  • [3] S. Brenner and R. Scott, The Mathematical Theory of Finite Element Methods, Springer, 3rd edition, 2008.
  • [4] A. Chorin, Numerical solution for the Navier-Stokes equations, Math. Comp., 22 (1968), 745-762.
  • [5] P. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, Philadelphia, 2002.
  • [6] Y. Feng, D. Owen, and D. Peric, A block conjugate gradient method applied to linear systems with multiple right hand sides, Comp. Meth. Appl. Mech. Eng., 127 (1995), 203-215.
  • [7] V. Girault and P.-A. Raviart, Finite Element Approximation of the Navier-Stokes Equations, Lecture Notes in Mathematics, Vol. 749, Springer, Berlin, 1979.
  • [8] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations - Theory and Algorithms, Springer, Berlin, 1986.
  • [9] M. Gunzburger, Finite Element Methods for Viscous Incompressible Flows - A Guide to Theory, Practice, and Algorithms, Academic Press, London, 1989.
  • [10] M. Gunzburger, N. Jiang and M. Schneier, An ensemble-proper orthogonal decomposition method for the nonstationary Navier-Stokes Equations, SIAM J. Numer. Anal., 2016, to appear, arXiv:1603.04777v2 [math.NA].
  • [11] M. Gunzburger, N. Jiang and M. Schneier, A higher-order ensemble/proper orthogonal decomposition method for the nonstationary Navier-Stokes Equations, submitted, 2016.
  • [12] J. Holman, Experimental methods for engineers 8th ed., McGraw-Hill series in mechanical engineering, 2011
  • [13] N. Jiang and H. Tran, Analysis of a stabilized CNLF method with fast slow wave splittings for flow problems, Comput. Meth. Appl. Math., 15 (2015), 307-330.
  • [14] N. Jiang, A higher order ensemble simulation algorithm for fluid flows, J. Sci. Comput., 64 (2015), 264-288.
  • [15] N. Jiang, A second-order ensemble method based on a blended backward differentiation formula timestepping scheme for time-dependent Navier-Stokes equations, Numer. Meth. Partial. Diff. Eqs., 33 (2017), 34-61.
  • [16] N. Jiang and W. Layton, An algorithm for fast calculation of flow ensembles, Int. J. Uncertain. Quantif., 4 (2014), 273-301.
  • [17] N. Jiang and W. Layton, Numerical analysis of two ensemble eddy viscosity numerical regularizations of fluid motion, Numer. Meth. Partial. Diff. Eqs., 31 (2015), 630-651.
  • [18] N. Jiang, S. Kaya, and W. Layton, Analysis of model variance for ensemble based turbulence modeling, Comput. Meth. Appl. Math., 15 (2015), 173-188.
  • [19] V. John and W. Layton, Analysis of numerical errors in large eddy simulation, SIAM J. Numer. Anal., 40 (2002), 995-1020.
  • [20] W. Layton, Introduction to the Numerical Analysis of Incompressible Viscous Flows, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2008.
  • [21] M. Mohebujjaman and L. Rebholz, An efficient algorithm for computation of MHD flow ensembles, Comput. Meth. Appl. Math., 2016, in press, DOI: 10.1515/cmam-2016-0033.
  • [22] D. O’Leary, The block conjugate gradient algorithm and related methods, Lin. Alg. Appl. 29 (1980), 292-322.
  • [23] C. Powell and D. Silvester, Preconditioning steady-state Navier–Stokes equations with random data, SIAM J. Sci. Comput., 34(5), 2012, A2482-A2506.
  • [24] B. Sousedík and H. Elman, Stochastic Galerkin methods for the steady-state Navier–Stokes equations, J. Comput. Phys., 316 (2016), 435-452.
  • [25] D. Tafti, Comparison of some upwind-biased high-order formulations with a second order central-difference scheme for time integration of the incompressible Navier-Stokes equations, Comput. & Fluids, 25 (1996), 647-665.
  • [26] A. Takhirov, M. Neda and Jiajia Waters, Time relaxation algorithm for flow ensembles, Numer. Meth. Partial. Diff. Eqs., 32 (2016), 757-777.

A Proof of Theorem LABEL:th:First-Order

Setting and in (LABEL:First-Order-h) and then adding two equations, we obtain

Applying Young’s inequality to the terms on the RHS yields, for ,


Both and on the RHS of (LABEL:ineq:_tri) need to be absorbed into on the LHS. To this end, we minimize by selecting so that (LABEL:ineq:_tri) becomes