Unstable systems have finite predictability horizons, Lorenz [1,2]. The chaotic nature of fluid flow and the uncertainties in initial conditions limit predictability. Under different initial conditions, the trajectories of the flow spread. Small errors in the (uncertain) initial conditions can grow exponentially until O (1), resulting in a loss of prediction ability [3].
Ensemble methods address the uncertainty in problem data by conducting numerical simulations with various initial and boundary conditions, external forces, viscosity, and other model parameters, Kalnay [4]. Ensemble mean takes into account the combined information from multiple simulations, effectively smoothing out individual variations and providing a more robust prediction, hence it is considered the best estimate. Leith [5] showed that using a sample size of eight with Monte Carlo sampling can obtain adequate accuracy for the ensemble mean.
Ensemble methods address data uncertainty, thus improving predictability. One remaining issue is that solving J separate linear systems of equations for an ensemble of size J can be computationally expensive. Our method resolves this by introducing a shared coefficient matrix, allowing the system to be assembled once and solved efficiently for different right-hand side vectors.
We develop a penalty-based ensemble method for the Navier-Stokes equations (NSE) to reduce computational cost. This enables larger ensemble size, which produces more reliable results. The method uses a shared coefficient matrix with different right-hand side vectors and relaxes the incompressibility condition to reduce the space complexity of the model while maintaining accuracy. Further, eliminating the J pressure variables saves memory and operations. In this report, we derive a stability proof and an error estimate and conduct three numerical tests to validate the method. In addition, we extend the method from deterministic NSE to NSE with random variables, following an approach similar to that of Luo and Wang [6] for random parabolic partial differential equations.
The incompressible NSE is given by ∂u ∂t + u • ∇u -ν u + ∇ p = f (x, t), and ∇ • u = 0, u = 0 on the boundary.
(
where u denotes the flow velocity and p denotes the flow pressure. The viscosity is denoted by ν, and f is the body force. In (1.1), the pressure is a Lagrange multiplier to enforce the incompressibility constraint, E and Liu [7]. The penalty method relaxes incompressibility by replacing ∇ • u = 0 with ∇ • u + p = 0, for 0 < 1.
For the continuous problem, the pressure p can be replaced by p = -1 ∇ • u to replace p, which uncouples u and p and yields the penalized NSE:
We adopt the ensemble approach of Nan and Layton [8] to the penalized NSE, using a shared coefficient matrix with different right-hand sides. We suppress the spatial discretization to present the idea. We define the ensemble mean and fluctuation at the timestep t n :
u ,n j , and U ,n j := u ,n ju n , where u ,n j is the penalized velocity for the j th ensemble member. We use an implicitexplicit time discretization, which allows the coefficient matrix to be independent of the ensemble member. The method is to find u ,n+1 j in the velocity space and p ,n+1 j in the pressure space:
Here, is the same for all ensemble members to ensure a shared coefficient matrix. The ensemble mean drives the flow. We can eliminate the pressure by setting
to reduce the memory.
Epstein [9] introduced the first forecasting method that explicitly accounted for the uncertainty in atmospheric model predictions, known as the stochastic-dynamics forecasting method, in 1969. Leith [5] later proposed using ensemble forecasting with multiple members instead of a single realization. He showed that the ensemble mean from Monte Carlo ensembles can achieve accurate results without linear regression. Luo and Wang [6] studied an ensemble algorithm for the deterministic and random parabolic partial differential equations, which led to a single discrete system with multiple right-hand side vectors. Temam [10] first introduced the penalty method with a modified nonlinear term to ensure energy dissipation. He proved in [10] that lim →0 (u , p ) = (u, p). Penalty methods have been widely studied, including Falk [11], Shen [12] and He [13], He and Li [14]. We can speed up the calculation by eliminating the pressure by p = -1 ∇ •u , Heinrich and Vionnet [15]. The velocity error depends on the penalty parameter , as shown by Bercovier and Engelman [16]. The condition number of the penalized system was studied in Layton and Xu [17], Hughes and Liu and Brooks [18]. Adapting penalty parameters, exploiting -sensitivity, can help with ill-conditioning and provide better accuracy [19][20][21][22], and pressure recovery in [23]. Preliminary tests of the penalty-based ensemble method are studied in Fang [24].
Let D ⊂ R d be an open regular domain, where d = 2 or 3. The L 2 (D) norm is denoted as • , and the inner product is denoted as (•, •). Similarly, we define the L p (D) norms
• L p , and the Sobolev W k p (D) norms • W k p . We denote the Sobolev space W k 2 (D) with norm • k as H k (D). We define the norms for the functions v(x, t) defined on
We denote the discrete-time equivalents of the norms as follows:
Let ( , F, P) be a complete probability space, where is the set of outcomes, F ⊂ 2 is the σ -algebra of events, and P : F → [0, 1] is a probability measure. We denote the set of all integrable functions for the probability measure P as the L 1 P ( ). Suppose a random variable Y such that Y ∈ L 1 P ( ), we define the expected value of Y as follows:
The stochastic Sobolev spaces are denoted by
W k p contains stochastic functions v : × D → R, that are measurable with respect to the product σ -algebra F ⊗ B(D), where B is a Borel set. W k p is equipped with the averaging norms
When p = 2, the above space is a Hilbert space and we have
Lemma 2.1 (See Layton [25], p. 28, p. 29
Then, there is a positive constant C P F such that
Thus, ∇v and v are equivalent norms on H 1 0 (D).
The space H -k (D) is the dual space of bounded linear functionals on H k 0 (D). A norm for H -1 (D) is given by
Let X be the velocity space and Q be the pressure space:
where
We denote the conforming velocity and pressure finite element spaces as follows:
We assume that (X h , Q h ) satisfies the following approximation properties and the Ladyzhenskaya-Babushka-Brezzi condition (L B B h ). For u ∈ H m+1 (D) d and p ∈ H m (D), inf
(2.7)
where β h is bounded away from zero uniformly in h.
We assume the mesh with quasi-uniform triangulation and finite element spaces satisfy the inverse inequality:
(2.9) Lemma 2.4 (See Ladyshenskaya [26]) For any vector function u : R d → R d with compact support and with finite L p norms:
.10) Lemma 2.5 (A discrete Gronwall lemma, see Lemma 5.1, p. 369, [27]) Let t, B, a n , b n , c n , d n be nonnegative numbers such that for l ≥ 1:
then for all t > 0,
(2.12) Lemma 2.6 (See Layton [25], p. 7, Hölder's and Young's inequalities) For any ξ > 0, 1 ≤ p < ∞, and 1 p + 1 q = 1, the Hölder and Young's inequalities:
The generalization for three functions,
Lemma 2.9 (See Layton [25] and Girault and Raviart [28]) b * (u, v, w) satisfies the following bounds:
for all u, v, w ∈ X.
We define the final time T and timestep size at the n th step t n . The total number of steps N is given by N = T / t. The fully-discrete approximation is then given
for all
2) is satisfied, we proceed to the next step. Otherwise, we halve the timestep size and repeat the current step. The penalty-based ensemble method is a one-step method, which allows us to assume a constant t in the following analysis for simplicity. However, the time step can be adapted as needed in numerical tests.
We can replace the pressure with p j,h = -
This reformulation eliminates the pressure variable, reducing computational complexity and speeding up calculations in numerical tests.
Let the difference between the ensemble member j and the ensemble average be denoted as
In Theorem 3.1, we prove the method is nonlinearly and long-time energy stable under the CFL condition:
Theorem 3.1 Suppose the following timestep condition holds:
It yields that for any N ≥ 1:
(3.5)
Proof Taking the inner product of the momentum equation with
Set v h = u ,n+1 j,h . Multiply t to both sides of (3.6) and apply the polarization identity:
Next, we treat the trilinear term with the help of inverse inequalities and interpolation,
(3.9)
Combine terms, we have
(3.10)
Add and subtract ν t 4 ∇u ,n j,hfoot_0 , we have
With the CFL condition in (3.4), (3.11) reduces to:
2 )
(3.12) Sum over all n from 0 to N -1, we have the final result. (3.4), the following holds:
The fluctuation of the flow as defined in Jiang, Kaya, and Layton [29], Definition 2.1, is:
Sum from n = 0 to n = N , and multiply by t:
is bounded by a finite number. By Theorem 3.1, for j = 1, . . . , J , we have
.15) Proposition 3.4 (See John [30], p. 164, Lemma 4.43) Let the domain D be bounded with polyhedral and Lipschitz continuous boundary and (u, p) ∈ (X , Q). Suppose L B B h Condition 2.3 holds, then it yields
Denote the error of the j th simulation at time t n , e ,n j := u ,n ju ,n j,h . Here, u ,n j is the solution of the penalized NSE at time t n and u ,n j,h is the fully discretized solution of penalty-based ensemble method. Theorem 3.5 Consider the method in (1.3) and assume the condition in (3.4) holds for all n:
then there are positive constants C and C 0 independent of h and t such that:
where
Proof We evaluate the continuous penalty-based NSE (1.2) at time t = t n+1 . For any v h ∈ X h , and
Subtract equation (3.1) from (3.18). We have
(3.20)
Let ũ ∈ X h and q ∈ Q h satisfy the Stokes projection:
(3.21)
and q h = p ,n+1 j,h -q, then apply the polarization identity. We have
We bound the terms on the right-hand side.
By the integral form of Taylor's theorem, we have
t n η j,t ds.
Divided by t on both sides, and take the L 2 norm on D,
By Fubini's theorem, we have
Thus, we have
(r ,n+1 j , φ ,n+1 j,h ) term:
. By the integral form of Taylor's theorem:
j,t -t n+1 t n u j,tt (t ns) ds, r ,n+1 j = 1 t t n+1 t n u j,tt (st n ) ds. (3.24) 123 r ,n+1 j 2 = D 1 t t n+1 t n u j,tt (st n ) ds 2 dx ≤ D t n+1 t n u j,tt ds 2 dx ≤ D t n+1 t n |u j,tt | ds 2 dx ≤ D t n+1 t n 1 ds t n+1 t n |u j,tt | 2 ds dx = t D t n+1 t n |u j,tt | 2 ds dx. (3.25) By Fubini's theorem, we have r ,n+1 j 2 ≤ t t n+1 t n D |u j,tt | 2 dx ds. Hence (r ,n+1 j , φ ,n+1 j,h ) ≤ C(ν)( t) 2 t n+1 t n D |u j,tt | 2 dx ds + ν 44 ∇φ ,n+1 j,h 2 . (3.26)
Last we bound the trilinear forms, i.e. b * (•, •, •). Denote
First, we add and subtract b * (u ,n j,h , u ,n+1 j,h , φ ,n+1 j,h ), and add b * (u ,n j,h ,φ ,n+1 j,h ,φ ,n+1 j,h ) = 0. We have
We add and subtract b * (u ,n j , u n+1 j , φ ,n+1 j,h ),
(3.27)
Denote
We estimate A i , where i = 1, . . . , 4, as follows. First, we bound A 1 .
(3.28)
We bound A 2 .
(3.30) Now we bound A 3 .
Last, we bound A 4 .
(3.34)
(3.38)
Combine terms,
(3.39)
By the CFL condition, we have
for some constant C 0 > 0. Recall (3.39), multiply by 2 t and organize terms:
(3.40)
Take the sum of (3.40) from n = 0 to n = N -1, we have
By Lemma 2.5, we have
By Proposition 3.2, we conclude that
By Proposition 3.4, we have
where
Apply interpolation inequalities in (2.7),
Recall that e ,n j = η ,n j -φ ,n j,h . Using the triangle inequality, we have
We complete the proof using the previous bounds for the η j terms.
Combining Theorem 3.5 with the result of Shen [12], Theorem 4.1, p. 395, and applying the triangle inequality,
We have the following corollaries. Corollary 3.6 Assume the regular solutions, under the CFL condition in (3.4), we have the following optimal estimates:
The error between the average of true solution and the average of penalized finite element approximations is
By the Cauchy Schwarz inequality,
By Corollary 3.6,
Thus,
Hence, we have
We consider the NSE with random body forces and initial conditions. We find random functions u : × D × [0, T ] → R d , and
We choose a set of random samples for the random body force f j ≡ f (ω j , •, •), initial condition u 0 j ≡ u 0 (ω j , •, •) for j = 1, . . . , J . Note that the corresponding solutions u(ω j , •, •) are independent, identically distributed (i.i.d).
The penalty-based ensemble Monte Carlo is defined as follows. Denote u ,n j,h = u h (ω j , x, t n ) and p ,n j,h = p (ω j , x, t n ). For the j th ensemble member and for 0 ≤ n ≤ N -1, find (u ,n+1 j,h , p ,n+1 j,h ) ∈ (X h , Q h ) satisfying:
for all
Theorem 3.1 together with the property of expectation leads to the following stability analysis for the finite element solution u ,n j,h .
Theorem 4.1 Suppose the following timestep condition holds:
Then for any N ≥ 1:
The fully discrete penalty-based ensemble Monte Carlo approximation is defined to be
We estimate E[u (t n )] -n h in averaged norms. We write
where
is the statistical error controls by the ensemble size.
then there are positive constant C and C 0 independent of h and t such that
Proof The conclusion follows Theorem 3.5 after applying the expectation on (3.17 Then for any N ≥ 1:
Proof Herein, we present the estimate
The last equality is due to the fact u ,n j,h for j = 1, . . . , J are i.i.d.. When i = j, the expectation of
The other terms involving E[ N S 2 ], E[ ∇ N S ] and E[ n+1 S -n S ] can be treated similarly. The statistical error from sampling n S is O( 1 √ J
). Combining Theorem 3.5 with the result of Shen [12], Theorem 4.1, p. 395, and using the triangle inequality, we will have the following corollary.
We present the results of three numerical tests to illustrate our theory. In the first test, we calculate the convergence rates using exact solutions with an ensemble size of two. Then, we construct a chaotic Lagrangian flow on a cylinder with perturbed body forces. In the third test, we extend the penalty-based ensemble algorithm to include the Coriolis force for a larger ensemble size for the benchmark test problem of flow past a cylinder. In these tests, we calculate various flow statistics to evaluate the flow dynamics:
viscous dissipation rate := ν ∇u 2 , numerical dissipation rate from backward Euler (BE)
numerical dissipation rate from penalizing incompressibility :=
We use a second-order polynomial to approximate the velocity field in the following tests. The unstructured mesh is generated by GMSH [31].
We verify the convergence rates for the method in (3.1) with a test from [32], described as follows. In D = (0, 1) 2 , the exact solution is given by u(x, y, t) = (exp(t) cos(y), exp(t) sin(x)) , p(x, y, t) = (xy)(1 + t).
We calculate the body force f by substituting u and p in the NSE. We impose the Dirichlet boundary conditions, where u h = u true on the boundary. We perturb the initial conditions as follows:
2 ) 1 • 27 0.00169 1.91 0.00639 1.91 ( 3 2 ) 2 • 27 0.00076 1.95 0.0029 1.95 ( 3 2 ) 3 • 27 0.00033 1.98 0.00127 1.98 ( 3 2 ) 4 • 27 0.00015 1.99 0.00057 1.99 We define a smooth bridging function, smooth_bridge(t), which provides a smooth transition near the edges of the active period:
The amplitude function, A(t), is periodic with period P and includes an active phase of length L within each cycle. It is defined as:
This construction ensures smooth transitions between active and inactive phases. We introduce a larger period, P large = 2P, to incorporate alternating positive and negative oscillations. The function is defined as:
2 . The left and right amplitudes are derived by phase-shifting the oscillations:
We have u(x, y) = 5A left/right (t)(y, -x) on ∂ D.
We set P = 5 and L = 2. The outer circle remains stationary. We choose mesh size h = 0.05, the final time T = 10, timestep t = 0.001, ν = 1/50 and Re = 1/ν. The penalty parameter = t. Flow is at rest at the beginning with exact boundary conditions. We perturbed the Dirichlet boundary conditions:
where σ 1 = 0.01, σ 2 = -0.02. For stability, we ensure t h ∇U ,n j ≤ 10 5 Re . This upper bound may not be necessary but sufficient. The solution u 0 is driven by the averaged Dirichlet boundary conditions across the ensemble members:
We define the ensemble spread as:
This definition captures the relative difference between the ensemble members normalized by their average. Figure 2a shows that the ensemble spread changes periodically, with the peak of the spread approximately at 0.6. We calculate the standard deviations considering u 0 and u ave . Figure 2b shows that the standard deviations for u 0 and u ave are similar. It indicates that the velocity is not chaotic. In Fig. 3a, we plot the numerical dissipation rates caused by penalizing the incompressibility condition and the BE time discretization. We compare them with the viscous dissipation rate. The numerical dissipation rate is much smaller than the viscous dissipation rate. In Fig. 3b, the numerical dissipation rates have similar magnitudes and vary over time.
We observe changes in kinetic energy, velocity divergence, angular momentum, and enstrophy as we activate and deactivate the spinning of the left and right circles over time. The flow statistics of u 0 , u 1 , u 2 , and u ave are closely aligned in Fig. 4 and indicate the velocity field is not chaotic, even though the trajectories of fluid particles exhibit chaotic behavior.
Everything on Earth is rotating even without our noticing. The rotation changes the airflow and affects the climate, as discussed in Lee, Ryi, and Lim [37]. The NSE with the Coriolis force is defined as follows:
Here Q is a skew-symmetric matrix with a matrix norm equal to one, and ω is the Coriolis coefficient. We extend the penalty-based ensemble method to the NSE with the Coriolis force. We evaluate this method using the benchmark 2D test of flow past a cylinder, as described in [38]. The inlet flow velocity is u(x, y, t) = 6y(0.41y) 0.41 2 , 0 .
We apply no-slip boundary conditions at the walls and on the obstacle. We generate second-order quadrilateral elements. We choose J = 10, T = 10, t = 0.002, ν = 0.001, and = t. The flow is at rest at t = 0. We perturb the inlet flow velocity Fig. 4 Flow statistics for u 0 , u 1 , u 2 and u 0 Fig. 5 The normalized standard deviation of the ensembles for different Coriolis coefficients for ensemble members: u j (x, y, t) = (1 + σ j sin(2π y))u, where j = 1, . . . , 10.
Here, σ j is randomly sampled from -0.1 to 0.1. We first set ω = 10. For stability, we ensure t h ∇U ,n j ≤ 10 5 Re . Figure 5a shows the spaghetti plot of the relative error of each single ensemble member to the mean flow. The normalized standard deviation for ω = 10 is around 0.15 after t = 2, as shown in Fig. 5b. We calculate the angular momentum, enstrophy, kinetic energy, and velocity divergence for all ensemble members and the mean flow, as shown in Fig. 6.
We set the Coriolis coefficient ω = 1, 10, 100, and 1000 to study the effect of the Coriolis force. We calculate the normalized standard deviation for different values of the Coriolis coefficient, as shown in Fig. 5b. For smaller ω values (ω = 1, 10, and 100), the standard deviations are similar, around 0.15. Increasing ω to 1000 significantly amplifies the rotational force, leading to a smaller standard deviation. Which suggests that the flow exhibits characteristics of rigid body rotation. We also observe much Fig. 6 Flow statistics for all ensemble members and the mean flow at ω = 10 larger magnitudes of angular momentum, enstrophy, kinetic energy, and divergence of the velocity for the ensemble mean when ω = 1000, as shown in Fig. 7.
Turbulent flows exhibit chaotic behavior, which inherently limits the predictability of numerical models to a finite horizon. The predictability relies on the accuracy of the initial conditions. Small imperfections in the initial conditions can lead to losing predictive skill. While ensemble methods effectively address this issue, they can be computationally costly. This report introduces a penalty-based ensemble method that reduces the computational cost of ensembles while preserving accuracy. The method uses a shared coefficient matrix for all ensemble members. It relaxes the incompressibility condition, uncoupling the flow velocity and pressure, reducing model complexity, and allowing for a larger ensemble size. We presented the stability and error estimates of the penalty-based ensemble method in (3.1). We extended the method to the NSE with random body forces and initial conditions with Monte Carlo sampling in Section 4. In Section 5.1, we verified the convergence rates with numerical experiments. In addition, we conducted a numerical experiment on chaotic advection, where the trajectories of the flow particles are chaotic, in Section 5.2. Furthermore, we performed a benchmark test for flow past a cylinder with the Coriolis force using large ensemble sizes in Section 5.3.
Open problems include extending penalty-based ensemble methods to turbulence models [39,40], adapting penalty parameters for ensemble methods, and exploring alternative techniques for effectively perturbing initial conditions.
-∇ u h n 2 .