 # Numerical approximation of statistical solutions of the incompressible Navier-Stokes Equations

Statistical solutions, which are time-parameterized probability measures on spaces of square-integrable functions, have been established as a suitable framework for global solutions of incompressible Navier-Stokes equations (NSE). We compute numerical approximations of statistical solutions of NSE on two-dimensional domains with non-periodic boundary conditions and empirically investigate the convergence of these approximations and their observables. For the numerical solver, we use Monte Carlo sampling with an H(div)-FEM based deterministic solver. Our numerical experiments for high Reynolds number turbulent flows demonstrate that the statistics and observables of the approximations converge. We also develop a novel algorithm to compute structure functions on unstructured meshes.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The flow of a viscous, incompressible Newtonian fluid is described by the (incompressible) Navier-Stokes equations:

 ∂tu+divx(u⊗u)−νΔxu+∇xp=f,divxu=0.

Here is the velocity, is the pressure, which plays the role of a Lagrange multiplier to enforce the divergence-free constraint on the velocity, is the kinematic viscosity and represents the effects of external forces, for example, gravity and buoyancy. This system also needs to be supplied with initial and boundary conditions.

Navier-Stokes equations play a fundamental role in many applications and have been studied extensively for more than a century. The existence of global-in-time weak solutions, both in two and three spatial dimensions, can be traced back to the pioneering works of Leray Le1934 and Hopf Ho1951 . In two spatial dimensions, the uniqueness of weak solutions of NSE has been established. However, in three spatial dimensions, the question of uniqueness remains unanswered and has been designated as a Millennium Prize Problem by the Clay Mathematics Institute.

Keeping the uniqueness question aside, it is known that, for small viscosity (or, equivalently, large Reynolds number

), i.e. when inertial forces are much stronger than viscous diffusion, the fluid flow is very sensitive to initial conditions and characterized by chaotic motions. Therefore, in a deterministic framework, measurement errors in the problem data could have a drastic affect on the solutions of NSE, and so, it is not physically meaningful to describe these turbulent flows as individual solutions. On the other hand, there is ample experimental evidence that the statistical observables, for example, mean and variance, can be inferred reliably for turbulent flows. Hence, for practical applications, it is more useful to consider the evolution of NSE under uncertainties on the problem data in a suitable probabilistic framework. In their seminal work

Fo1972 ; FPr1976

, Foiaş and Prodi proposed the framework of statistical solutions of incompressible Navier-Stokes equations, in which, given a probability distribution on the initial velocity, velocity ensembles are evolved according to the NSE and described by a time-parameterized family of probability measures on the function space of initial velocity. The existence of such statistical solutions, both in two and three spatial dimensions, has been well established

Fo1972 ; FMRTe2001 . The global uniqueness in three spatial dimensions is an open problem. In two dimensions, statistical solutions are unique and they are defined as the push-forward of the probability measure on the initial velocity data. The authors in FMRTe2001 show that several results of the conventional theory of turbulence for NSE attributed to the ground-breaking work of Kolmogorov Ko1941a ; Ko1941b can be recovered with statistical solutions, thus, providing evidence of the importance of this solution framework.

The computation of statistical solutions of NSE is a special case of Uncertainty Quantification (UQ) in Computational Fluid Dynamics (CFD), for which different methods are available in the literature (see BLMSc2013 and the references therein). Loosely speaking, the different methods used to solve UQ problems can be categorized into two classes: stochastic Galerkin methods and stochastic collocation methods. Stochastic Galerkin methods consider the Ritz-Galerkin formulation of the underlying PDEs in the stochastic parameter space, thus, these methods are highly intrusive and not amenable for implementation from a practical viewpoint (Su2015, , Chap. 12). On the other hand, stochastic collocation methods are well-suited for large scale applications as they are non-intrusive (Su2015, , Chap. 13). For applications, Monte Carlo (MC) and Quasi-Monte Carlo sampling can be seen as a subclass of stochastic collocation methods. Multi-level Monte Carlo methods have also received a lot of attention in research in recent times (see LMSc2016 and the references therein).

In the discussion above, the non-intrusive methods require a deterministic solver that numerically approximates the solution of the underlying system of PDEs (in our case the NSE). There are many numerical methods available in the literature that can efficiently approximate the NSE. For problems with periodic boundary, spectral methods are appealing because of their efficiency and high resolution Le2018 . Finite difference methods with Leray projection were first introduced in KLa1966 ; Ch1967 ; Ch1969 . Later, the authors in BCGl1989 proposed a finite-volume scheme that was an efficient variant of Ch1969 . In LMSc2016 , the authors develop a finite-difference scheme for the vorticity formulation of the NSE with uniformly stable bounds with respect to viscosity; this scheme is a variant of the one proposed in LTa1997 for incompressible Euler equations. Finite element methods are another class of methods to solve NSE GRa1986 . While standard mixed FEM are not pressure-robust, i.e. pressure approximation influences the velocity approximation, (weakly) divergence-free mixed methods are (see JLMNRe2017 and the references therein). In these weakly divergence-free mixed methods, if the velocity field belongs to the Hilbert space

, then the a priori error bounds are non-uniform with respect to viscosity. However, if H(div)-conforming spaces are used instead, then the a priori velocity error estimates are robust with respect to viscosity

SLe2018 . Moreover, the divergence-free constraint on the velocity field is satisfied point-wise in the spatial domain. The authors in LSc2016 proposed a computationally efficient hybridized variant of H(div) FEM, called the H(div)-HDG method, where they introduce hybridization for the tangential velocity components. This method was further improved in LLSc2018 and used to study turbulent flows in FKLLSc2019 ; SJLLLSc2019 .

### 1.1 Contributions

We restrict ourselves to two-dimensional open, bounded spatial domains and choose the framework of statistical solutions for NSE. Given the discussion above, to numerically approximate these solutions, we use Monte Carlo sampling for the stochastic space and H(div)-FEM as the deterministic solver. We can easily use QMC instead of MC, but it is not necessary as our main focus is on the presence or absence of convergence of the approximate solutions and their statistics, not the rates. We choose the H(div) scheme and not the more efficient H(div)-HDG scheme because of the ease of implementation. In the deterministic solver, we use the implicit Euler method for time integration as there is no CFL restriction on the time-step size (see MPRo1990 ; LSc2016 for higher-order time integration methods). To the best of my knowledge, this is the first computational effort in this direction.

Structure functions are an important observable in turbulence literature Fr1995 . In the context of incompressible Euler equations, a uniform decay of the structure functions ensures the convergence of approximate statistical solutions in the very recent work LMPa2021a . The behaviour of structure functions is also crucial in the study of the vanishing viscosity limit of incompressible Navier-Stokes in LMPa2021b . In this light, a novel contribution of this work is the development of an algorithm for approximating structure functions on unstructured meshes. We compute the structure functions only at the final time, which is a good indicator of their behaviour over the whole time interval as seen in LMPa2021a .

### 1.2 Outline

After covering some preliminaries, we state the initial boundary value problem (IBVP) and introduce the concept of weak solutions for NSE. In §4, we describe statistical solutions of NSE and define structure functions and Wasserstein distances. We present the deterministic H(div)-conforming scheme in §5. In §6, we describe the algorithm for approximating structure functions on unstructured meshes. We present the results of our numerical experiments in §7 and conclusions in §8.

## 2 Preliminaries

Consider an open time interval , where denotes a finite time horizon, and let be an open and bounded polygon with Lipschitz boundary and with outward unit normal . Define the space-time cylinder as the Cartesian product of the spatial domain and the time domain , i.e. . Let the scalar

and the vector

denote the time coordinate and the spatial Cartesian coordinates, respectively, then the space-time Cartesian coordinates are .

We highlight vector fields by bold face fonts and tensor fields by bold face fonts with an underline.

Given and a vector , the norm , we also denote the Euclidean norm by .

Given and two matrices with entries , and , we denote by the Frobenius inner product of the two matrices.

### 2.1 Function spaces

For any open, bounded domain , we use the Lebesgue spaces for scalar-valued functions with the associated norm , for . We also use the standard Hilbert spaces with associated norms and seminorms , for . For vector-valued functions of size , we indicate these spaces as and . Spaces and norms for tensor-valued functions are indicated with bold font. We use the standard notation for the space of continuous functions on and , , for the space of functions on that are -times differentiable with continuous -th derivative. The space

 H(div;D) :={τ∈L2(D)2:∇x⋅τ∈L2(D)}. (2.1)

For a Hilbert space and , we use the standard notation for Bochner spaces and define .

### 2.2 Probability

Let be a probability space Sc2005

. For any random variable

, its mean and variance are given by

 E[X] =∫ΩX(ω)dP(ω), (2.2a) Var[X] =E[(X−E[X])2]=E[X2]−E[X]2. (2.2b)

We assume that is complete.

Given a topological space , we denote by the Borel -algebra of and by the space of all probability measures on .

### 2.3 Meshes, mesh faces, averages and jumps

We choose a mesh of domain such that it is shape-regular and does not allow hanging nodes. We allow both simplicial meshes and meshes with quadrilateral elements, which we also refer to as triangular and quadrilateral meshes, respectively. The mesh-size of is denoted by .

We define the following sets and unions for mesh faces:

• denotes the set of all faces of ,

• denotes the set of all interior faces of ,

• denotes the set of all boundary faces of ,

• the mesh skeleton ,

• the interior mesh skeleton ,

• the boundary mesh skeleton .

Further, we divide into two disjoint sets and that correspond to the boundaries with Dirichlet and outflow boundary conditions, respectively, in NSE (3.1). Using these sets, we also define the unions and .

For , we denote the outward-pointing unit normal vector on its boundary by . Let be distinct. Then, for a vector field and a tensor field , which are element-wise continuous on , we define the averages and jumps on the interior mesh face as follows:

 {{w}} :=w|K1+w|K22, (2.3) {{τ––}} :=τ––|K1+τ––|K22, [[w⊗nF]] :=w|K1⊗nK1+w|K2⊗nK2, [[w⋅nF]] :=w|K1⋅nK1+w|K2⋅nK2, [[[w]]] :=w|K1−w|K2.

Here is the unit normal at face and we set by convention. We extend the definitions of the jumps and averages for boundary faces. Let such that , then we define

 {{w}} :=w|K, {{τ––}} :=τ––|K, (2.4) [[w⊗nF]] :=w|K⊗nK, [[w⋅nF]] :=w|K⋅nK.

We frequently drop the subscript in for convenience.

### 2.4 Raviart-Thomas spaces

We revisit the definition of Raviart-Thomas spaces and some important results associated with them, which are presented in (Ga2014, , Chapter 3) and (BBFo2013, , Chapter 2). We use these spaces to construct element-wise polynomial subspaces of in our numerical scheme.

Given , we denote by the space of polynomials of total degree at most , and we define the discrete polynomial space

 (2.5)

We denote Raviart-Thomas spaces of degree by , they are defined as follows:

###### Definition 2.1 (Raviart-Thomas spaces).

For any element , the local Raviart-Thomas space of degree is defined as

 RTk(K):=Pk(K)2+Pk(K)x. (2.6)

Let denote the global Raviart-Thomas space, it is defined as

 RTk(Th):={v∈H(div;D):v|K∈RTk(K) ∀K∈Th}. (2.7)
###### Lemma 2.2 ((BBFo2013, , §2.5.2)).

Any function is continuous in the normal direction at all interior faces of the mesh, i.e.

 [[v⋅n]]=0, ∀F∈Fint. (2.8)

Moreover, the following holds true:

 divx(RTk(Th))=Qk(Th). (2.9)

Consider the equation , for and for all , which appears in the numerical scheme we use in §5.1. From Lemma 2.2, we deduce that in . Thus, the divergence-free constraint in NSE can be imposed easily with the pair of FE spaces and .

## 3 Navier-Stokes equations

We assume that the boundary can be divided into Dirichlet boundary and outflow boundary , where either or may be empty. We also assume that and consist of entire segments of .

The IBVP for the time-dependent incompressible Navier-Stokes equations in the divergence form reads as follows:

 ut+divx(u⊗u)−νΔxu+∇xp =f in Q, (3.1a) divxu =0 in Q, (3.1b) u(⋅,0) =u0 in D, (3.1c) u =g on ∂DD×I, (3.1d) (ν∇xu−pI–)⋅nD =0 on ∂Dout×I, (3.1e)

where is the velocity, is the pressure, is the kinematic viscosity, is the initial velocity, is the forcing, is the boundary data and is the identity tensor.

###### Remark 3.1.

For very large Reynolds numbers, a stabilized outflow boundary condition is needed to ensure numerical stability in the case of back-flow at the outlet, we refer to BCBBr2018 for a review of the different stabilization methods available in the literature and in particular to DKCh2014 ; Do2015 for some popular choices. In the present work, see §7, our numerical experiments for problems with outflow preclude any backward-flow at the outlet. Therefore, the outflow condition (3.1e) is sufficient for our purpose.

### 3.1 Weak solutions

We assume homogeneous Dirichlet BCs, i.e. and . We choose the functional setting used in (FMRTe2001, , §II.5) and define the following Hilbert spaces:

 H:=Hnsp :={v∈L2(D)2:divxv=0, v⋅nD|∂D=0}, (3.2) V:=Vnsp :={v∈H1(D)2:divxv=0, v|∂D=0}. (3.3)

On these spaces, we define the inner products

 (v,w)H :=∫Dv⋅w dx, for v,w∈H, (3.4) (v,w)V (3.5)

and the associated norms

 (3.6)

The weak formulation reads as: find such that, ,

 ddt(u,v)H+ν(u,v)V+(u⋅∇xu,v)L2(D)=(f,v)H, (3.7)

with the viscosity , the initial velocity and the forcing . This weak formulation can be traced back to the seminal work of Leray Le1933 ; Le1934 .

It is well-known that in two-dimensional spatial domains with sufficiently smooth boundaries, for any , any initial velocity and any forcing , weak solutions of NSE exist and they are unique. We state this result in the following (cf. (FMRTe2001, , Chapter II, Theorem 7.3 and Remark 7.2)):

###### Theorem 3.2 (Existence and uniqueness of weak solutions of NSE).

Assume that the spatial domain has boundary , and that

 u0∈H,f∈L2(I;H). (3.8)

Then, for every , there exists a unique solution of (3.7) such that

 u ∈L2(I;V)∩C(¯¯¯¯I;H), (3.9)

and, for all , it satisfies the energy equation

Moreover, there exists a continuous solution operator such that .

## 4 Statistical solutions

The discussion in this section is based on the monograph (FMRTe2001, , Chapter V) and it partially uses the notation prescribed in (LMSc2016, , §2.2). We consider statistical solutions of NSE (3.1) with fixed, no-slip boundary, i.e. and , where the boundary is of class .

Given a probability distribution on the initial velocity data, i.e.

 P({u0∈E})=μ0(E), for E∈B(H),

the main idea of statistical solutions, as introduced by Foiaş-Prodi Fo1972 ; FPr1976 , is to describe the evolution of velocity ensembles by a time-parameterized family of probability measures , such that (up to null sets), and for every , . This solution framework includes the individual weak solutions of NSE as a special case with the probability measure , where is the Dirac measure (also known as the unit mass) at , cf. (Sc2005, , §4.7).

From Theorem 3.2, we know that the individual weak solutions of NSE are unique and the solution operator exists. As a result, for any time-independent forcing , unique solutions are given by transporting the initial probability measure under the solution operator , i.e.

 μνt(E)=μ0((Sν(t))−1E), ∀E∈B(H). (4.1)

Concretely, we state this existence and uniqueness result in the following (cf. (FMRTe2001, , Chapter V, Theorems 1.1 and 1.2)):

###### Theorem 4.1.

Let be a probability measure on with finite kinetic energy, i.e.

Assume that the forcing term . Then, for every and homogeneous Dirichlet boundary conditions, there exists a statistical solution of the incompressible Navier-Stokes equations (3.1) on .

Moreover, if has bounded support in and is independent of time, then the statistical solution is unique and is given by , where is the solution operator of the incompressible Navier-Stokes equations, cf. Theorem 3.2.

### 4.1 Structure functions

For a solution , other than its common statistics like mean and variance, its structure functions are of great interest to us. The structure functions are defined as follows (we refer to Ly2020 ):

###### Definition 4.2 (Structure functions).

Let , and . Let be a random field such that , for any , and be the probability distribution of at time . Then, the structure function associated with is given by:

 Spr,t(μt)=(∫Lp(D)2∫D\fintBr(x)∥u(⋅;x,t)−u(⋅;y,t)∥ppdy dx dμt(u))1p, (4.2)

where denotes the ball of radius centered at a point .

We use the structure functions as a measure of the regularity of the velocity field , and we investigate the scaling behaviour of with respect to the offset in our numerical experiments.

###### Remark 4.3.

If is Lipschitz-continuous in , then straightforward calculations show that . We expect the same scaling behaviour for the velocity field in NSE due to regularity (3.9) under the assumptions of Theorem 3.2.

### 4.2 Wasserstein distances

In the next section, we describe a numerical method to approximate statistical solutions with ensembles of approximate individual velocity solutions. To study the convergence of these velocity ensembles, it is important to define a notion of distance between two ensembles. For this purpose, we use Wasserstein distance (cf. (FLMi2017, , Definition 2.2)):

###### Definition 4.4.

Let be a separable Banach space, and let be probability measures on with finite

-th moments, i.e.

and . Then, the -Wasserstein distance between and is defined as

 Wp(μ,ρ):=(infγ∈Γ(μ,ρ)∫X2|x−y|p dγ(x,y))1p,

where the infimum is taken over the set of all transport plans from to , i.e. those that satisfy

 ∫X2(F(x)+G(y))dγ(x,y)=∫XF(x)dμ(x)+∫XG(y)dρ(y),∀F,G∈Cb(X).

Here denotes the space of bounded, continuous, real-valued functionals on .

Given , we can define -Wasserstein distance for the -point distribution at and the -point correlation at , respectively, as:

 Wp(μ(x),ρ(x)) :=(infγ∈Γ(μ,ρ)∫H2∥u(x)−v(x)∥p2dγ(u,v))1p, (4.3) Wp(μ(x,y),ρ(x,y)) :=(infγ∈Γ(μ,ρ)∫H2∥∥∥[u(x)u(y)]−[v(x)v(y)]∥∥∥p2dγ(u,v))1p. (4.4)

For the whole domain , we can define

 W1p(μ,ρ) :=∫x∈DWp(μ(x),ρ(x)) dx, (4.5) W2p(μ,ρ) :=∫(x,y)∈D2Wp(μ(x,y),ρ(x,y)) dydx. (4.6)

In this work, we always use -Wasserstein distances and , so we simply refer to them as Wasserstein distances and drop the subscripts.

## 5 Numerical approximation of statistical solutions

In this section, our goal is to approximate the statistical solutions described in the previous section. To this end, we assume that initial probability measure is given as the law of a random field defined on the underlying probability space . This allows us to draw random samples from the initial distribution, which can then be evolved using a numerical solution operator. We estimate the statistical solutions using (5.1). Concretely, this Monte Carlo type method is described in Algorithm 1.

Here is a numerical approximation of the operator such that .

To approximate an individual solution, we choose an H(div)-conforming method based on the work of Sequeira et al. GSSe2017 and Cockburn et al. CKSc2007 ; this scheme has also been considered in SLe2018 ; SLu2018 .

For spatial discretisation §5.1, we use the H(div)-conforming scheme designed in GSSe2017 for the incompressible Euler equations along with the symmetric interior penalty discretisation for the viscous diffusion given in CKSc2007 and (DEr2012, , §4.2.2). This semi-discrete scheme is pressure-robust and -semi-robust, we refer to SLe2018 for details, which means that the velocity error bounds have no explicit dependence on the pressure approximation and the Reynolds number. To obtain a fully discrete scheme §5.2, we combine the H(div) spatial discretisation with implicit Euler time-stepping as formulated in GSSe2017 .

### 5.1 H(div)-conforming spatial discretisation

We use the function space , which is defined as

 Vk:=Vk(Th):={v∈RTk(Th):v⋅nD=0 on ∂DD}, (5.2)

for the velocity field. Here the space is given by (2.7). In the case of purely Dirichlet boundary, we denote by .We use the function space (2.5) for the pressure field.

###### Remark 5.1.

For purely Dirichlet boundary, i.e. , if is a solution of (3.1) such that , then , for any constant , is also a solution. Therefore, an additional constraint is needed to ensure a unique pressure solution. A common choice is to impose vanishing mean for the pressure (cf., for example, GSSe2017 ; JLMNRe2017 ) using the function space

 Q0,k:=Q0,k(Th) :={q∈L20(D):q|K∈Pk(K) ∀K∈Th}, (5.3a) where L20(D) :={q∈L2(D):∫Dq dx=0}. (5.3b)

Given initial velocity data , let denote the -orthogonal projection of into . Then, the H(div)-conforming spatial discretisation of the IBVP (3.1) reads as follows:

Given the initial data , the forcing and the boundary data , find and such that, for all and for all ,

 =ℓh(vh), (5.4a) Bh(uh,qh) =0, (5.4b)

where the mass bilinear form

 Mh(u,v):=∫Thu⋅v dx,for u,v∈Vk, (5.5)

the convection trilinear form with upwind flux

 Cuph(w;u,v) := ∑K∈Th∫K(w⋅∇xu)⋅v dx−∫Finth(w⋅n)[[[u]]]⋅{{v}} dx (5.6) +∫Finth|w⋅n|[[[u]]]⋅[[[v]]] dx,for w,u,v∈Vk,

where (cf. (2.8)), the symmetric interior penalty diffusion bilinear form

 Asiph(u,v) :=∑K∈Th∫K∇xu:∇xv dx (5.7) −∑F∈Fint∪Fbdr,D∫F[[v⊗n]]:{{∇xu}} dx −∑F∈Fint∪Fbdr,D∫F[[u⊗n]]:{{∇xv}} dx +∑F∈Fint∪Fbdr,Dσh−1F∫F[[u⊗n]]:[[v⊗n]] dx,for u,v∈Vk,

where is the jump penalization parameter and is the size of face , the divergence bilinear form

 Bh(u,q):=∑K∈Th∫Kqdivx(u) dx,for u∈Vk, q∈Qk, (5.8)

and the linear form

 ℓh(v) :=∫Thf⋅v dx (5.9) −ν∑F∈Fbdr,D∫F(g⊗n):∇xv ds +ν σ∑F∈Fbdr,Dh−1F∫F(g⊗n):(v⊗n) dx,for v∈Vk.
###### Remark 5.2.

In the formulation (5.4), the outflow boundary conditions are imposed weakly through the boundary integrals on that appear in the bilinear forms and .

###### Remark 5.3.

For homogeneous Dirichlet boundary conditions, the scheme (5.4) is the same as the one analyzed in SLe2018 , except for the slight difference in the convection trilinear term that does not affect the error estimates. Therefore, for , from Theorem 5.3, Theorem 5.6 and Corollary 5.9 in SLe2018 , we deduce that the error for the discrete solution of (5.4) is bounded by , , for with sufficient regularity.

### 5.2 Fully-discrete scheme

As mentioned, we combine the semi-discrete formulation (5.4) with implicit Euler time-stepping to obtain a fully discrete scheme §5.2, by following the procedure described in GSSe2017 .

For , let be a set of points in the time interval such that with . Given the solution of the IBVP (3.1), we define the notation and , and we denote by and the numerical approximations of and , respectively.

The Taylor series expansion of at is given by

 u(⋅,tn−1) =u(⋅,tn)−Δt∂tu(⋅,tn)+(Δt)22∂ttu(⋅,tn)+… =u(⋅,tn)+O(Δt), (5.10)

which yields the time derivative

 ∂tu(⋅,tn) =u(⋅,tn)−u(⋅,tn−1)Δt+E0(tn), (5.11)

with the truncation error

 E0(tn) :=Δt2∂ttu(⋅,tn)+…=O(Δt). (5.12)

Substituting (5.11) in (3.1), at any discrete time , we obtain

 un+Δt(un⋅∇xun−νΔxun+∇xpn) =un−1+Δt f(⋅,tn)+Δt E0(tn), divxun =0,

which can be written equivalently as

 =un−1+Δt f(⋅,tn) (5.13) +Δt((un−1−un)⋅∇xun+E0(tn))=: ~E0(tn), divxun =0.

Using (5.2) and (5.12), we deduce that the term . Dropping the truncation term in (5.13), the remaining equation is linear with respect to the variable . We introduce H(div)-conforming spatial discretisation in the truncated system (5.13) to arrive at a fully discrete method, which reads as follows:

Given the initial data , the forcing and the boundary data , find and such that, for all and for all ,

 Mh(unh,vh) =Mh(un−1h,vh)+Δt ℓh(vh), (5.14a) Bh(unh,qh) =0, (5.14b)

for , with .

###### Remark 5.4.

The discrete formulation (5.14) satisfies the following:

 ∥∥unh∥∥2L2(D)2+ Δt(∣∣unh∣∣2un−1h,up+νCσ∥∥unh∥∥2e) ≤(un−1h,unh)L2(D)2+Δt(f,unh)L2(D)2,

for . Then, for the forcing , we have

 ∥∥unh∥∥L2(D)2≤∥∥un−1h∥∥L2(D)2.

Therefore, the scheme is -stable.

###### Remark 5.5.

For viscosity , to the best of my knowledge, error estimates of our numerical scheme are not available, and it is beyond the scope of this work to derive them. However, based on the reasoning in the proofs of (SLu2018, , Lemma 5.5, Theorem 5.6) and (SLe2018, , Lemma 3.4, Theorem 3.6), we expect the velocity error to be bounded by .

## 6 Implementation details

In this section, we discuss important details about our implementation of the Monte Carlo algorithm 1, which we refer to as MC-FEM solver; here all the samples are evolved in parallel using the self-scheduling algorithm (GLSk2014, , §3.6). In particular, we describe the deterministic solver that is used for a single Monte Carlo sample and present the algorithms used to approximate structure functions and Wasserstein distances.

We denote vectors by a lower-case font with an underline and matrices by upper-case font or bold font with an underline, only for the discussion on the deterministic solver. Given bases of and , we denote the corresponding vectors of basis coefficients by and . Then, the finite element assembly for the discrete formulation (5.14) yields the following block-structured linear system of equations:

 [Mh+Δt(Ch+νAh)ΔtB⊤hBh0––][u––p–]=[f––h0–], (6.1)

where the matrices , , and are assembled from the mass form , the convection form , the diffusion form and the divergence form , respectively. Here the right-hand side vector is assembled from the linear form . Note that we must use the Piola transformation in our finite element assembly.

We need to assemble the matrices , and only once, but the matrix has to be re-assembled for each discrete time because of its dependence on the velocity solution of the previous time step . Using direct solvers for the linear system (6.1) would involve factorizing the re-assembled system matrix at each time step, which is computationally expensive. Therefore, we use suitable iterative solvers with preconditioning. Re-scaling the pressure solution as