Partitioned second order method for magnetohydrodynamics in Elsässer variables

Magnetohydrodynamics (MHD) studies the dynamics of electrically conducting fluids, involving Navier-Stokes equations coupled with Maxwell equations via Lorentz force and Ohm's law. Monolithic methods, which solve fully coupled MHD systems, are computationally expensive. Partitioned methods, on the other hand, decouple the full system and solve subproblems in parallel, and thus reduce the computational cost. This paper is devoted to the design and analysis of a partitioned method for the MHD system in the Elsasser variables. The stability analysis shows that for magnetic Prandtl number of order unity, the method is unconditionally stable. We prove the error estimates and present computational tests that support the theory.

The MHD flows entail two distinct physical processes: the motion of fluid is governed by hydrodynamics equations, and the magnetic field is governed by Maxwell equations. One approach to solve the coupled problem is by monolithic methods, or implicit (fully coupled) algorithms (e.g., [36]). In these methods, the globally coupled problem is assembled at each time step and then solved iteratively. Thus although robust and stable, they are quite demanding in computational time and resources.
In contrast, partitioned methods are semi-implicit algorithms that treat subphysics/subdomain problems implicitly and the coupling terms explicitly. Hence such methods are able to solve subproblems in parallel and significantly reduce the computational complexity. Partitioned methods are widely used in oceanatmosphere models, see e.g., [6]. However, there has been much less work on timedependent MHD. To the best of our knowledge, such methods are proposed in [32], [35] and [22]. The first two papers developed unconditionally stable, first order and second order partitioned methods for full MHD based on decoupling Elsässer variables respectively, while the last paper presented such methods for reduced MHD.
In this report we propose a partitioned method for decoupling the MHD in Elsässer variables. It is a two step, second order method that adopts implicit discretization on the subproblem terms and explicit discretization on coupling terms. The stability analysis shows that the method is unconditionally stable if the magnetic Prandtl number, P r m , satisfies 1/2 < P r m < 2, and is conditionally stable otherwise. In addition, the algorithm is shown to be long-time stable, i.e., the energy is bounded uniformly in time. We also perform numerical tests to verify the theory.
To specify the problem considered, we describe the full MHD equations below (see [4,28] for more details). Given a bounded domain Ω ⊂ R d , d = 2 or 3, and time T > 0, the fluid velocity field u, the magnetic field B (rescaled to give it dimensions of a velocity), and the total pressure p (kinetic and magnetic) satisfy

1)
∂B ∂t + u · ∇B − B · ∇u − ν m ∆B = 0, ∇ · B = 0, (1.2) where ν is the kinematic viscosity, ν m is the magnetic resistivity, and f is a body force. An important dimensionless parameter in MHD is the magnetic Prandtl number P r m := ν/ν m . In practice it may vary considerably depending on the medium, for instance P r m = 7 for water, 0.7 for air, and ∼ 10 −5 in the liquid core of the Earth. In a number of laboratory simulation, however, the magnetic Prandtl number is taken to be unity, or of order unity, see e.g., [4,17,25,31] and references therein.
The magnetic field can be split in two parts, B = B 0 +B, where B 0 andB are mean and fluctuation, respectively. The Elsässer variables [11] merge the physical properties of the Navier-Stokes and Maxwell equations. By adding (1.1) and (1.2), and subtracting (1.2) from (1.1), we obtain the momentum equations in Elsässer variables: The interesting property of the Elsässer variables is that there is no self-coupling in the nonlinear term in (1.3), but only cross-coupling of z + and z − . This is the basis of the Alfvén effect, which describes a fundamental interaction process, see [20,21,8,24,33,29,14,15,34]. From the point of computational view, this property may suggest the use of partitioned methods.
which is continuous and coercive. Due to the divergence-free condition, we write the nonlinear term (u · ∇v, w) in a trilinear form, b(·, ·, ·) : Note that there exists a generic constant C = C(Ω) such that (see e.g., [28]) The variational formulation for the continuous problem (1.3) is: find (z + , z − , p) : where ν ± = (ν ± ν m )/2. Denote X h and Q h the finite element spaces for X and Q respectively, built on a conforming, edge to edge triangulation with the maximum triangle parameter denoted by a subscript "h". Likewise, the discrete divergence-free space is denoted by We assume that the finite element spaces satisfy the inverse inequality: The semi-discrete approximation of (2.3) is to find 3. The partitioned method. The method we propose and analyze herein is a combination of a two-step implicit method with the coupling terms treated by an explicit discretization. Due to the symmetry of the Elsässer variables, we shall use the same time step ∆t in both subproblems. The method is described as follows: Note that the momentum equations in z + h,n+1 and z − h,n+1 are decoupled, however, the corresponding momentum equations of fluid velocity field u and magnetic field B are not: 3.1. Long-time stability of the partitioned method. The goal of this section is to demonstrate the stability of the method (3.1). It will be shown that the stability depends on the magnetic Prandtl number. More specifically, the partitioned method is unconditionally stable for 1/2 < P r m < 2; otherwise it is conditionally stable under the Courant-Friedrichs-Lewy (CFL) condition (3.13).
Lemma 3.1. If 1/2 < P r m < 2, the solution to the partitioned method (3.1) satisfies the following energy identity for any N ≥ 1, where w ± h,n = (z ± h,n , z ± h,n−1 ) T and each P n+1 is a positive term
Theorem 3.1. Assume that 1/2 < P r m < 2, then we have the following results. 1 • . The partitioned method (3.1) is unconditionally stable, i.e., for any N ≥ 1 there holds , the solution is uniformly bounded for all time:

7)
where Proof. 1 • . Due to the Poincaré inequality and Young's inequality, the forcing terms in (3.2) are bounded by Then (3.6) follows from (3.9) and by dropping the positive term P N .
and (3.9) gives 2 to both sides of the above inequality, we obtain which is equivalent to

YONG LI AND CATALIN TRENCHEA
Applying the Poincaré inequality and the equivalence of G-norm and L 2 -norm (2.1), Thus, by setting C 1 = min which, by induction, implies Setting we complete the second part (3.7).
Lemma 3.2. If P r m ≥ 2 or P r m ≤ 1/2, the solution to the partitioned method (3.1) satisfies the following energy identity for any N ≥ 1, where eachP n+1 is a positive term Proof. If P r m ≥ 2 or P r m ≤ 1/2 there no longer holds ν + > 3ν − . Thus in this case, the term associated with ν − in (3.3) is equivalent to Applying the above equality to (3.3) gives Therefore, (3.11) follows by summing (3.12) from n = 1 to N − 1. (3.13) More precisely, for any N ≥ 1 there holds If f ∈ L ∞ (0, T ; L 2 (Ω)), then the solution is uniformly bounded for all time. In addition, if f ≡ 0, then Proof. Using the inverse inequality (2.4), we have The forcing terms are bounded due to the Poincaré inequality and Young's inequality Combining (3.12), (3.16), (3.17) and droppingP n+1 yields Note that the numerical dissipation term, in the second line above, is positive under the CFL condition (3.13). By a telescope sum, we obtain (3.14).
Next, we add ∆t(ν + −|ν − |) to both sides of (3.18) and let Then (3.12) implies The rest of the proof follows analogously as in Theorem 3.1, i.e., utilizing the Poincaré inequality and the equivalence of G-norm and L 2 -norm (2.1).
Remark 3.1. The time step condition (3.13) is reasonably large for P r m of order unity. To this end, we write the time step condition as follows,

Error analysis.
In this section, we study the convergence of the method (3.1), where spatial discretization is effected using finite element methods. Recall that our finite element spaces satisfy the discrete inf-sup conditions. To establish the optimal error estimates for the approximation, we assume that the true solutions satisfy regularity conditions The errors are denoted by e ± n = z ± n − z ± h,n . Similar to the stability analysis, the error estimate depends on the values of magnetic Prandtl number in the sense of time step restriction. Nevertheless, the convergence rate is the same with respect to the mesh size and time step in both situations. For the sake of readability, the proofs are given in the Appendix.
Corollary 4.1. Suppose that (X h , Q h ) is given by P 2 -P 1 Taylor-Hood approximation elements, i.e., piecewise quadratic finite elements for z ± h and piecewise linear finite elements for p ± h . Then there is a positive constant C 0 such that 5. Numerical tests.

2D electrically conducted traveling wave problem.
In this section we verify the rate of convergence of the method (3.1) on an electrically conducted two-dimensional traveling wave problem (see [35]). The true solutions (in Elsässer variables) are ) sin(2π(y − t))e −8π 2 νt − 1 10 (y + 1) 2 e νmt − 1 4 sin(2π(x − t)) cos(2π(y − t))e −8π 2 νt − 1 10 (x + 1) 2 e νmt , p = − 1 64 (cos(4π(x − t)) + cos(4π(y − t))) e −16π 2 νt , defined on the domain Ω = [0.5, 1.5] 2 . The kinematic viscosity and magnetic resistivity are set to ν = ν m = 2.5 × 10 −4 so that P r m = 1. The time interval is 0 ≤ t ≤ 1. We adopt piecewise quadratic finite elements for z ± h and piecewise linear finite elements for p ± h . The initial data and source terms are chosen to correspond the exact solutions. According to the convergence analysis (4.4), the errors are second order with respect to the mesh size h and time step ∆t. Therefore we take ∆t = h to easily observe the convergence. Table 1 presents and confirms the rate of convergence provided by Corollary 4.1, where · ∞ = · L ∞ (0,T ;L 2 (Ω)) and · 2 = · L 2 (0,T ;L 2 (Ω)) . Figure 1 shows the log-log plot of the error for backward-Euler-forward-Euler (BEFE) (see [32]), spectral deferred correction (SDC) (see [35]) and the algorithm (3.1). Interestingly, the rate of convergence of the SDC method is slightly less than two.  Once the problem is re-written in Elsässer variables, we adopt quadratic finite elements for z ± h and piecewise linear finite elements for p ± h . Table 2 shows the L 2 error for the algorithm (3.1) at the final time with respect to several time steps. 6. Conclusion. The evolutionary coupled MHD studies the dynamics of the electrically conducting fluids and the electromagnetic fields. When solving a fully coupled MHD, it is usually computationally expensive to use a monolithic method. In contrast, a partitioned method is an attractive approach due to the decoupling of the subproblems and thus, the ability to solve them in parallel. However, to the best of our knowledge, there has been less work dedicated to partitioned methods. Perhaps, this is because of the complex (self-and cross-) couplings of the fluid velocity field and the magnetic field, which may impose a restrictive stability condition on the partitioned methods compared with monolithic methods. The Elsässer variables merges the physical properties of the Navier-Stokes and Maxwell equations. An interesting property of the MHD equations in the Elsässer variables is that it contains only the cross-coupling between z + and z − , which may suggest the use of partitioned methods. In fact, it turns out that a partitioned method in Elsässer variables is no longer a partitioned method in the fluid velocity field and magnetic field.
In this paper, we proposed and analyzed a such method (3.1) applied on the Elsässer variables, aiming to reduce the computational complexity. We presented a complete analysis on the long-time stability and error estimate. Depending on the magnetic Prandtl number, the algorithm may or may not be unconditionally stable. Yet the convergence of the error coincides in both situations.
Many open problems remain, such as developing more stable partitioned methods for large or small magnetic Prandtl number (see [18] for a recent study), and preserving the divergence-free condition of the magnetic field on the discrete level.
7. Appendix A. Proof of Theorem 4.1. The true solution (z ± , p) at time t n+1 satisfies ∂t .
Let e ± n+1 = z ± n+1 − z ± h,n+1 denote the error. We decompose it as , and We bound the first three terms on the right-hand side of (7.2) as follows, and 2 . (7.5) By adding and subtracting Using (2.2), each term above is estimated as Using the a priori bound from the stability analysis, i.e., z ∓ h,n ≤ C, we have and similarly Next, we add and subtract ν − a (2z ∓ n − z ∓ n−1 ), ξ ± n+1 to the linear term M ± n+1 so that . These terms will be estimated as follows: and (7.14) Combine (7.3)-(7.14) with ε = (ν + − 3|ν − |)/34, then (7.2) becomes YONG LI AND CATALIN TRENCHEA We then sum up (7.15) from n = 1 to n = N − 1 and multiply by 2∆t. Applying the Gronwall inequality results in  where we used the standard interpolation error estimates. As a consequence, there exists a positive constant C 0 such that To complete the error estimates, we add both sides of (7.16) with and apply the triangle inequality for the left-hand side. Noticing that the upcoming new terms are already contained in the right hand side of (7.16), we conclude the proof.
Proof of Theorem 4.2. In this case we would add and subtract ν − a z ∓ h,n+1 , ξ ± n+1 to the linear term M ± n+1 , which becomes The last terms is bounded by where we used the inverse inequality at the last step. Treating all the other terms analogously as in the previous proof, we obtain (4.3).