An efficient, partitioned ensemble algorithm for simulating ensembles of evolutionary MHD flows at low magnetic Reynolds number

Introduction

Magnetohydrodynamics (MHD) studies the dynamics of electrically conducting fluids in the presence of a magnetic field. It has many applications in astrophysics, planetary science, plasma physics and metallurgical industries, such as MHD turbulence in accretion disks , geodynamo simulations , plasma containment in fusion reactors and magnetic damping of jets and vortices . In a typical laboratory or industrial process, liquid-metal MHD usually has a modest conductivity ($\sim 10^{6}\,\Omega^{-1}m^{-1}$) and low velocity ($\sim 1m/s$), which makes the induced current densities rather modest. When this modest current density is spread over a small area ($\sim 0.1m$ in a laboratory), the induced magnetic field is usually found to be negligible by comparison with the imposed magnetic field, . Such flows, i.e. MHD flows that occur at low magnetic Reynolds number, can be modeled by the following reduced MHD system, .

Let $\Omega$ be a bounded Lipschitz domain in $R^{d}$ $(d=3)$. The governing equations of the reduced MHD system are: Given known body force $f(x,t)$ and imposed static magnetic field $B(x)$, find the fluid velocity $u(x,t)$, the pressure $p(x,t)$ and the electric potential $\phi(x,t)$ such that

Nonlinear dynamical systems such as the MHD system are sensitive to small changes in initial conditions, boundary conditions, body forces and many other input parameters. It is important to understand and quantify the limits of predictability of the system, and to develop computational approaches to reduce simulation time and computational cost while preserving a certain degree of accuracy. Most approaches to represent the uncertainties are ensemble based. Specifically, an ensemble of samples are generated to represent possible events, and then individual simulations are run for each sample. These computations are usually very expensive, and even prohibitive, especially if the size of the ensemble is large. Recently a new ensemble algorithm was proposed for fast calculation of an ensemble of the Navier-Stokes equations , which constructs linear systems with the same coefficient matrix for all realizations at each time step and thus allows the use of the either direct methods such as the LU factorization or iterative methods such as block CG , block GMRS for fast solving the linear systems. In this report, we extend the ensemble algorithm studied in to the reduced MHD system.

Herein we consider computing the reduced MHD system $J$ times with different initial conditions and/or body forces. The solution $(u_{j},p_{j},\phi_{j})$ of $j$-th realization, which corresponds to the initial condition $u^{0}_{j}(x)$ and body force $f_{j}(x,t)$, satisfies, for $j=1,2,...,J$,

Two aspects need to be considered to construct an efficient ensemble algorithm to solve the above coupled nonlinear system. The first is to use a partitioned method to uncouple the problem into two separate subproblems. This reduces solving a large linear system to solving two much smaller linear systems, which reduces the computational time and memory storage required. Furthermore, uncoupling the system also makes possible the use of highly optimized legacy code for each sub-physics problem, which reduces the main computational complexity. The other aspect is to design an ensemble algorithm for the reduced MHD system such that all ensemble members share one coefficient matrix at each time step.

To start, we first define the ensemble mean of the velocity $u_{j}^{n}$ and the electric potential $\phi_{j}^{n}$ respectively

where $u_{j}^{n}=u_{j}(x,t_{n})$, $\phi_{j}^{n}=\phi_{j}(x,t_{n})$ and $t_{n}=n\Delta t$ ($n=0,1,2,...$).

We then propose a first order, partitioned, ensemble algorithm given by

Sub-problem 1: Given $u_{j}^{n}$ and $\phi_{j}^{n}$, find $u_{j}^{n+1}$ and $p_{j}^{n+1}$ satisfying

Sub-problem 2: Given $u_{j}^{n}$, find $\phi_{j}^{n+1}$ satisfying

In Sub-problem 1, moving all the known quantities (at time level $t_{n}$) to the right hand side, one can see all ensemble members $u_{j}$ have the same coefficient matrix. Sub-problem 2 is a linear problem for $\phi_{j}$ that results in one common constant coefficient matrix for all realizations. Sub-Problem 1 and 2 are fully uncoupled at each time step and can be run in parallel. Naturally, if the ensemble is large, it can be divided into several subgroups and then one can apply the algorithm to each subgroup.

This paper is organized into four sections. In Section $2$ we establish the notation and give a weak formulation of the reduced MHD system. In Section $3$ we prove the long-time stability of the proposed algorithm under a timestep condition. In Section $4$ we present the convergence analysis of the algorithm. Several numerical examples are presented in Section $5$ to describe the implementation of the algorithm and to demonstrate its efficiency.

The ensemble method was first proposed by Jiang and Layton in to efficiently compute ensembles of Navier-Stokes equations with low/modest Reynolds numbers. For high Reynolds number flows, two ensemble eddy viscosity regularization methods were studied in , and a time relaxation algorithm in . Higher order ensemble methods can be found in . To further reduce the computation cost, incorporating reduced order modeling techniques with the ensemble algorithm was investigated in . An ensemble algorithm for computing flows with varying model parameters were developed in . The ensemble method has also been extended for computing full MHD flows in Els$\ddot{s}$sser variables in .

Notation and preliminaries

Throughout this paper the $L^{2}(\Omega)$ norm of scalars, vectors, and tensors will be denoted by $\|\cdot\|$ with the usual $L^{2}$ inner product denoted by $(\cdot,\cdot)$. $H^{k}(\Omega)$ is the Sobolev space $W_{2}^{k}(\Omega)$, with norm $\|\cdot\|_{k}$. For functions $v(x,t)$ defined on $(0,T)$, we define the norms, for $1\leq m<\infty$,

The norm on the dual space of $X$ is defined by

A weak formulation of the reduced MHD equations is: Find $u:[0,T]\rightarrow X$, $p:[0,T]\rightarrow Q$, and $\phi:[0,T]\rightarrow S$ for a.e. $t\in(0,T]$ satisfying

We will use the discrete Gronwall inequality (Lemma 2 below) in the error analysis, see for proof.

Let $D\geq 0$ and $\kappa_{n},A_{n},B_{n},C_{n}\geq 0$ for any integer $n\geq 0$ and satisfy

Suppose that for all n, $\Delta t\kappa_{n}\leq 1,$ and set $g_{n}=(1-\Delta t\kappa_{n})^{-1}$. Then,

We denote conforming velocity, pressure, potential finite element spaces based on an edge to edge triangulation ($d=2$) or tetrahedralization ($d=3$) of $\Omega$ with maximum element diameter $h$ by

We also assume the finite element spaces ($X_{h}$, $Q_{h}$) satisfy the usual discrete inf-sup /$LBB^{h}$ condition for stability of the discrete pressure, see for more on this condition. Taylor-Hood elements, e.g., , , are one such choice used in the tests in Section $5$. We further assume the finite element spaces satisfy the approximation properties of piecewise polynomials on quasiuniform meshes

where the generic constant $C>0$ is independent of mesh size $h$. An example for which the $LBB_{h}$ stability condition and the approximation properties are satisfied is the finite elements pair ($P^{k+1}$–$P^{k}$–$P^{k+1}$), $k\geq 1$. For finite element methods see for more details.

The discretely divergence free subspace of $X_{h}$ is

We assume the mesh and finite element spaces satisfy the standard inverse inequality

that is known to hold for standard finite element spaces with locally quasi-uniform meshes . We also define the standard explicitly skew-symmetric trilinear form

The full discretization of the proposed partitioned ensemble algorithm is

Sub-problem 1: Given $u_{j,h}^{n}\in X_{h}$ and $\phi_{j,h}^{n}\in S_{h}$, find $u_{j,h}^{n+1}\in X_{h}$ and $p_{j,h}^{n+1}\in Q_{h}$ satisfying

Sub-problem 2: Given $u_{j,h}^{n}\in X_{h}$, find $\phi_{j,h}^{n+1}\in S_{h}$ satisfying

Stability of the method

In this section, we prove Algorithm (3) is long time, nonlinearly stable under a CFL like time step condition.

Consider the method with a standard spacial discretization with mesh size $h$. Suppose the following time step conditions hold

Set $v_{h}=u_{j,h}^{n+1}$ in (14) and multiply through by $\Delta t$. This gives

Set $\psi_{h}=\phi_{j,h}^{n+1}$ in (15) and multiply through by $\Delta t$. This gives

The following equality will be used in the next step.

Adding (18) and (19) and using equality (3) gives

Applying Cauchy-Schwarz and Young’s inequality on the right hand side of the equation gives

Next, we bound the trilinear terms using (12), (10) and Young’s inequality.

With this bound, combining like terms, (22) becomes,

With the time step restriction (16) assumed, we have

Summing up (25) and multiplying through by 2 gives

Error Analysis

where $v^{n}=v(t_{n})$ and $t_{n}=n\Delta t$.

To analyze the rate of convergence of the approximation, we assume that the following regularity for the exact solutions:

Let $e_{u,j}^{n}=u_{j}^{n}-u_{j,h}^{n}$ and $e_{\phi,j}^{n}=\phi_{j}^{n}-\phi_{j,h}^{n}$ denote the approximation error of the $j$-th simulation at the time instance $t_{n}$. We then have the following error estimates.

For all $j=1,\ldots,J$, if the following time step conditions hold

then, there exists a positive constant $C$ independent of the time step such that

In particular, if Taylor-Hood elements ($k=2$, $s=1$) are used, i.e., the $C^{0}$ piecewise-quadratic velocity space $X_{h}$ and the $C^{0}$ piecewise-linear pressure space $Q_{h}$, and $P_{2}$ element ($m=2$) is used for $S_{h}$, we then have the following estimate.

Assume that $\|e_{u,j}^{0}\|$, $\|\nabla e_{u,j}^{0}\|$ and , $\|\nabla e_{\phi,j}^{0}\|$ are all $O(h^{2})$ accurate or better. Then, if $(X_{h},Q_{h},S_{h})$ is chosen as the $(P_{2},P_{1},P_{2})$ elements, we have

The true solution $(u_{j},p_{j},\phi_{j})$ of the reduced MHD system (2) satisfies

where \text{Intp}(u_{j}^{n+1};v_{h})=\frac{1}{N}\big{(}\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}-u_{j,t}(t^{n+1}),v_{h}\big{)}.

where $I_{h}u_{j}^{n}\in V_{h}$ is an interpolant of $u_{j}^{n}$ in $V_{h}$, and $I_{h}\phi_{j}^{n}\in S_{h}$ is an interpolant of $\phi_{j}^{n}$ in $S_{h}.$

Setting $v_{h}=U_{j,h}^{n+1}\in V_{h}$ and $\psi_{h}=\Phi_{j,h}^{n+1}\in S_{h}$, rearranging the nonlinear terms and multiply (38) by $2$, we have

Adding (39) and (40) and using equality (3) gives

We bound the terms on the right hand side of (39) as follows.

Next we analyze the nonlinear terms in (39) one by one. For the first nonlinear term, we have

Using inequality (11) and Young’s inequality, we have the following estimates.

Because $b^{*}(\cdot,\cdot,\cdot)$ is skew-symmetric, we have

For the last nonlinear term in (43), we have

Next, we bound the last two nonlinear terms on the RHS of (39) as follows:

With the assumption $u_{j}^{n+1}\in L^{\infty}(0,T;H^{1}(\Omega))$, we have

Using the inequality (13), Young’s inequality, and $u_{j}^{n+1}\in L^{\infty}(0,T;H^{1}(\Omega))$, we get

where we set $\alpha=\frac{C_{0}N}{CM^{2}}$ and $\delta=\frac{8C_{0}^{2}N^{2}}{C^{2}M^{4}}$. By Young’s inequality, inequality (13), and the result (17) from the stability analysis, i.e., $\|u_{j,h}^{n}\|^{2}\leq C$, we also have

For the pressure term in (41), because $U_{j,h}^{n+1}\in V_{h}$, for $\forall q_{j,h}^{n+1}\in Q_{h}$ we have

Combining (39)-(61), and taking $C_{0}=\frac{1}{24},C_{1}=\frac{1}{8}$, we have

By the convergence condition (28), we have

Then, after rearranging terms, (62) reduces to

Summing (63) and multiplying both sides by $2N\Delta t$ gives

Using the interpolation inequality (6) and the result (17) from the stability analysis, i.e., $\Delta t\sum_{l=0}^{n-1}\|\nabla u_{j,h}^{l+1}\|^{2}\leq CM$, we have

Applying the interpolation inequalities (5), (6), and (7) gives

We now add the following terms to both sides of (67).

Using the triangle inequality on the error equation gives

Applying the interpolation inequalities (5), (6), and (7) and absorbing constants into a new constant $C$ yields

Numerical Experiments

In this section we present numerical experiments for Algorithm 3 demonstrating the convergence and stability theorems proven in the previous sections. For all examples we will use the finite element triplet ($P^{2}$–$P^{1}$–$P^{2}$) and the finite element software package FEniCS .

For our first test problem we verify the convergence rates proven in section 4 using a variation of the test problem used in . Take the time interval $0\leq t\leq 1$, M = 16, N = 20, $\Omega=[0,\pi]^{2}$, and the imposed magnetic field $B=(0,0,1)$. We consider the true solution $(u,p,\phi)$ given by

where $\epsilon$ is a given perturbation. For this problem we will consider two perturbations $\epsilon_{1}=10^{-3}$ and $\epsilon_{2}=-10^{-3}$. The boundary conditions are taken to be $u_{h}=u_{\epsilon}$ and $\phi_{h}=\phi_{\epsilon}$ on $\partial\Omega$. The initial conditions and source terms are chosen to correspond with the exact solution for the given perturbation. As can be seen in tables 1 2 3 and 4 we achieve the expected convergence rates.

2 Efficiency Test

For our second experiment we will consider the same setting as the first numerical experiment except we will use $11$ perturbations $\epsilon_{i}=10^{-2}-.0009*i,i=0,\ldots,10$. In order to measure the efficiency of the ensemble method we compare the CPU time measured in seconds and accuracy of Algorithm 3 versus the non-ensemble IMEX version of Algorithm 3 in terms of the averages $\bar{u}^{n}$ and $\bar{\phi}^{n}$. For both algorithms we will use the direct LU solver MUMPS . We see in tables 5 and 6 that the ensemble algorithm is able to achieve similar accuracy to the non-ensemble algorithm with significant cost savings.

3 Stability Test

In this experiment we test the time step restriction for the stability of our algorithm by using a variation on the test for liquid aluminum performed in . Let $0\leq t\leq 1$, M = 12255, N = 347, $\Omega=[0,10^{-1}]^{2}$, and the imposed magnetic field $B=(0,0,1)$. We take $f$ and the boundary conditions equal to and the initial conditions to be equal to

for which we will consider the two perturbations $\epsilon_{1}=10^{-1}$ and $\epsilon_{2}=10^{-2}$. Due to the fact that there is no external energy exchange or body forces the energy in the system should decay to over time assuming the algorithm is stable. For $h=\frac{1}{10}$ we compute the average energy $E^{n}=\frac{1}{2}\|\bar{\phi}^{n}\|^{2}+\frac{1}{2}\|\bar{u}^{n}\|^{2}$ over a number of different time steps. As we can see in figure 1 our method is unstable for $\Delta t=\frac{1}{10},\frac{1}{100}$, but becomes stable with $\Delta t=\frac{1}{1000}$ .

References