An efficient algorithm for simulating
ensembles of parameterized flow problems
Abstract
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 righthand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and longterm stable under certain conditions on the timestep size and a parameter deviation ratio. Rigorous numerical error estimate shows the scheme is of firstorder accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results.
Key words. NavierStokes 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 lowdimensional surrogates for the PDE solution such as sparsegrid 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 largescale problems because each individual realization requires a largescale 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 NavierStokes 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 righthand sides (RHS), reducing both the storage requirements and computational costs of the solution process. The algorithm of [16] is firstorder accurate in time; it is extended to higherorder 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 reducedorder 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 ensemblebased 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 noslip boundary conditions for which, for , an individual member solves the system
\hb@xt@.01(1.1) 
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 semiimplicit schemes are invariably used for time integration. For a semiimplicit 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 ensemblebased semiimplicit scheme
For clarity, we temporarily suppress the spatial discretization and only consider the ensemblebased implicitexplicit temporal integration scheme
\hb@xt@.01(1.2) 
where and are the ensemble means of the velocity and viscosity coefficient, respectively, defined as
After rearranging the system, we have, at time ,
\hb@xt@.01(1.3) 
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 righthand 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 subensembles so as to balance memory, communication, and computational costs and then (LABEL:FirstOrder) can be applied to each subensemble.
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:FirstOrder). In §LABEL:sec:err, we derive an error estimate for the fullydiscrete 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 infsup (or ) condition required for the stability of the finite element approximation and that the finite element spaces satisfy the approximation properties
\hb@xt@.01(1.4)  
\hb@xt@.01(1.5)  
\hb@xt@.01(1.6) 
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 TaylorHood –, , element pairs. For details concerning finite element methods see [5] and see [7, 8, 9, 20] for finite element methods for the NavierStokes 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
\hb@xt@.01(1.7) 
that is known to hold for standard finite element spaces with locally quasiuniform meshes [3]. We also define the standard explicitly skewsymmetric trilinear form
that satisfies the bounds [20]
\hb@xt@.01(1.8)  
\hb@xt@.01(1.9) 
We also denote the exact and approximate solutions at as and , respectively.
2 Stability analysis
The fullydiscrete finite element discretization of (LABEL:FirstOrder) is given as follows. Given , for , find and satisfying
\hb@xt@.01(2.1) 
We begin by proving the conditional, nonlinear, longtime stability of the scheme (LABEL:FirstOrderh) under a timestep condition and a parameter deviation condition.
Theorem 2.1 (Stability)
For all , if for some , , and some , , the following timestep condition and parameter deviation condition both hold
\hb@xt@.01(2.2)  
\hb@xt@.01(2.3) 
then, the scheme (LABEL:FirstOrderh) 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:CFLh1) that the upper bound in the timestep condition increases as decreases. As , the bound approaches . Because the upper bound for the relative deviation of viscosity coefficient in (LABEL:ineq:CFLh2) is bounded by , the two stability conditions are oppositional to each other.
Remark 2.3
Noting that the condition (LABEL:ineq:CFLh1) only depends on known quantities such as the solution at and that the scheme (LABEL:FirstOrderh) is a onestep method, (LABEL:ineq:CFLh1) 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 timestep 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:FirstOrderh) is equivalent to: Given , for , find such that
\hb@xt@.01(3.1)  
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:FirstOrderh))
For all , if for some , , and some , , the following timestep condition and parameter deviation condition both hold
\hb@xt@.01(3.2)  
\hb@xt@.01(3.3) 
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 TaylorHood elements (, ) are used, i.e., the piecewisequadratic velocity space and the piecewiselinear pressure space , we have the following estimate.
Corollary 3.2
Assume that and are both accurate or better. Then, if is chosen as the TaylorHood element pair, we have
4 Numerical experiments
In this section, we illustrate the effectiveness of our proposed method (LABEL:FirstOrder) and the associated theoretical analyses in §LABEL:sec:stab and §LABEL:sec:err by considering two examples: a GreenTaylor 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 GreenTaylor vortex problem
The GreenTaylor vortex flow is commonly used for testing convergence rates, e.g., see [1, 2, 4, 13, 19, 13, 25]. The GreenTaylor vortex solution given by
\hb@xt@.01(4.1)  
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:CFLh2) 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 
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 
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:CFLh2) 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 
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 
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
The flow is driven by a counterclockwise rotational body force
with noslip 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 reinteracts with that circle and with itself, generating complex flow patterns. We consider multiple numerical simulations of the flow with different viscosity coefficients using the ensemblebased algorithm (LABEL:FirstOrderh). For spatial discretization, we apply the TaylorHood 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:CFLh2) 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 timestep size.
The comparison of the energy evolution of ensemblebased 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:CFLh2) 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 .
5 Conclusions
In this paper, we consider a set of NavierStokes 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 righthand 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 GreenTaylor flow and serves to illustrate the firstorder accuracy in time of the ensemblebased 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 higherorder accurate schemes for the flow ensemble simulations.
References
 [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), 4474.
 [2] L. Berselli, On the large eddy simulation of the TaylorGreen vortex, J. Math. Fluid Mech., 7 (2005), S164S191.
 [3] S. Brenner and R. Scott, The Mathematical Theory of Finite Element Methods, Springer, 3rd edition, 2008.
 [4] A. Chorin, Numerical solution for the NavierStokes equations, Math. Comp., 22 (1968), 745762.
 [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), 203215.
 [7] V. Girault and P.A. Raviart, Finite Element Approximation of the NavierStokes Equations, Lecture Notes in Mathematics, Vol. 749, Springer, Berlin, 1979.
 [8] V. Girault and P.A. Raviart, Finite Element Methods for NavierStokes 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 ensembleproper orthogonal decomposition method for the nonstationary NavierStokes Equations, SIAM J. Numer. Anal., 2016, to appear, arXiv:1603.04777v2 [math.NA].
 [11] M. Gunzburger, N. Jiang and M. Schneier, A higherorder ensemble/proper orthogonal decomposition method for the nonstationary NavierStokes Equations, submitted, 2016.
 [12] J. Holman, Experimental methods for engineers 8th ed., McGrawHill 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), 307330.
 [14] N. Jiang, A higher order ensemble simulation algorithm for fluid flows, J. Sci. Comput., 64 (2015), 264288.
 [15] N. Jiang, A secondorder ensemble method based on a blended backward differentiation formula timestepping scheme for timedependent NavierStokes equations, Numer. Meth. Partial. Diff. Eqs., 33 (2017), 3461.
 [16] N. Jiang and W. Layton, An algorithm for fast calculation of flow ensembles, Int. J. Uncertain. Quantif., 4 (2014), 273301.
 [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), 630651.
 [18] N. Jiang, S. Kaya, and W. Layton, Analysis of model variance for ensemble based turbulence modeling, Comput. Meth. Appl. Math., 15 (2015), 173188.
 [19] V. John and W. Layton, Analysis of numerical errors in large eddy simulation, SIAM J. Numer. Anal., 40 (2002), 9951020.
 [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/cmam20160033.
 [22] D. O’Leary, The block conjugate gradient algorithm and related methods, Lin. Alg. Appl. 29 (1980), 292322.
 [23] C. Powell and D. Silvester, Preconditioning steadystate Navier–Stokes equations with random data, SIAM J. Sci. Comput., 34(5), 2012, A2482A2506.
 [24] B. Sousedík and H. Elman, Stochastic Galerkin methods for the steadystate Navier–Stokes equations, J. Comput. Phys., 316 (2016), 435452.
 [25] D. Tafti, Comparison of some upwindbiased highorder formulations with a second order centraldifference scheme for time integration of the incompressible NavierStokes equations, Comput. & Fluids, 25 (1996), 647665.
 [26] A. Takhirov, M. Neda and Jiajia Waters, Time relaxation algorithm for flow ensembles, Numer. Meth. Partial. Diff. Eqs., 32 (2016), 757777.
A Proof of Theorem LABEL:th:FirstOrder
Setting and in (LABEL:FirstOrderh) and then adding two equations, we obtain
Applying Young’s inequality to the terms on the RHS yields, for ,
\hb@xt@.01(A.1)  
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
\hb@xt@.01(A.2)  