High order efficient algorithm for computation of MHD flow ensembles

Numerical simulations of realistic flows are significantly affected by input data, e.g., initial conditions, boundary conditions, forcing functions, viscosities, etc, which involve uncertainties. As a result, uncertainty quantification (UQ) plays an important role in the validation of simulation methodologies and helps in developing rigorous methods to characterize the effect of the uncertainties on the final quantities of interest. A popular approach for dealing with uncertainties in the data is the computation of an ensemble average of several realizations. In many fluid dynamics applications e.g. ensemble Kalman filter approach, weather forecasting, and sensitivity analyses of solutions [14, 39, 42, 43, 44, 53] require multiple numerical simulations of a flow subject to $J$ different input conditions (realizations), which are then used to compute means and sensitivities.

Recently, the study of MHD flows has become important due to applications in e.g. engineering, physical science, geophysics and astrophysics [6, 9, 16, 19, 28, 54], liquid metal cooling of nuclear reactors [5, 24, 57], process metallurgy [15, 55], and MHD propulsion[41, 45]. For the time dependent, viscous and incompressible magnetohydrodynamic (MHD) flow simulations, this leads to solving the following $J$ separate nonlinearly coupled systems of PDEs [8, 15, 37, 47]:

Here, $u_{j}$, $B_{j}$, $p_{j}$, and $\lambda_{j}$ denote the velocity, magnetic field, pressure, and artificial magnetic pressure solutions, respectively, of the $j$-th member of the ensemble with slightly different initial conditions $u_{j}^{0}$ and $B_{j}^{0}$, and forcing functions $f_{j}$ and $\nabla\times g_{j}$ for all $j=1,2,\cdots,J$. The $\Omega\subset\mathbb{R}^{d}(d=2\hskip 2.84526pt\text{or}\hskip 2.84526pt3)$ is the convex domain, $\nu$ is the kinematic viscosity, $\nu_{m}$ is the magnetic diffusivity, $s$ is the coupling number, and $T$ is the simulation time. The artificial magnetic pressure $\lambda_{j}$ are Lagrange multipliers introduced in the induction equations to enforce divergence free constraints on the discrete induction equations but in continuous case $\lambda_{j}=0$. All the variables above are dimensionless. The magnetic diffusivity $\nu_{m}$ is defined by $\nu_{m}:=Re_{m}^{-1}=1/(\mu_{0}\sigma)$, where $\mu_{0}$ is the magnetic permeability of free space and $\sigma$ is the electric conductivity of the fluid. For the sake of simplicity of our analysis, we consider homogeneous Dirichlet boundary conditions for both velocity and magnetic fields. For periodic boundary conditions or inhomogeneous Dirichlet boundary conditions, our analyses and results will still work after minor modifications.

To obtain an accurate, even classical Navier Stokes (NSE) simulation for a single member of the ensemble, the required number of degrees of freedom (dofs) is very high, which is known from Kolmogorov’s 1941 results [38]. Thus, for a single member of MHD ensemble simulation, where velocity and magnetic fields are nonlinearly coupled, is computationally very expensive with respect to time and memory. As a result, the computational cost of the above coupled system (1)-(6) will be approximately equal to $J\times$(cost of one MHD simulation) and will generally be computationally infeasible. Our objective in this paper is to construct and study an efficient and accurate algorithm for solving the above ensemble systems. It has been shown in recent works [1, 26, 47, 59] that by using Elsässer variables formulation, efficient stable decoupled MHD simulation algorithms can be created. That is, at each time step, instead of solving a fully coupled linear system, two separate Oseen-type problems need to be solved.

Defining $v_{j}:=u_{j}+\sqrt{s}B_{j}$, $w_{j}:=u_{j}-\sqrt{s}B_{j}$, $f_{1,j}:=f_{j}+\sqrt{s}\nabla\times g_{j}$, $f_{2,j}:=f_{j}-\sqrt{s}\nabla\times g_{j}$, $q_{j}:=p_{j}+\sqrt{s}\lambda_{j}$ and $r_{j}:=p_{j}-\sqrt{s}\lambda_{j}$ produces the Elsässer variable formulation of the ensemble systems:

together with the initial and boundary conditions.

To reduce the ensemble simulation cost, a breakthrough idea was presented in [33] to find a set of $J$ solutions of the NSEs for different initial conditions and body forces. The fundamental idea is that, at each time step, each of the $J$ systems shares a common coefficient matrix, but the right-hand vectors are different. Thus, the global system matrix needs to be assembled and the preconditioner needs to be built only once per time step, and can reuse for all $J$ systems. Also, the algorithm can save storage requirements and take advantage of block linear solvers. This breakthrough idea has been implemented in heat conduction [18], Navier-Stokes simulations [30, 31, 34, 50], magnetohydrodynamics [35, 47], parameterized flow problems [23], and turbulence modeling [32]. In our earlier works in [46, 47], we adopted the idea and considered the first-order scheme with a stabilization term to compute MHD flow ensemble average subject to different initial conditions and forcing functions. In which the optimal convergence in 2D obtained under a mild criteria but in 3D the convergence was proven to be suboptimal. The objective of this paper, is to improve the temporal accuracy of the MHD ensemble scheme.

It has been shown that in Elsässer formulation simulations, the second-order approximation of certain viscous terms causes instability unless there is an unusual data restriction $1/2<\nu/\nu_{m}<2$ [40]. To overcome this issue, a practically second order $\theta-$scheme is proposed in [26], where a convex combination of the first and second order approximations of such a viscous term is taken. We extend this idea together with the efficient ensemble algorithm described above to propose a practically second order timestepping scheme for MHD flow ensemble simulation.

We consider a uniform timestep size $\Delta t$ and let $t_{n}=n\Delta t$ for $n=0,1,\cdots$., for simplicity, we suppress the spatial discretization momentarily. Then computing the $J$ solutions independently, takes the following form:
Step 1: For $j$=1,…,$J$,

Step 2: For $j$=1,…,$J$,

Where $v_{j}^{n},w_{j}^{n},q_{j}^{n}$, and $r_{j}^{n}$ denote approximations of $v_{j}(\cdot,t^{n}),w_{j}(\cdot,t^{n}),q_{j}(\cdot,t^{n})$, and $r_{j}(\cdot,t^{n})$, respectively. The ensemble mean and fluctuation about the mean are denoted by $<u>$, and $u_{j}^{{}^{\prime}}$ respectively, and they are defined as follows:

The key features to the efficiency of the above algorithm: (1) The MHD system is decoupled in a stable way into two identical Oseen problems and can be solved simultaneously if the computational resources are available. (2) At each time step, the coefficient matrices of (10)-(11) and (12)-(13) are independent of $j$. Thus, for each sub-problem, all the $J$ members in the ensemble share the same coefficient matrix. That is, at each time step, instead of solving $J$ individual systems, we only need to solve a single linear system with $J$ different right-hand-side constant vectors. (3) It provides second-order temporal accuracy when $|\nu-\nu_{m}|\ll 1$. (4) No data restriction is needed on $\nu$, and $\nu_{m}$ to avoid instability.

We give rigorous proofs that the decoupled scheme is conditionally stable and the ensemble average of the $J$ different computed solutions converges to the ensemble average of the $J$ different true solutions, as the timestep size and the spatial mesh width tend to zero. To the best of our knowledge, this second order timestepping scheme to the MHD flow ensemble averages is new.

This paper is organized as follows. Section 2 presents notation and mathematical preliminaries which are necessary for a smooth presentation and analysis to follow. In section 3, we present and analyze a fully discrete and decoupled algorithm corresponding to (10)-(13), and prove it’s stability and convergent theorems. Numerical tests are presented in section 4, and finally, conclusions and future directions are given in section 5.

Let $\Omega\subset\mathbb{R}^{d}\ (d=2,3)$ be a convex polygonal or polyhedral domain in $\mathbb{R}^{d}(d=2,3)$ with boundary $\partial\Omega$. The usual $L^{2}(\Omega)$ norm and inner product are denoted by $\|.\|$ and $(.,.)$, respectively. Similarly, the $L^{p}(\Omega)$ norms and the Sobolev $W_{p}^{k}(\Omega)$ norms are $\|.\|_{L^{p}}$ and $\|.\|_{W_{p}^{k}}$, respectively for $k\in\mathbb{N},\hskip 2.84526pt1\leq p\leq\infty$. Sobolev space $W_{2}^{k}(\Omega)$ is represented by $H^{k}(\Omega)$ with norm $\|.\|_{k}$. The natural function spaces for our problem are

Recall the Poincare inequality holds in $X$: There exists $C$ depending only on $\Omega$ satisfying for all $\phi\in X$,

The divergence free velocity space is given by

We define the trilinear form $b:X\times X\times X\rightarrow\mathbb{R}$ by

and recall from [21] that $b(u,v,v)=0$ if $u\in V$, and

The conforming finite element spaces are denoted by $X_{h}\subset X$ and $Q_{h}\subset Q$, and we assume a regular triangulation $\tau_{h}(\Omega)$, where $h$ is the maximum triangle diameter. We assume that $(X_{h},Q_{h})$ satisfies the usual discrete inf-sup condition

where $\beta$ is independent of $h$.

The space of discretely divergence free functions is defined as

For simplicity of our analysis, we will use Scott-Vogelius (SV) [56] finite element pair $(X_{h},Q_{h})=((P_{k})^{d},P_{k-1}^{disc})$, which satisfies the inf-sup condition when the mesh is created as a barycenter refinement of a regular mesh, and the polynomial degree $k\geq d$ [4, 61]. Our analysis can be extended without difficulty to any inf-sup stable element choice, however, there will be additional terms that appear in the convergence analysis if non-divergence-free elements are chosen.

We have the following approximation properties in $(X_{h},Q_{h})$: [13]

where $|\cdot|_{r}$ denotes the $H^{r}$ seminorm.

We will assume the mesh is sufficiently regular for the inverse inequality to hold, and with this and the LBB assumption, we have approximation properties

where $P_{L^{2}}^{V_{h}}(u)$ is the $L^{2}$ projection of $u$ into $V_{h}$.

The following lemma for the discrete Gronwall inequality was given in [27].

Now we present and analyze an efficient, fully discrete, decoupled, and practically second-order timestepping scheme for computing MHD flow ensemble. Like other BDF2 schemes, two initial conditions should be known; if the first initial conditions are known, using the linearized backward Euler scheme in [47] without the ensemble eddy viscosity terms on the first timestep, we can get the second initial condition without affecting the stability or accuracy. The scheme is defined below.