Ensemble Timestepping Algorithms for Natural Convection

# Ensemble Timestepping Algorithms for Natural Convection

J. A. Fiordilino University of Pittsburgh, Department of Mathematics, Pittsburgh, PA 15260    S. Khankan11footnotemark: 1
###### Abstract

This paper presents two algorithms for calculating an ensemble of solutions to laminar natural convection problems. The ensemble average is the most likely temperature distribution and its variance gives an estimate of prediction reliability. Solutions are calculated by solving two coupled linear systems, each involving a shared coefficient matrix, for multiple right-hand sides at each timestep. Storage requirements and computational costs to solve the system are thereby reduced. Stability and convergence of the method are proven under a timestep condition involving fluctuations. A series of numerical tests, including predictability horizons, are provided which confirm the theoretical analyses and illustrate uses of ensemble simulations.

## 1 Introduction

Ensemble calculations are essential in predictions of the most likely outcome of systems with uncertain data, e.g., weather forecasting [12], ocean modeling [14], turbulence [11], etc. Ensemble simulations classically involve J sequential, fine mesh runs or J parallel, coarse mesh runs of a given code. This leads to a competition between ensemble size and mesh density. We develop linearly implicit timestepping methods with shared coefficient matrices to address this issue. For such methods, it is more efficient in both storage and solution time to solve J linear systems with a shared coefficient matrix than with J different matrices.

Prediction of thermal profiles is essential in many applications [1, 7, 17, 16]. Herein, we extend [6] from isothermal flows to temperature dependent natural convection. We consider two natural convection problems enclosed in mediums with: non-zero wall thickness [3] and zero wall thickness; Figure 1 illustrates a typical setup. The latter problem is often utilized as a thin wall approximation.

Consider the Thick wall problem. Let be polyhedral domains in with boundaries and , respectively, such that dist(,) . The boundary is partitioned such that with and . Given and for , let , , and satisfy

 ut+u⋅∇u−PrΔu+∇p =PrRaγT+finΩf, (1) ∇⋅u =0inΩf, (2) Tt+u⋅∇T−∇⋅(κ∇T) =ginΩ, (3) u=0on∂Ωf,u=0inΩ−Ωf,T =0onΓ1andn⋅∇T=0onΓ2. (4)

Here denotes the usual outward normal, denotes the unit vector in the direction of gravity, is the Prandtl number, is the Rayleigh number, and in and in is the thermal conductivity of the fluid or solid medium. Further, and are the body force and heat source, respectively.

Let and . To present the idea, suppress the spatial discretization for the moment. We apply an implicit-explicit time-discretization to the system (LABEL:s1) - (LABEL:s1f), while keeping the coefficient matrix independent of the ensemble members. This leads to the following timestepping method:

 un+1−unΔt+n⋅∇un+1+u′n⋅∇un−Pr△un+1+∇pn+1 =PrRaγTn+1+fn+1, (5) ∇⋅un+1 =0, (6) Tn+1−TnΔt+n⋅∇Tn+1+u′n⋅∇Tn−κΔTn+1 =gn+1, (7)

Consider the Thin wall problem. The main difference is a “” term on the r.h.s of the temperature equation (LABEL:s2T) absent in (LABEL:s1T). This apparently small difference in the model produces a significant difference in the stability of the approximate solution. In particular, a discrete Gronwall inequality is used which allows for the loss of long-time stability; see Section 4 below. Consider:

 ut+u⋅∇u−PrΔu+∇p =PrRaγT+finΩ, (8) ∇⋅u =0inΩ, (9) Tt+u⋅∇T−∇⋅(κ∇T) =u1+ginΩ, (10) u=0on∂Ω,T=0onΓ1,n⋅∇T =0onΓ2, (11)

where is the first component of the velocity. If we again momentarily disregard the spatial discretization, our timestepping method can be written as:

 un+1−unΔt+n⋅∇un+1+u′n⋅∇un−Pr△un+1+∇pn+1 =PrRaγTn+fn+1, (12) ∇⋅un+1 =0, (13) Tn+1−TnΔt+n⋅∇Tn+1+u′n⋅∇Tn−κΔTn+1 =un1+gn+1, (14)

By lagging both and the coupling terms in the method, the fluid and thermal problems uncouple and each sub-problem has a shared coefficient matrix for all ensemble members.

Remark: The formulation (LABEL:d1a) - (LABEL:d2a) arises, e.g., in the study of natural convection within a unit square or cubic enclosure with a pair of differentially heated vertical walls. In particular, the temperature distribution is decomposed into , where is the linear conduction profile and satisfies homogeneous boundary conditions on the corresponding pair of vertical walls.

In Section 2, we collect necessary mathematical tools. In Section 3, we present algorithms based on (LABEL:d1a) - (LABEL:d2a) and (LABEL:d1) - (LABEL:d2). Stability and error analyses follow in Section 4. We end with numerical experiments and conclusions in Sections 5 and 6. In particular, two stable, convergent ensemble algorithms are presented. These algorithms can be used to efficiently compute an ensemble of solutions to (LABEL:s1) - (LABEL:s1f) and (LABEL:s2) - (LABEL:s2f) and estimate predictability horizons. The ensemble average is shown to produce a better estimate of the energy in the system, for a test problem, than any member of the ensemble.

## 2 Mathematical Preliminaries

The inner product is and the induced norm is . Define the Hilbert spaces,

 X :=H10(Ω)d={v∈H1(Ω)d:v=0on∂Ω},Q:=L20(Ω)={q∈L2(Ω):∫Ωqdx=0}, W :={S∈H1(Ω):S=0onΓ1},V:={v∈X:(q,∇⋅v)=0∀q∈Q}.

The explicitly skew-symmetric trilinear forms are denoted:

 b(u,v,w) =12(u⋅∇v,w)−12(u⋅∇w,v)∀u,v,w∈X, b∗(u,T,S) =12(u⋅∇T,S)−12(u⋅∇S,T)∀u∈X,∀T,S∈W.

They enjoy the following continuity results and properties.

###### Lemma 1.

There are constants and such that for all u,v,w X and T,S , and satisfy

 b(u,v,w) =∫Ωu⋅∇v⋅wdx+12∫Ω(∇⋅u)v⋅wdx, b∗(u,T,S) =∫Ωu⋅∇TSdx+12∫Ω(∇⋅u)TSdx, b(u,v,w) ≤C1∥∇u∥∥∇v∥∥∇w∥, b(u,v,w) ≤C2√∥u∥∥∇u∥∥∇v∥∥∇w∥, b∗(u,T,S) ≤C3∥∇u∥∥∇T∥∥∇S∥, b∗(u,T,S) ≤C4√∥u∥∥∇u∥∥∇T∥∥∇S∥, b(u,v,w) ≤C5∥∇u∥∥∇v∥√∥w∥∥∇w∥, b∗(u,T,S) ≤C6∥∇u∥∥∇T∥√∥S∥∥∇S∥.
###### Proof.

The proof of the first two identities is a calculation. The next four results follow from applications of Hölder and Sobolev embedding inequalities; see Lemma 2.2 on p. 2044 of [13]. We will prove the last two results for ; for they are improvable. For all u,v,w X,

 |(u⋅∇v,w)| ≤C∥u∥L6∥∇v∥∥w∥L3 ≤C∥∇u∥∥∇v∥√∥w∥∥∇w∥,

where Hölder, Ladyzhenskaya and Gagliardo-Nirenberg inequalities were used, respectively. Using the above result and inequalities and the first identity in Lemma LABEL:l1,

 |b(u,v,w)| =|(u⋅∇v,w)+12∫Ω(∇⋅u)v⋅wdx| ≤|(u⋅∇v,w)|+|12∫Ω(∇⋅u)v⋅wdx| ≤C∥∇u∥∥∇v∥√∥w∥∥∇w∥+C∥∇⋅u∥∥v∥L6∥w∥L3 ≤C∥∇u∥∥∇v∥√∥w∥∥∇w∥+C∥∇u∥∥∇v∥√∥w∥∥∇w∥ ≤C∥∇u∥∥∇v∥√∥w∥∥∇w∥.

In similar fashion, there is a such that

 |b∗(u,T,S)| ≤|(u⋅∇T,S)|+|12∫Ω(∇⋅u)TSdx| ≤C∥∇u∥∥∇T∥√∥S∥∥∇S∥+C∥∇⋅u∥∥T∥L6∥S∥L3 ≤C∥∇u∥∥∇T∥√∥S∥∥∇S∥+C∥∇u∥∥∇T∥√∥S∥∥∇S∥ ≤C∥∇u∥∥∇T∥√∥S∥∥∇S∥.

The weak formulation of system (LABEL:s1) - (LABEL:s1f) is: Find , , for a.e. satisfying for :

 (ut,v)+b(u,u,v)+Pr(∇u,∇v)−(p,∇⋅v) =PrRa(γT,v)+(f,v)∀v∈X, (15) (q,∇⋅u) =0∀q∈Q, (16) (Tt,S)+b∗(u,T,S)+κ(∇T,∇S) =(g,S)∀S∈W. (17)

Similarly, the weak formulation of system (LABEL:s2) - (LABEL:s2f) is: Find , , for a.e. satisfying for :

 (ut,v)+b(u,u,v)+Pr(∇u,∇v)−(p,∇⋅v) =PrRa(γT,v)+(f,v)∀v∈X, (18) (q,∇⋅u) =0∀q∈Q, (19) (Tt,S)+b∗(u,T,S)+κ(∇T,∇S) =(u1,S)+(g,S)∀S∈W. (20)

### 2.1 Finite Element Preliminaries

Consider a regular, quasi-uniform mesh of with maximum triangle diameter length . Further, for the system (LABEL:s1) - (LABEL:s1f), suppose that and lie along the meshlines of the triangulation of . Let , , and be conforming finite element spaces consisting of continuous piecewise polynomials of degrees j, l, and j, respectively. Moreover, assume they satisfy the following approximation properties :

 ≤Chk+1|u|k+1, (21) infqh∈Qh∥p−qh∥ ≤Chm|p|m, (22) ≤Chk+1|T|k+1, (23)

for all , , and . Furthermore, we consider those spaces for which the discrete inf-sup condition is satisfied,

 infqh∈Qhsupvh∈Xh(qh,∇⋅vh)∥qh∥∥∇vh∥≥β>0, (24)

where is independent of . The space of discretely divergence free functions is defined by

 Vh:={vh∈Xh:(qh,∇⋅vh)=0,∀qh∈Qh}.

The space , dual to , is endowed with the following dual norm

 ∥w∥V∗h:=supvh∈Vh(w,vh)∥∇vh∥.

The discrete inf-sup condition implies that we may approximate functions in well by functions in ,

###### Lemma 2.

Suppose the discrete inf-sup condition (LABEL:infsup) holds, then for any

 infvh∈Vh∥∇(v−vh)∥≤C(β)infvh∈Xh∥∇(v−vh)∥.

###### Proof.

See Chapter 2, Theorem 1.1 on p. 59 of [8].

We will also assume that the mesh and finite element spaces satisfy the standard inverse inequality [5]:

 ∥∇χ1,2∥≤Cinv,1,2(αmin)h−1∥χ1,2∥∀χ1∈Xh,∀χ2∈Wh,

where denotes the minimum angle in the triangulation. A discrete Gronwall inequality will play a role in the upcoming analysis.

###### Lemma 3.

(Discrete Gronwall Lemma). Let , H, , , , and be finite nonnegative numbers for n 0 such that for N 1

 aN+ΔtN∑0bn ≤ΔtN−1∑0dnan+ΔtN∑0cn+H,

then for all and N 1

 aN+ΔtN∑0bn ≤exp(ΔtN−1∑0dn)(ΔtN∑0cn+H).

###### Proof.

See Lemma 5.1 on p. 369 of [10].

The discrete time analysis will utilize the following norms :

 :=max0≤n≤N∥vn∥k,|||v|||p,k:=(ΔtN∑n=0∥vn∥pk)1/p.

## 3 Numerical Scheme

Denote the fully discrete solutions by , , and at time levels , , and . Given , find satisfying, for every , the fully discrete approximation of the Thick wall problem:

 (un+1h−unhΔt,vh)+b(n,un+1h,vh)+b(u′nh,unh,vh)+Pr(∇un+1h,∇vh)−(pn+1h,∇⋅vh)=PrRa(γTn+1h,vh)+(fn+1,vh)∀vh∈Xh, (25)
 (qh,∇⋅un+1h)=0∀qh∈Qh, (26)
 (Tn+1h−TnhΔt,Sh)+b∗(n,Tn+1h,Sh)+b∗(u′nh,Tnh,Sh)+κ(∇Tn+1h,∇Sh)=(gn+1,Sh)∀Sh∈Wh. (27)

Thin wall problem:

 (un+1h−unhΔt,vh)+b(n,un+1h,vh)+b(u′nh,unh,vh)+Pr(∇un+1h,∇vh)−(pn+1h,∇⋅vh)=PrRa(γTnh,vh)+(fn+1,vh)∀vh∈Xh, (28)
 (qh,∇⋅un+1h)=0∀qh∈Qh, (29)
 (Tn+1h−TnhΔt,Sh)+b∗(n,Tn+1h,Sh)+b∗(u′nh,Tnh,Sh)+κ(∇Tn+1h,∇Sh)=(un1,Sh)+(gn+1,Sh)∀Sh∈Wh. (30)

Remark: The treatment of the nonlinear terms in the time discretizations (LABEL:d1a) - (LABEL:d2a) and (LABEL:d1) - (LABEL:d2) leads to a shared coefficient matrix independent of the ensemble members.

## 4 Numerical Analysis of the Ensemble Algorithm

We present stability results for the aforementioned algorithms under the following timestep condition:

 C†Δthmax1≤j≤J∥∇u′nh∥2≤1, (31)

where . In Theorems LABEL:t1 and LABEL:t2, the nonlinear stability of the velocity, temperature, and pressure approximations are proven under condition (LABEL:c1) for the thick wall (LABEL:scheme:one:velocity) - (LABEL:scheme:one:temperature) and thin wall problems (LABEL:scheme:two:velocity) - (LABEL:scheme:two:temperature), respectively.

Remark: Stability of the numerical approximations can also be proven under: . If , then can be replaced with .

### 4.1 Stability Analysis

###### Theorem 4.

Consider the Thick wall problem (LABEL:scheme:one:velocity) - (LABEL:scheme:one:temperature). Suppose , . If (LABEL:scheme:one:velocity) - (LABEL:scheme:one:temperature) satisfy condition (LABEL:c1), then

 ∥TNh∥2+∥uNh∥2+12N−1∑n=0(∥Tn+1h−Tnh∥2+∥un+1h−unh∥2)+κΔt∥∇TNh∥2+PrΔt∥∇uNh∥2≤2ΔtPrRa2C2PF,1N−1∑n=0(Δtκn∑k=0∥gk+1∥2−1+∥T0h∥2+κΔt∥∇T0h∥2)+2ΔtPrN−1∑n=0∥fn+1∥2−1+∥u0h∥2+PrΔt∥∇u0h∥2+∥T0h∥2+κΔt∥∇T0h∥2.

Further,

 βΔtN−1∑n=0∥pn+1h∥≤2ΔtN−1∑n=0(C1∥∇n∥∥∇un+1h∥+C1∥∇u′nh∥∥∇unh∥+Pr∥∇un+1h∥+PrRaCPF,1∥Tn+1h∥+∥fn+1∥−1).

###### Proof.

Let in equation (LABEL:scheme:one:temperature) and use the polarization identity. Multiply by on both sides and rearrange. Then,

 12{∥Tn+1h∥2−∥Tnh∥2+∥Tn+1h−Tnh∥2}+κΔt∥∇Tn+1h∥2=Δt(gn+1,Tn+1h)−Δtb∗(u′nh,Tnh,Tn+1h). (32)

Use Cauchy-Schwarz-Young on ,

 Δt(gn+1,Tn+1h)≤Δt2ϵ∥gn+1∥2−1+Δtϵ2∥∇Tn+1h∥2. (33)

Consider . Add and subtract , use skew-symmetry, Lemma LABEL:l1, the inverse inequality, and the Cauchy-Schwarz-Young inequality. Then,

 |−Δtb∗(u′nh,Tnh,Tn+1h)| =|−Δtb∗(u′nh,Tnh,Tn+1h−Tnh)| (34) ≤ΔtC6∥∇u′nh∥∥∇Tnh∥√∥Tn+1h−Tnh∥∥∇(Tn+1h−Tnh)∥ ≤ΔtC6C1/2inv,2h1/2∥∇u′nh∥∥∇Tnh∥∥Tn+1h−Tnh∥ ≤C26Cinv,2Δt2h∥∇u′nh∥2∥∇Tnh∥2+14∥Tn+1h−Tnh∥2.

Using (LABEL:stability:thick:estg) and (LABEL:stability:thick:estbstar) in (LABEL:stability:thick) leads to

 12{∥Tn+1h∥2−∥Tnh∥2+∥Tn+1h−Tnh∥2}+κΔt∥∇Tn+1h∥2≤Δt2ϵ∥gn+1∥2−1+Δtϵ2∥∇Tn+1h∥2+C26Cinv,2Δt2h∥∇u′nh∥2∥∇Tnh∥2+14∥Tn+1h−Tnh∥2.

Let , add and subtract to the l.h.s. Regrouping terms leads to

 12{∥Tn+1h∥2−∥Tnh∥2}+14∥Tn+1h−Tnh∥2+κΔt2{∥∇Tn+1h∥2−∥∇Tnh∥2}+κΔt2∥∇Tnh∥2[1−2C26Cinv,2Δtκh∥∇u′nh∥2]≤Δt2κ∥gn+1∥2−1.

By hypothesis, . Thus,

 12{∥Tn+1h∥2−∥Tnh∥2}+14∥Tn+1h−Tnh∥2+κΔt2{∥∇Tn+1h∥2−∥∇Tnh∥2}≤Δt2κ∥gn+1∥2−1.

Sum from to and put all data on the right hand side. This yields

 12∥TNh∥2+14N−1∑n=0∥Tn+1h−Tnh∥2+κΔt2∥∇TNh∥2≤Δt2κN−1∑n=0∥gn+1∥2−1+12∥T0h∥2+κΔt2∥∇T0h∥2. (35)

Therefore, the l.h.s. is bounded by data on the r.h.s. The temperature approximation is stable.

We follow an almost identical form of attack for the velocity as we did for the temperature. Let in (LABEL:scheme:one:velocity) and use the polarization identity. Multiply by on both sides and rearrange terms. Then,

 12{∥un+1h∥2−∥unh∥2+∥un+1h−unh∥2}+ΔtPr∥∇un+1h∥2=−Δtb(u′nh,unh,un+1h)+ΔtPrRa(γTn+1h,un+1h)+Δt(fn+1,un+1h). (36)

Use the Cauchy-Schwarz-Young inequality on and and note that ,

 ΔtPrRa(γTn+1h,un+1h) ≤ΔtPr2Ra2C2PF,12ϵ∥Tn+1h∥2+Δtϵ2∥∇un+1h∥2, (37) Δt(fn+1,un+1h) ≤Δt2ϵ∥fn+1∥2−1+Δtϵ2∥∇un+1h∥2. (38)

Using skew-symmetry, Lemma LABEL:l1, the inverse inequality, and the Cauchy-Schwarz-Young inequality on leads to

 |−Δtb(u′nh,unh,un+1h)|≤C25Cinv,1Δt2h∥∇u′nh∥2∥∇unh∥2+14∥un+1h−unh∥2. (39)

Using (LABEL:stability:thick:estT), (LABEL:stability:thick:estf), and (LABEL:stability:thick:estb) in (LABEL:stability:thicku) leads to

 12{∥un+1h∥2−∥unh∥2+∥un+1h−unh∥2}+PrΔt∥