## A Kernel Based High Order 'Explicit' Unconditionally Stable Constrained

## Transport Method for Ideal Magnetohydrodynamics

Andrew Christlieb , Firat Cakir 1 2 , and Yan Jiang 3

## Abstract.

The ideal Magnetohydrodynamics (MHD) equations are challenging because one needs to maintain the divergence free condition, ∇· B = 0. Many numerical methods have been developed to enforce this condition. In this work, we further our work on mesh aligned constrained transport by developing a new kernel based approach for the vector potential in 2D and 3D. The approach for solving the vector potential is based on the method of lines transpose and is A-stable, eliminating the need for diffusion limiters needed in our previous work in 3D. The work presented here is an improvement over the previous method in the context of problems with strong shocks due to the fact that we could eliminate the diffusion limiter that was needed in our previous version of constrained transport. The method is robust and has been tested on the 2D and 3D cloud shock, blast wave and field loop problems.

Key Words: Magnetohydrodynamics; Constrained transport; Kernel based scheme; High order accuracy; Plasma physics.

1 Department of Computational Mathematics, Science and Engineering, Department of Mathematics and Department of Electrical Engineering, Michigan State University, East Lansing, MI, 48824. E-mail: christli@msu.edu

2 Department of Mathematics, Michigan State University, East Lansing, MI, 48824. E-mail: cakirfir@msu.edu

3 School of Mathematical Sciences, University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China. E-mail: jiangy@ustc.edu.cn

## 1 Introduction

The ideal Magnetohydrodynamics (MHD) equations are one of the most important classical models of plasma physics explaining the macroscopic phenomena of a quasi-neutral plasma system. The ideal MHD equations are the set of transport evolution equations for the quantities of mass, momentum, and energy density as well as the magnetic field in a conducting fluid. Mathematically, the MHD equations are a system of nonlinear hyperbolic conservation laws with divergence-free magnetic field condition. Satisfying this condition in numerical simulations is the main challenge for the numerical methods. Many standard schemes fail to guarantee ∇ · B = 0. There have been overall four different dominant approaches overcoming the difficult in the literature: 8-wave formulation [41, 42]; projection methods [4, 48]; hyperbolic divergence cleaning methods [24]; and constrained transport methods [2, 7, 22, 27, 28, 31, 38, 39, 43, 44, 23, 48, 47, 21].

The first methodology used for divergence free condition is developed in [9] utilizing the classical projection method. They solve a Poisson equation to project the incorrect magnetic field to a divergence free subspace using the Hodge decomposition. However, it is difficult to extend the method to an Adaptive Mesh Refinement (AMR) approach, since this method has to solve a Poisson equation on each time step. The AMR idea is to refine the mesh based on an accuracy criterion, providing increased resolution at shock and increased computational efficiency. Zachary [51] presented a Riemann solver which removes negative pressures and densities and used the Projection method to get a divergence free magnetic field. The 8-wave scheme is the second approach, developed by Powell in [41]. He adds an extra source term to update the magnetic field which satisfies the divergence free condition. By this treatment, the ideal MHD equations, with constraint, becomes a 8 × 8 hyperbolic system with a source term. This scheme is robust and can easily be extended to an AMR framework. However, it has been shown in [48] that this method can create inaccurate jumps in a discontinuous example, such as rotated shock tube problem, since this procedure is non-conservative because of the source term. The third method introduced by Dedner [24] is the hyperbolic divergence-cleaning method. This method is a similar to the projection method. In this method, hyperbolic and parabolic corrections are combined together and solved for the magnetic field divergence error. This method is fully explicit, efficient and fast. However, this method has two tunable parameters, the speed of propagation of the error and the damped divergence error rate, which needed to be adjusted to ensure good solutions. This method gives a

damped hyperbolic equation for the divergence error and it does not exactly satisfy divergence free condition.

The main focus in this work is the constrained transport (CT) method for correcting the magnetic field to satisfy the divergence free constraint. This methodology was invented by Evans and Hawley in [27]. The original CT method is considered to be a modification of the Yee method [50] from electromagnetics for the ideal MHD equations. In the original constrained transport methodology, staggered electric and magnetic fields are used to create appropriate finite difference operators. These operators eventually lead to a globally divergence-free magnetic field.

In the literature, various modification of the constrained transport method have been presented. In particular, high resolution shock capturing schemes have been a primary focus. DeVore [25] presented an application of a flux corrected transport approach satisfying a divergence free magnetic field. There are a range of different approaches for building the electric field using Ohm's law in a constrained transport methodology including those presented by Balsara and Spicer [7], Dai and Woodard [22], and Ryu et al. [44]. Londrillo and Zanna [38, 39] constructed one of the first high order upwind scheme based on the work of Evans and Hawley. De Sterck [23] introduced a similar constrained transport scheme on unstructured triangle grids based on multidimensional upwind advection schemes. Balsara [1] described a divergence free adaptive mesh refinement (AMR) method utilizing a constrained transport approach. T´ oth [48] compared several of these schemes maintaining divergence free condition and also showed a staggered magnetic field is unnecessary as well as developing some unstaggered CT frameworks.

Unstaggered Constrained Transport schemes are getting increased attention over the last few years, since it is easy to implement and to extend to adaptive mesh refinement (AMR) methods. For example , Fey and Torrilhon [28] developed an unstaggered upwind method satisfying the divergence free constraint. Rossmanith [43] designed an ustaggered wave propagation scheme for MHD flows based on the algorithms in [36] using a constrained transport method to keep the divergence free magnetic field. Helzel et al. [31, 32] generalized the 2D unstaggered CT work to 3D MHD equations so that the method is applicable on both Cartesian and rectangular mapped grids.

In resent years, several high order methods have been developed for the ideal MHD equations using various discretization approaches. Balsara [3] designed a third order divergence-free weighted

essentially non-oscillatory (WENO) methods for MHD equations with third order Runge-Kutta time integration using a staggered magnetic field. Balsara et al. [5, 6] introduced a high accuracy ADER-WENO schemes for divergence free magnetohydrodynamics on structured meshes, again using a staggered magnetic field. Li et al. [37, 29] and Cheng et al. [15] presented high order central discontinuous Galerkin schemes satisfying the divergence free constraint globally for ideal MHD simulations on two overlapping meshes, called primal and dual meshes, by utilizing different discretization for magnetic induction equations. Kawai [35] introduced a divergence-free high order accurate finite difference scheme which has an effective shock capturing capability for the MHD equations by constructing artificial diffusion terms to capture numerical discontinuous in the magnetic field.

In our previous paper [21], we introduced a high order FD-WENO for the ideal MHD in 2D and 3D by developing a high order unstaggered constrained transport methodology for Hamilton Jacobi equations using a version of FD-WENO. In that work, WENO formulation is applied to the central derivative of the solution A instead of the flux values on grid points to approximate the one-sided partial derivative terms A -x and A + x appear in HJ equation. If explicit time stepping is used for HJ equation, then it looks like a convection reaction equation and it tends to be unstable. That's why we needed to add artificial resistivity terms for the 3D case to stabilize it and to control unphysical oscillations in the magnetic field.

In this current work, we further our work [21] on mesh aligned constrained transport by developing a new kernel-based approach for the magnetic vector potential to solve the ideal MHD equations in 2D and 3D. The approach is based on a kernel-based numerical scheme [17, 18], which is derived from the Method of Lines Transpose (MOL T ) [14, 11, 13, 12, 16, 10]. The equation is first discretized in time with an explicit method and transformed to a boundary value problem (BVP) at discrete time levels. The spatial operators are converted from local representations, i.e.. derivatives, to global representations using convolutions with kernels, i.e. Green's functions, and is similar in spirt of taking a Fast Fourier transform. Thus our method is fully implicit since the partial derivative terms are represented by global operators using convolution integrals which are communicating all the previous and future data. By this methodology, we update the predicted magnetic field obtained from the base scheme by a corrected divergence free magnetic field. The approach for solving the vector potential is derived from the method of lines transpose and is A-

stable, eliminating the need of diffusion limiters introduced in our previous work in 3D [21] since the kernel based method is a fully implicit method. The corrected magnetic field is computed by 4th-order accurate central finite difference operators that approximates the curl of the magnetic vector potential. The magnetic vector potential is made to satisfy a weakly hyperbolic system using the Weyl gauge condition. This system is solved using our kernel based scheme developed for Hamilton-Jacobi equations. Our solver is coupled with the 5th-order FD-WENO scheme of Jiang and Shu [34] as the base scheme for ideal MHD equations. Third order explicit strong-stabilitypreserving (SSP) Runge-Kutta (RK) is used for the time discretization. The work presented here is an improvement over the previous method in the context of problems with strong shocks due to the fact that, using the kernel-based approach, we could eliminate the diffusion limiter that was needed in our previous version of constrained transport. This method is robust and has been tested on the 2D and 3D cloud shock, blast wave and field loop problems.

The rest of this paper is organized as follows: in section 2, we briefly review the MHD equations and the evolution of the magnetic vector equations; in section 3, we give a brief outline of the CT algorithm; in section 4, we present our novel numerical scheme for 1D Hamilton-Jacobi equations; we introduce the multidimensional solver in section 5; the resulting 2D and 3D schemes are tested on several numerical problems in section 6.

## 2 The Ideal MHD Equations

In this section we present a brief review of the ideal MHD equations, which is a first order hyperbolic system of conservation laws. The conservative form of the ideal MHD equations can be written as

$$\partial _ { t } \begin{bmatrix} \rho \\ \rho \bar { u } \\ \mathcal { E } \\ \mathcal { B } \end{bmatrix} + \nabla \cdot \begin{bmatrix} \rho \bar { u } \\ \bar { u } \otimes \bar { u } + p _ { t o t } \mathbb { I } - B \otimes \mathcal { B } \\ \bar { u } \otimes \mathcal { B } - B \otimes \bar { u } \end{bmatrix} = 0,$$

$$\nabla \cdot { \mathbf B } = 0,$$

with the equation of state as

$$\mathcal { E } = \frac { p } { \gamma - 1 } + \frac { \rho \left \| \mathfrak { u } \right \| ^ { 2 } } { 2 } + \frac { \left \| \mathbf B \right \| ^ { 2 } } { 2 }.$$

where the total mass ρ , the momentum ρ u = ( ρu , ρu , ρu 1 2 3 ) T , the energy densities E of the plasma system, and the magnetic field vector B = ( B ,B ,B 1 2 3 ) T are all conserved variables. The velocity

u and the total pressure p tot = p + 1 2 ‖ B ‖ 2 , together with the hydrodynamic pressure, which is given by the ideal gas law as

$$p = ( \gamma - 1 ) ( \mathcal { E } - \frac { 1 } { 2 } \| \mathbf B \| ^ { 2 } - \frac { 1 } { 2 } \rho \| \mathbf u \| ^ { 2 } ),$$

are derived quantities. γ = 5 3 is the ideal gas constant and the notation / ‖ · ‖ is used for the purpose of the Euclidean vector norm. The MHD equations (2.1)-(2.2) are derived and discussed in many standard plasma textbooks (e.g., [40]).

## 2.1 Hyperbolicity of the governing equations

The system (2.1) along with the equation of state (2.3) comprise a system of hyperbolic conservation laws

$$\mathbf q _ { t } + \nabla \cdot \mathbf F ( \mathbf q ) = 0.$$

This set of equations describes the time evolution of all eight conserved variables, q = ( ρ, ρ u , E , B ). We will denote the Jacobian of the hyperbolic system as M ( q ) = ∂ F /∂ q , which is a diagonalizable matrix with real eigenvalues. The eigenvalues of the Jacobian matrix represent the wave speeds of the eight waves in MHD system. In particular, in some arbitrary direction n ( ‖ n ‖ = 1), they can be written as

$$\lambda ^ { 1, 8 } & = \mathbf u \cdot \mathbf n \mp c _ { f } & \quad \text{fast magnetosonic waves, } \quad \text{(2.6a)}$$

$$\lambda ^ { 2, 7 } & = \mathbf u \cdot \mathbf n \mp c _ { a } & \text{Alvén waves},$$

$$\lambda ^ { 3, 6 } & = \mathbf u \cdot \mathbf n \mp c _ { s } & \text{ slow magnetosonic waves, } \text{ \quad \ \ } ( 2. 6 c )$$

$$\lambda ^ { 4 } & = u \cdot n & & \text{entropy waves}, & & ( 2. 6 d )$$

$$\lambda ^ { 5 } & = u \cdot n & & \text{divergence waves}, & & ( 2. 6 e )$$

$$a \equiv \sqrt { \frac { \gamma p } { \rho } }$$

$$c _ { a } \equiv \sqrt { \frac { ( { \mathbf B } \cdot { \mathbf n } ) ^ { 2 } } { \rho } }$$

$$c _ { f } \equiv \begin{array} { c } \cdot \\ \frac { 1 } { 2 } \left [ a ^ { 2 } + \frac { \| \mathbf B \| ^ { 2 } } { \rho } + \sqrt { \left ( a ^ { 2 } + \frac { \| \mathbf B \| ^ { 2 } } { \rho } \right ) ^ { 2 } - 4 a ^ { 2 } \frac { ( \mathbf B \cdot \mathbf n ) ^ { 2 } } { \rho } } \right ] \right \} ^ { 1 / 2 } \end{array}$$

$$c _ { s } \equiv \left \{ \frac { 1 } { 2 } \left [ a ^ { 2 } + \frac { \| { \mathbf B } \| ^ { 2 } } { \rho } - \sqrt { \left ( a ^ { 2 } + \frac { \| { \mathbf B } \| ^ { 2 } } { \rho } \right ) ^ { 2 } - 4 a ^ { 2 } \frac { ( { \mathbf B } \cdot { \mathbf n } ) ^ { 2 } } { \rho } } \right ] \right \} ^ { 1 / 2 }$$

where

The eight eigenvalues are well ordered as

$$\lambda ^ { 1 } \leq \lambda ^ { 2 } \leq \lambda ^ { 3 } \leq \lambda ^ { 4 } \leq \lambda ^ { 5 } \leq \lambda ^ { 6 } \leq \lambda ^ { 7 } \leq \lambda ^ { 8 }$$

where the fast and slow magnetosonic waves are nonlinear while the rest of the waves are linearly degenerate.

## 2.2 Magnetic potential in 3D

There have been many numerical methods presented in the literature for numerically solving the MHD system, but they have faced the main challenge of satisfying the divergence free condition on the magnetic field. Here, we will derive magnetic vector potential equation from the magnetic field equation given in (2.1)-(2.2) system. This will serve as the foundation of our constrained transport framework.

Since the magnetic field is divergence free, it can always be written as the curl of a magnetic vector potential

$$B = \nabla \times \mathbf A.$$

The key step of the constrained transport scheme is to solve the magnetic potential for correcting the magnetic field. The evolution of the magnetic field in (2.1) can be written in the following form

$$B _ { t } + \nabla \times ( B \times \mathfrak { u } ) = 0,$$

using the relation

$$\nabla \cdot ( \mathbf u \otimes \mathbf B - \mathbf B \otimes \mathbf u ) = \nabla \times ( \mathbf B \times \mathbf u ).$$

Since B is divergence free, we set B = ∇× A and rewrite the evolution equation (2.10) as

$$\nabla \times \left \{ \mathbf A _ { t } + ( \nabla \times \mathbf A ) \times \mathbf u \right \} = 0.$$

This implies that there exists a scalar function ψ such that

$$\mathbf A _ { t } + ( \nabla \times \mathbf A ) \times \mathbf u = - \nabla \psi.$$

There are various choices of the gauge conditions depending on how we chose the ψ . Helzel et al. [31] showed that using the Weyl gauge, i.e., setting ψ ≡ 0, one can achieve stable solutions. This condition results in the evolution equation for the magnetic vector potential as

$$\partial _ { t } \mathbf A + ( \nabla \times \mathbf A ) \times \mathbf u = 0,$$

which can be written as:

$$\partial _ { t } \mathbf A + N _ { 1 } \mathbf A _ { x } + N _ { 2 } \mathbf A _ { y } + N _ { 3 } \mathbf A _ { z } = 0,$$

with

$$N _ { 1 } = \begin{bmatrix} 0 & - u ^ { 2 } & - u ^ { 3 } \\ 0 & u ^ { 1 } & 0 \\ 0 & 0 & u ^ { 1 } \end{bmatrix}, \quad N _ { 2 } = \begin{bmatrix} u ^ { 2 } & 0 & 0 \\ - u ^ { 1 } & 0 & - u ^ { 3 } \\ 0 & 0 & u ^ { 2 } \end{bmatrix}, \quad N _ { 3 } = \begin{bmatrix} u ^ { 3 } & 0 & 0 \\ 0 & u ^ { 3 } & 0 \\ - u ^ { 1 } & - u ^ { 2 } & 0 \end{bmatrix}.$$

The resulting system is only weakly hyperbolic since the matrix of right eigenvectors of the flux Jacobian doesn't have full rank in certain directions. We begin with the flux Jacobian matrix in some arbitrary direction n = ( n , n 1 2 , n 3 ) to show the weakly hyperbolicity:

$$M ( \mathbf n, \mathbf u ) \coloneqq n ^ { 1 } N _ { 1 } + n ^ { 2 } N _ { 2 } + n ^ { 3 } N _ { 3 } = \begin{bmatrix} n ^ { 2 } u ^ { 2 } + n ^ { 3 } u ^ { 3 } & - n ^ { 1 } u ^ { 2 } & - n ^ { 1 } u ^ { 3 } \\ - n ^ { 2 } u ^ { 1 } & n ^ { 1 } u ^ { 1 } + n ^ { 3 } u ^ { 3 } & - n ^ { 2 } u ^ { 3 } \\ - n ^ { 3 } u ^ { 1 } & - n ^ { 3 } u ^ { 2 } & n ^ { 1 } u ^ { 1 } + n ^ { 2 } u ^ { 2 } \end{bmatrix}, \quad \text{(2.13)}$$

which has real eigenvalues for all ‖ n ‖ = 1 as

$$\lambda ^ { 1 } = 0, \quad \lambda ^ { 2 } = \lambda ^ { 3 } = \mathbf n \cdot \mathbf u,$$

and the right eigenvectors matrix is

$$R = \left [ r ^ { ( 1 ) } \Big | r ^ { ( 2 ) } \Big | r ^ { ( 3 ) } \right ] = \begin{bmatrix} n ^ { 1 } & n ^ { 2 } u ^ { 3 } - n ^ { 3 } u ^ { 2 } & u ^ { 1 } ( \mathbf u \cdot \mathbf n ) - n ^ { 1 } \| \mathbf n \| ^ { 2 } \\ n ^ { 2 } & n ^ { 3 } u ^ { 1 } - n ^ { 1 } u ^ { 3 } & u ^ { 2 } ( \mathbf u \cdot \mathbf n ) - n ^ { 2 } \| \mathbf n \| ^ { 2 } \\ n ^ { 3 } & n ^ { 1 } u ^ { 2 } - n ^ { 2 } u ^ { 1 } & u ^ { 3 } ( \mathbf u \cdot \mathbf n ) - n ^ { 3 } \| \mathbf n \| ^ { 2 } \end{bmatrix}.$$

If we assume that ‖ u ‖ glyph[negationslash] = 0 and ‖ n ‖ = 1, then the determinant of the matrix of the right eigenvector R can be found as

$$d e t ( R ) = - \| \mathbf u \| ^ { 3 } \cos ( \alpha ) \sin ^ { 2 } ( \alpha )$$

where α is the angle between the vectors n and u . Hence, given any nonzero velocity vector u , there exist four degenerate directions α = 0 , π/ 2 , π, 3 π/ 2, where the eigenvectors are incomplete due to det R ( ) = 0. Therefore, the system (2.12) is only weak hyperbolic.

## 2.3 Magnetic potential in 2D

In two dimensional case (e.g., in the xy-plane), each conserved variable, q = ( ρ, ρ u , E , B ), depends only on t , x , and y independent variables. Thus, the divergence free condition can be simplified as

$$\nabla \cdot { \mathbf B } = B _ { x } ^ { 1 } + B _ { y } ^ { 2 }.$$

There is no point to define B 3 since any value of B 3 satisfies the divergence free condition in 2D. The evolution of the magnetic field involves only the third component of the magnetic vector potential, so that magnetic field components, B 1 and B 2 , can be found as

$$B ^ { 1 } = A _ { y } ^ { 3 } \text{ \ and \ } B ^ { 2 } = - A _ { x } ^ { 3 }.$$

Since we only need the third component of the magnetic vector potential to define magnetic field, the constrained transport scheme in 2D is reduced to the magnetic vector potential to a magnetic scalar potential equation as

$$A _ { t } ^ { 3 } + u ^ { 1 } A _ { x } ^ { 3 } + u ^ { 2 } A _ { y } ^ { 3 } = 0$$

This scalar advection equation is strongly hyperbolic by contrast with the 3D case.

## 3 Outline of the constrained transport methodology

The major challenge when dealing with the numerical solution of the system (2.1)-(2.3) is to satisfy the divergence-free condition (2.2). The constrained transport methodology is one of the dominant approach in the literature to overcome this challenge. In this section we will give the main idea of this method and our perspective how we are using high order CT framework. The idea of the constrained transport methodology is based on updating the conserved variables q = ( ρ, ρ u , E , B ). One of the early works on this idea of using the magnetic vector potential equations for CT solution to the MHD equations is presented by Wilson in [49], as well as Dorfi [26]. However, modern shock capturing strategies were not used in those works and hence caused strong numerical diffusion. Londrillo and Del Zanna [38] used the magnetic potential solutions in the context of shock capturing methods, along with De Sterck [23], Londrillo and Del Zanna [39], and Rossmanith [43]. Helzet et al. [32] developed a framework for unstaggered CT scheme coupling a conservative finite volume hyperbolic scheme for the MHD equations with a non conservative finite volume scheme for the magnetic vector potential equation. In our previous work [21], we developed a high order finite difference constrained transport scheme. In this paper, we present a kernel based high order unconditionally stable CT method. We correct the conserved variables, q , using the magnetic vector potential formulation with respect to relation B = ∇× A . Here, we give an outline of the general unstaggered CT framework listing all important steps. Consider a semi discrete system of

ordinary differential equations for MHD equations (2.1)

$$\mathfrak { q } _ { m h d } ^ { \prime } ( t ) = \mathcal { L } _ { 1 } ( \mathfrak { q } _ { m h d } ( t ) )$$

where q mhd ( ) t represents the grid function at time t consisting of all point-wise values of the conserved quantities in the ideal MHD system q mhd = ( ρ, ρ u , E , B ). The details about L 1 ( q mhd ( )) t were presented in [21].

And also consider the update for the magnetic potential equation ((2.17) for 2D case and (2.12) for 3D case) on the mesh,

$$\mathbf q _ { A } ^ { \prime } ( t ) = \mathcal { L } _ { 2 } ( \mathbf q _ { A } ( t ), \mathbf q _ { m h d } ( t ) ),$$

where q A denotes vector potential A in 3D or the scalar potential A 3 in 2D.

The key steps advancing the solution from its current time step t = t n (or the initial condition at t 0 ) to its new time step t n +1 are listed below:

- · step 0: Start with the given current time step q n mhd = ( ρ , ρ n u n , E n , B n T ) and q n A .
- · step 1: Obtain q ∗ mhd and q n +1 A separately, where

$$\mathbf q _ { m h d } ^ { * } = \left ( \rho ^ { n + 1 }, \rho ^ { n + 1 } \mathbf u ^ { n + 1 }, \mathcal { E } ^ { * }, \mathbf B ^ { * } \right ).$$

Here, E ∗ and B ∗ are given with a ∗ superscript instead of n +1 to indicate that the predicted B and E will be corrected by a predictor-corrector Constrained Transpose method before the end of the time step.

- · step 2: Replace B ∗ to B n +1 by a discrete curl of q n +1 A .

$$\mathbf B ^ { n + 1 } = \nabla \times \mathbf q _ { A } ^ { n + 1 }.$$

- · step 3: Set the corrected total energy density value E n +1 based on one of the following options:

Option 1: Keep the total energy conserved

$$\mathcal { E } ^ { n + 1 } = \mathcal { E } ^ { * }.$$

Option 2: Keep the pressure the same after updating the magnetic field

$$\mathcal { E } ^ { n + 1 } = \mathcal { E } ^ { * } + \frac { 1 } { 2 } ( \left \| \mathbf B ^ { n + 1 } \right \| ^ { 2 } - \left \| \mathbf B ^ { * } \right \| ^ { 2 } )$$

(Second option sometimes helps to prevent negative pressure).

In Section 5, we show that our approach leads to a divergence free solution. We now describe how we constraint the update for q n +1 A using our kernel based approach.

## 4 Hamilton-Jacobi equations

In this section we introduce the main ideas in our kernel-based method [17, 18], which is derived from the Method of Lines Transpose (MOL T ) [14, 11, 13, 12, 16, 10]. The simplest way to describe MOLT is: we start by discretizing the problem in time; we then use a global approximation for the inherently local term (the derivative). By doing this, we are able to make an explicit approximation unconditionally stable. As for the global approximation, think FFT, but with a different convolution kernel. The form of the approximation is what facilitates the stability of the method. The method is an O(N) because the convolution with the kernels can be evaluated using a three term recreation.

## 4.1 1D Hamilton-Jacobi equations

Consider a 1D Hamilton-Jacobi equation

$$A _ { t } + H ( A _ { x } ) = 0, \quad A ( x, 0 ) = A ^ { 0 } ( x ), \quad x \in [ a, b ],$$

where H is the Hamiltonian flux. For the time discretization purpose to evolve the solution from time t n to t n +1 , we use the classical explicit SSP RK schemes [30]. In this work, we propose to use the following SSP RK schemes such as the first order forward Euler scheme

$$A ^ { n + 1 } = A ^ { n } - \Delta t \hat { H } ( A _ { x } ^ { n, - }, A _ { x } ^ { n, + } ) ;$$

the second order SSP RK scheme

$$A ^ { ( 1 ) } & = A ^ { n } - \Delta t \hat { H } ( A _ { x } ^ { n, - }, A _ { x } ^ { n, + } ), \\ A ^ { n + 1 } & = \frac { 1 } { 2 } A ^ { n } + \frac { 1 } { 2 } \left ( A ^ { ( 1 ) } - \Delta t \hat { H } ( A _ { x } ^ { ( 1 ), - }, A _ { x } ^ { ( 1 ), + } ) \right ) ;$$

and the third order SSP RK scheme

$$A ^ { ( 1 ) } & = u ^ { n } - \Delta t \hat { H } ( A _ { x } ^ { n, - }, A _ { x } ^ { n, + } ), \\ A ^ { ( 2 ) } & = \frac { 3 } { 4 } A ^ { n } + \frac { 1 } { 4 } \left ( A ^ { ( 1 ) } - \Delta t \hat { H } ( A _ { x } ^ { ( 1 ), - }, A _ { x } ^ { ( 1 ), + } ) \right ), \\ A ^ { n + 1 } & = \frac { 1 } { 3 } A ^ { n } + \frac { 2 } { 3 } \left ( A ^ { ( 2 ) } - \Delta t \hat { H } ( A _ { x } ^ { ( 2 ), - }, A _ { x } ^ { ( 2 ), + } ) \right ).$$

where ∆ t denotes the time step and ˆ H is the the numerical Hamiltonian, e.g., the local LaxFriedrichs Hamiltonian flux [33]:

$$\hat { H } ( \phi ^ { - }, \phi ^ { + } ) = H ( \frac { \phi ^ { - } + \phi ^ { + } } { 2 } ) - c ( \phi ^ { - }, \phi ^ { + } ) \frac { ( \phi ^ { + } - \phi ^ { - } ) } { 2 },$$

with c φ ( -, φ + ) = max φ ∈ [min( φ -,φ + ) max( , φ -,φ + )] | H φ ′ ( ) . | Here, A -x and A + x are one-side derivatives with left -biased and right -biased methods, respectively, to approximate A x . Finding these derivatives has the dominant role in our work and we will present the details of the construction of them in the following subsection.

## 4.2 Approximation of the first order derivative ∂ x

Approximating the partial spatial derivative terms with kernel based scheme is the major part of this work. In this section we briefly review the construction of the ∂ x derivative approximation using kernel based formulation established in [17]. Let's define operators L L and L R , and their inverse operators

$$\mathcal { L } _ { L } = \mathcal { I } + \frac { 1 } { \alpha } \partial _ { x } \ \Rightarrow \ \mathcal { L } _ { L } ^ { - 1 } [ v, \alpha ] ( x ) = \alpha \int _ { a _ { \iota } } ^ { x } e ^ { - \alpha ( x - y ) } v ( y ) d y + A _ { L } e ^ { - \alpha ( x - a ) },$$

$$\mathcal { L } _ { R } = \mathcal { I } - \frac { 1 } { \alpha } \partial _ { x } \ \Rightarrow \ \mathcal { L } _ { R } ^ { - 1 } [ v, \alpha ] ( x ) = \alpha \stackrel { \mu \mu } { \int _ { x } ^ { b } e ^ { - \alpha ( y - x ) } v ( y ) d y } + B _ { R } e ^ { - \alpha ( b - x ) },$$

where I is the identity operator and α is a positive constant. We use I L and I R to denote the convolution integral as

$$I ^ { L } = \alpha \int _ { a } ^ { x } e ^ { - \alpha ( x - y ) } v ( y ) d y, \quad I ^ { R } = \alpha \int _ { x } ^ { b } e ^ { - \alpha ( y - x ) } v ( y ) d y.$$

We can see that I L depends on the function values of v from left end point a to x , as well as I R is for right end point b to x . Also note that A L and B R are determined by the boundary conditions. For instance, if we assume periodic boundary conditions, i.e.,

$$\mathcal { L } _ { L } ^ { - 1 } [ v, \alpha ] ( a ) = \mathcal { L } _ { L } ^ { - 1 } [ v, \alpha ] ( b ), \quad \text{and}, \quad \mathcal { L } _ { R } ^ { - 1 } [ v, \alpha ] ( a ) = \mathcal { L } _ { R } ^ { - 1 } [ v, \alpha ] ( b ).$$

then we obtain

$$A _ { L } = \frac { I ^ { L } [ v, \alpha ] ( b ) } { 1 - \mu }, \quad \text{and} \quad B _ { R } = \frac { I ^ { R } [ v, \alpha ] ( a ) } { 1 - \mu },$$

where µ = e -α b ( -a ) .

Now, let's define new operators D L and D R ,

$$\mathcal { D } _ { L } = \mathcal { I } - \mathcal { L } _ { L } ^ { - 1 }, \text{ \ and \ } \mathcal { D } _ { R } = \mathcal { I } - \mathcal { L } _ { R } ^ { - 1 }.$$

Then 1 α ∂ x can be represented using the infinite series of these operators as

$$\frac { 1 } { \alpha } \partial _ { x } = \mathcal { L } _ { L } - \mathcal { I } = \mathcal { L } _ { L } ( \mathcal { I } - \mathcal { L } _ { L } ^ { - 1 } ) = \mathcal { D } _ { L } ( \mathcal { I } - \mathcal { D } _ { L } ) ^ { - 1 } = \sum _ { p = 1 } ^ { \infty } \mathcal { D } _ { L } ^ { p }$$

$$\frac { 1 } { \alpha } \partial _ { x } = \mathcal { I } - \mathcal { L } _ { R } = \mathcal { L } _ { R } ( \mathcal { L } _ { R } ^ { - 1 } - \mathcal { I } ) = - \mathcal { D } _ { R } ( \mathcal { I } - \mathcal { D } _ { R } ) ^ { - 1 } = - \sum _ { p = 1 } ^ { \infty } \mathcal { D } _ { R } ^ { p }$$

Hence, if A is a periodic function, then we can obtain the approximations of the first derivatives A ± x using partial sums as

$$\lambda \cdot \\ A _ { x } ^ { - } & \approx \begin{cases} \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { L } ^ { p } [ A, \alpha ] ( x ), & k = 1, 2, \\ \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { L } ^ { p } [ A, \alpha ] ( x ) - \alpha \mathcal { D } _ { 0 } * \mathcal { D } _ { L } ^ { 2 } [ A, \alpha ] ( x ), & k = 3, \end{cases}$$

and

$$\lambda _ { x } ^ { + } & \approx \begin{cases} \ - \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { R } ^ { p } [ A, \alpha ] ( x ), & k = 1, 2, \\ - \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { R } ^ { p } [ A, \alpha ] ( x ) + \alpha \mathcal { D } _ { 0 } * \mathcal { D } _ { R } ^ { 2 } [ A, \alpha ] ( x ), & k = 3. \\ \cdot \quad \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \real \null \end{cases}$$

For k = 3 case, there is an additional term D 0 which is defined as

$$\mathcal { D } _ { 0 } [ v, a ] = v ( a ) - \frac { \alpha } { 2 } \int _ { a } ^ { b } e ^ { - \alpha | x - y | } v ( y ) d y - A _ { 0 } e ^ { - \alpha ( x - a ) } - B _ { 0 } e ^ { - \alpha ( b - x ) }.$$

where A 0 and B 0 are determined by boundary condition. For example, if D 0 is a periodic function such that

$$\mathcal { D } _ { 0 } [ v, \alpha ] ( a ) = \mathcal { D } _ { 0 } [ v, \alpha ] ( b ), \text{ \ and \ } \partial _ { x } \mathcal { D } _ { 0 } [ v, \alpha ] ( a ) = \partial _ { x } \mathcal { D } _ { 0 } [ v, \alpha ] ( b ),$$

then we get

$$A _ { 0 } = \frac { I ^ { 0 } [ v, \alpha ] ( b ) } { 1 - \mu }, \quad B _ { 0 } = \frac { I _ { 0 } [ v, \alpha ] ( a ) } { 1 - \mu },$$

with I 0 [ v, α ]( x ) = α 2 ∫ b a e - | - | α x y v y dy ( ) .

Here, we take α = β/ c ( ∆ ), where c is the maximum wave speed and t β is a constant independent of the time step ∆ . Therefore, the accuracy for the approximation (4.11) to the t ∂ A x is O (∆ t k ). If β is chosen appropriately, then the scheme is A-stable and hence allows for large time step evolution. Additionally, the linear stability of the method which argued in the following theorem has been proven in [17].

Theorem 4.1. For the linear equation φ t + cφ x = 0 , (i.e. the Hamiltonian is linear) with periodic boundary conditions, we consider the k th order SSP RK method as well as the k th partial sum in (4.11) , with α = β/ ( | c | ∆ ) t . Then there exists a constant β k,max &gt; 0 for k = 1 , 2 , 3 , such that the scheme is A-stable provided 0 &lt; β ≤ β k, max . The constants β k,max for k = 1 , 2 , 3 are summarized in Table 4.1.

Table 4.1: β k, max in Theorem 4.1 for k = 1 , 2 , 3 .

| k       |   1 |   2 |     3 |
|---------|-----|-----|-------|
| β k,max |   2 |   1 | 1.243 |

## 4.3 Non periodic boundary conditions

Based on the boundary conditions, we need to modify the treatments for D p L , D p R , and D 0 for non periodic boundary conditions so that the boundary conditions specified on these operators are consistent with the boundary conditions imposed on A .

Firstly, we examine the operators D ∗ for non periodic boundary condition, where ∗ denotes 0, L , or R . Assume that we are given numbers C a and C b .

- · If we require

$$\mathcal { D } _ { L } [ v, \alpha ] ( a ) = C _ { a }, \text{ \ and \ } \mathcal { D } _ { R } [ v, \alpha ] ( b ) = C _ { b },$$

then, the boundary coefficients are obtained as

$$A _ { L } = v ( a ) - C _ { a }, \quad \text{and} \quad B _ { R } = v ( b ) - C _ { b }.$$

- · If we require

$$\mathcal { D } _ { 0 } [ v, \alpha ] ( a ) = C _ { a }, \text{ \ and \ } \mathcal { D } _ { 0 } [ v, \alpha ] ( b ) = C _ { b },$$

then, the boundary coefficients are obtained as

$$A _ { 0 } = \frac { 1 } { 1 - \mu ^ { 2 } } \left ( \mu \left ( I ^ { 0 } [ v, \alpha ] ( b ) - v ( b ) + C _ { b } \right ) - \left ( I ^ { 0 } [ v, \alpha ] ( a ) - v ( a ) + C _ { a } \right ) \right ),$$

$$B _ { 0 } = \frac { 1 } { 1 - \mu ^ { 2 } } \left ( \mu \left ( I ^ { 0 } [ v, \alpha ] ( a ) - v ( a ) + C _ { a } \right ) - \left ( I ^ { 0 } [ v, \alpha ] ( b ) - v ( b ) + C _ { b } \right ) \right ).$$

$$A _ { 0 } & = \frac { 1 } { 1 - \mu ^ { 2 } } \left ( \mu \left ( I ^ { 0 } [ v, \alpha ] ( b ) - v ( b ) + C _ { b } \right ) - \left ( I ^ { 0 } [ v, \alpha ] ( a ) - v ( a ) + C _ { a } \right ) \right ), \quad \quad ( 4. 1 5 \tt a ) \\ B _ { 0 } & = \frac { 1 } { 1 - \mu ^ { 2 } } \left ( \mu \left ( I ^ { 0 } [ v, \alpha ] ( a ) - v ( a ) + C _ { a } \right ) - \left ( I ^ { 0 } [ v, \alpha ] ( b ) - v ( b ) + C _ { b } \right ) \right ). \quad \quad ( 4. 1 5 \tt b )$$

Next, we introduce a modification of the partial sums (4.11) to achieve a higher order accuracy for the non-periodic case. Assume that we have some derivative values at boundary, i.e., ∂ m x A a ( ) and ∂ m x A b ( ), m ≥ 1. Using integration by parts one can derive the following modified partial sums for k ≤ 3 to deal with the non-periodic boundary conditions

$$A _ { x } ^ { - } ( x ) & \approx \mathcal { P } _ { k } ^ { - } [ A, \alpha ] = \begin{cases} & \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { L } [ A _ { 1, p }, \alpha ] ( x ), & k = 1, \, 2, \\ & \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { L } [ A _ { 1, p }, \alpha ] ( x ) - \alpha \mathcal { D } _ { 0 } [ A _ { 1, 3 }, \alpha ] ( x ), & k = 3, \end{cases}$$

$$A _ { x } ^ { + } ( x ) & \approx \mathcal { P } _ { k } ^ { + } [ A, \alpha ] = \begin{cases} \begin{array} { c c } - \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { R } [ A _ { 2, p }, \alpha ] ( x ), & k = 1, \, 2, \\ - \alpha \sum _ { p = 1 } ^ { k } \mathcal { D } _ { R } [ A _ { 2, p }, \alpha ] ( x ) + \alpha \mathcal { D } _ { 0 } [ A _ { 2, 3 }, \alpha ] ( x ), & k = 3. \end{array} \end{cases} \quad ( 4. 1 6 b )$$

And A 1 ,p and A 2 ,p are given as

$$\mathbf l A _ { 1, p } & \text{ and } A _ { 2, p } \text{ are given as} \\ \begin{cases} A _ { 1, 1 } = A, \\ A _ { 1, 2 } = \mathcal { D } _ { L } [ A _ { 1, 1 }, \alpha ] - \sum _ { m = 2 } ^ { k } \left ( - \frac { 1 } { \alpha } \right ) ^ { m } \partial _ { x } ^ { m } A ( a ) e ^ { - \alpha ( x - a ) }, \\ A _ { 1, 3 } = \mathcal { D } _ { L } [ A _ { 1, 2 }, \alpha ] + \sum _ { m = 2 } ^ { k } ( m - 1 ) \left ( - \frac { 1 } { \alpha } \right ) ^ { m } \partial _ { x } ^ { m } A ( a ) e ^ { - \alpha ( x - a ) }, \\ \end{cases} \\ \begin{cases} A _ { 2, 1 } = A, \\ A _ { 2, 2 } = \mathcal { D } _ { R } [ A _ { 2, 1 }, \alpha ] - \sum _ { m = 2 } ^ { k } \left ( \frac { 1 } { \alpha } \right ) ^ { m } \partial _ { x } ^ { m } A ( b ) e ^ { - \alpha ( b - x ) }, \\ A _ { 2, 3 } = \mathcal { D } _ { R } [ A _ { 2, 2 }, \alpha ] + \sum _ { m = 2 } ^ { k } ( m - 1 ) \left ( \frac { 1 } { \alpha } \right ) ^ { m } \partial _ { x } ^ { m } A ( b ) e ^ { - \alpha ( b - x ) }, \\ \end{cases} \\ \text{re the boundary conditions for the operators are imposed as}$$

where the boundary conditions for the operators are imposed as

$$\alpha \mathcal { D } _ { L } [ A _ { 1, 1 }, \alpha ] ( a ) = A _ { x } ( a ), \ \alpha \mathcal { D } _ { R } [ A _ { 2, 1 }, \alpha ] ( b ) = - A _ { x } ( b ),$$

$$\alpha \mathcal { D } _ { L } [ A _ { 1, p }, \alpha ] ( a ) = \alpha \mathcal { D } _ { R } [ A _ { 2, p }, \alpha ] ( b ) = 0, \text{ \ for } p \geq 2,$$

$$& x ( b ), \\ & 2,$$

$$\alpha \mathcal { D } _ { 0 } [ A _ { *, 3 }, \alpha ] ( a ) = \alpha \mathcal { D } _ { 0 } [ A _ { *, 3 }, \alpha ] ( b ) = 0, \ \ * \text{could be $1$ or $2$}.$$

The modified partial sum (4.16) agrees with the derivative values at the boundary

$$\mathcal { P } _ { k } ^ { - } [ A, \alpha ] ( a ) = A _ { x } ( a ), \ \mathcal { P } _ { k } ^ { + } [ A, \alpha ] ( b ) = A _ { x } ( b ).$$

Furthermore, we have the following theorem, which is a result of the Theorem 2.3 from [18].

Theorem 4.2. Suppose A ∈ C k +1 [ a, b ] . If we take α = β/ c ( ∆ ) t , then the modified partial sums (4.16) satisfy

$$\| A _ { x } - \mathcal { P } _ { k } ^ { - } [ A, \alpha ] \| _ { \infty } = \mathcal { O } ( \Delta t ^ { k } ), \ \ \| A _ { x } - \mathcal { P } _ { k } ^ { + } [ A, \alpha ] \| _ { \infty } = \mathcal { O } ( \Delta t ^ { k } ), \ \ k = 1, \, 2, \, 3.$$

As we have seen in (4.17), we need the derivatives ∂ m x A a ( ) and ∂ m x A b ( ), m ≥ 1, for non periodic boundary conditions. Here, we will only focus on the outflow boundary conditions. High order extrapolations are used to get all the derivatives at boundaries (that is, there are no physical boundary given). The details of the general non-periodic conditions can be found in our previous work [18].

## 4.4 Space discretization

According to the MOL T method, a boundary value problem is obtained by discretizing the equation in time first and an integral formulation is presented to solve the boundary value problem. A fully discrete numerical solution is then obtained by discretizing D L and D R operators in space. In this subsection we will give the details of the spatial discretization of the D L and D R operators and WENO-based quadrature formulation to approximate the convolution integrals appear in the D L and D R operators.

Suppose we divide the domain [ a, b ] with N +1 uniformly distributed grid points

$$\Delta x = ( b - a ) / N, \ \ x _ { i } = a + i \Delta x, \ \ i = 0, 1, \dots, N.$$

The convolution integrals I L i = I L [ v, α ]( x i ) and I R i = I R [ v, α ]( x i ) satisfy a recursive relation

$$I _ { i } ^ { L } = e ^ { - \alpha \Delta x _ { i } } I _ { i - 1 } ^ { L } + J _ { i } ^ { L }, \ \ i = 1, \dots, N, \ \ I _ { 0 } ^ { L } = 0,$$

$$I _ { i } ^ { R } = e ^ { - \alpha \Delta x _ { i + 1 } } I _ { i + 1 } ^ { R } + J _ { i } ^ { R }, \ \ i = 0, \dots, N - 1, \ \ I _ { N } ^ { R } = 0,$$

where,

$$J _ { i } ^ { L } = \alpha \int _ { x _ { i - 1 } } ^ { x _ { i } } v ( y ) e ^ { - \alpha ( x _ { i } - y ) } d y, \quad J _ { i } ^ { R } = \alpha \int _ { x _ { i } } ^ { x _ { i + 1 } } v ( y ) e ^ { - \alpha ( y - x _ { i } ) } d y.$$

In [17], we have presented a high order and robust framework to calculate J L i and J R i . Here we will provide a brief description of methodology.

Note that, to approximate J L i with k th order accuracy, we may choose the interpolation stencil S i ( ) = { x i -r , . . . , x i -r + k } , which contains x i -1 and x i . There is a unique polynomial p x ( ) of degree at most k that interpolates v x ( ) at the nodes in S i ( ). Then J L i is approximated by

$$J _ { i } ^ { L } & \approx \alpha \int _ { x _ { i - 1 } } ^ { x _ { i } } p ( y ) e ^ { - \alpha ( x _ { i } - y ) } d y.$$

Note that, the integral on the right hand side can be evaluated exactly. Similarly, we can approximate J R i by

$$J _ { i } ^ { R } & \approx \alpha \int _ { x _ { i } } ^ { x _ { i + 1 } } p ( y ) e ^ { - \alpha ( y - x _ { i } ) } d y,$$

with polynomial p x ( ) interpolating v x ( ) on stencil S i ( ) = { x i + r -k , . . . , x i + r } , which includes x i and x i +1 .

Since the quadratures with a fixed stencil for approximating the J L i and J L i may develop spurious oscillations which violates the entropy solutions, we use WENO-based quadrature formula and the nonlinear filter to control oscillations and capture the correct solution. Such a methodology is discussed in [17] and we will modify that method to solve HJ equations since both equations involve discontinuous derivatives for the solutions and hence, generate the entropy violating solutions.

To summarize, if A x is periodic function, then we will be using the following modified sums framework from [18] instead of (4.11) for approximation of A ± x at x i :

$$A _ { x, i } ^ { - } & = \alpha \mathcal { D } _ { L } [ A, \alpha ] ( x _ { i } ) + \alpha \sum _ { p = 2 } ^ { k } \sigma _ { i, L } ^ { p - 2 } \mathcal { D } _ { L } ^ { p } [ A, \alpha ] ( x _ { i } ), & ( 4. 2 3 \tt a )$$

$$A _ { x, i } ^ { + } = - \alpha \mathcal { D } _ { R } [ A, \alpha ] ( x _ { i } ) - \alpha \sum _ { p = 2 } ^ { k } \sigma _ { i, R } ^ { p - 2 } \mathcal { D } _ { R } ^ { p } [ A, \alpha ] ( x _ { i } ) ;$$

and if A x is non periodic function, then we we will use the following formulation instead of (4.16)

$$A _ { x, i } ^ { - } & = \alpha \mathcal { D } _ { L } [ A _ { 1, 1 }, \alpha ] ( x _ { i } ) + \alpha \sum _ { p = 2 } ^ { k } \sigma _ { i, L } ^ { p - 2 } \mathcal { D } _ { L } [ A _ { 1, p }, \alpha ] ( x _ { i } ), & ( 4. 2 4 a )$$

$$A _ { x, i } ^ { + } = - \alpha \mathcal { D } _ { R } [ A _ { 2, 1 }, \alpha ] ( x _ { i } ) - \alpha \sum _ { p = 2 } ^ { k } \sigma _ { i, R } ^ { p - 2 } \mathcal { D } _ { R } [ A _ { 2, p }, \alpha ] ( x _ { i } ).$$

We only use the WENO formulation when p = 1 while we use cheap high order linear formulation for the case p &gt; 1. The filters σ i,L and σ i,R are obtained based on the smoothness indicators from the WENO quadrature.

The WENO quadrature for J L i is presented here as an example. The related stencil is given in Figure 4.1. We choose the big stencil as S i ( ) = { x i -3 , . . . , x i +2 } and the three small stencils as S r ( ) = i { x i -3+ r , . . . , x i + r } , r = 0 1 2. We will only present the formulas for the case of a uniform , , mesh, i.e., ∆ x i = ∆ x for all i here, but the WENO methodology is still applicable to the case of a nonuniform mesh, see [46].

## 1. We approximate the integrals on each small stencils S r ( ) as follows i

$$J _ { i, r } ^ { L } & = \alpha \int _ { x _ { i - 1 } } ^ { x _ { i } } e ^ { - \alpha ( x _ { i } - y ) } p _ { r } ( y ) d x,$$

where p r ( x ) is the polynomial interpolating v x ( ) on nodes S r ( ). i

- 2. Similarly, on the big stencil S i ( ), we obtain

$$J _ { i } ^ { L } = \alpha \int _ { x _ { i - 1 } } ^ { x _ { i } } e ^ { - \alpha ( x _ { i } - y ) } p ( y ) d x = \sum _ { r = 0 } ^ { 2 } d _ { r } J _ { i, r } ^ { L },$$

Figure 4.1: The structure of the stencils in WENO integration.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-34fd52f3-f286-452a-96c7-f10d6fc7286a.jpg)

with the linear weights d r satisfying ∑ 2 r =0 d r = 1.

- 3. We develop the following nonlinear weights ω r using the linear weights d r

$$\omega _ { r } = \tilde { \omega } _ { r } / \sum _ { s = 0 } ^ { 2 } \tilde { \omega } _ { s }, \ \ r = 0, \ 1, \ 2,$$

with

$$\tilde { \omega } _ { r } = d _ { r } \left ( 1 + \frac { \tau _ { 5 } } { \epsilon + \beta _ { r } } \right ).$$

We take glyph[epsilon1] = 10 -6 as a small positive number, glyph[epsilon1] &gt; 0, in our numerical test problems to avoid zero at the denominator. The smoothness indicator β r is determined as

$$\beta _ { r } = \sum _ { l = 2 } ^ { 3 } \int _ { x _ { i - 1 } } ^ { x _ { i } } \Delta x _ { i } ^ { 2 l - 3 } \left ( \frac { \partial ^ { l } p _ { r } ( x ) } { \partial x ^ { l } } \right ) ^ { 2 } d x,$$

which is used to measure the relative smoothness of the function v x ( ) in the stencil S r ( ). i In particular, we have the expressions as

$$\beta _ { 0 } & = \frac { 1 3 } { 1 2 } ( - v _ { i - 3 } + 3 v _ { i - 2 } - 3 v _ { i - 1 } + v _ { i } ) ^ { 2 } + \frac { 1 } { 4 } ( v _ { i - 3 } - 5 v _ { i - 2 } + 7 v _ { i - 1 } - 3 v _ { i } ) ^ { 2 }, \\ \beta _ { 1 } & = \frac { 1 3 } { 1 2 } ( - v _ { i - 2 } + 3 v _ { i - 1 } - 3 v _ { i } + v _ { i + 1 } ) ^ { 2 } + \frac { 1 } { 4 } ( v _ { i - 2 } - v _ { i - 1 } - v _ { i } + v _ { i + 1 } ) ^ { 2 }, \\ \beta _ { 2 } & = \frac { 1 3 } { 1 2 } ( - v _ { i - 1 } + 3 v _ { i } - 3 v _ { i + 1 } + v _ { i + 2 } ) ^ { 2 } + \frac { 1 } { 4 } ( - 3 v _ { i - 1 } + 7 v _ { i } - 5 v _ { i + 1 } + v _ { i + 2 } ) ^ { 2 }.$$

Here, τ 5 states the absolute difference between β 0 and β 2

$$\tau _ { 5 } = | \beta _ { 0 } - \beta _ { 2 } |.$$

Furthermore, we introduce a parameter ξ i as

$$\xi _ { i } & = \frac { \beta _ { m i n } } { \beta _ { m a x } },$$

which will be use to create the nonlinear filter. Here,

$$\beta _ { m a x } = 1 + \left ( \frac { \tau _ { 5 } } { \epsilon + \min ( \beta _ { 0 }, \beta _ { 2 } ) } \right ) ^ { 2 }, \quad \beta _ { m i n } = 1 + \left ( \frac { \tau _ { 5 } } { \epsilon + \max ( \beta _ { 0 }, \beta _ { 2 } ) } \right ) ^ { 2 }.$$

Note that we have developed the nonlinear weights using the idea of the WENO-Z method proposed in [8], which has less dissipation and higher resolution comparing to the original WENO method.

## 4. Lastly, we obtain the approximation

$$J _ { i } ^ { L } = \sum _ { r = 0 } ^ { 2 } \omega _ { r } J _ { i, r } ^ { L }.$$

The filter σ i,L is determined as

$$\sigma _ { i, L } = \min ( \xi _ { i - 1 }, \xi _ { i } ).$$

The process to obtain J R i and σ i,R is mirror symmetric to that of J L i and σ i,L with respect to point x i .

## 4.5 2D magnetic potential equation

According to the Constrained Transport formulation described in Section 3, we must update the solution of the magnetic potential equation by solving a discrete version of the following equation

$$A _ { t } ^ { 3 } + u ^ { 1 } ( x, y ) A _ { x } ^ { 3 } + u ^ { 2 } ( x, y ) A _ { y } ^ { 3 } = 0$$

where the velocity components u 1 and u 2 are known from the previous time step due to the solution of the base part, HCL. Since the velocity functions are given, we can consider (4.34) as a HamiltonJacobi equation

$$A _ { t } ^ { 3 } + H ( A _ { x } ^ { 3 }, A _ { y } ^ { 3 } ) = 0$$

with Hamiltonian flux

$$H ( A _ { x } ^ { 3 }, A _ { y } ^ { 3 } ) = u ^ { 1 } ( x, y ) A _ { x } ^ { 3 } + u ^ { 2 } ( x, y ) A _ { y } ^ { 3 }.$$

We can directly apply a two dimensional version of the kernel-based framework presented in Section 4.1. The 2D semi discrete scheme can be written as

$$\frac { d A _ { i, j } ^ { 3 } ( t ) } { d t } = - \hat { H } ( A _ { x } ^ { 3 - } | _ { i, j }, A _ { x } ^ { 3 + } | _ { i, j }, A _ { y } ^ { 3 - } | _ { i, j }, A _ { y } ^ { 3 + } | _ { i, j } )$$

on each point ( x , y i j ), where ˆ H is a Lipschitz continuous Hamiltonian flux. We use the LaxFriedrichs flux and the equation (4.37) becomes

$$\frac { d A _ { i, j } ^ { 3 } ( t ) } { d t } = & - u _ { i, j } ^ { 1 } \left ( \frac { A _ { x } ^ { 3 - } + A _ { x } ^ { 3 + } } { 2 } \right ) | _ { i, j } - u _ { i, j } ^ { 2 } \left ( \frac { A _ { y } ^ { 3 - } + A _ { y } ^ { 3 + } } { 2 } \right ) | _ { i, j } \\ & + c _ { 1 } \left ( \frac { A _ { x } ^ { 3 + } - A _ { x } ^ { 3 - } } { 2 } \right ) | _ { i, j } + c _ { 2 } \left ( \frac { A _ { y } ^ { 3 + } - A _ { y } ^ { 3 - } } { 2 } \right ) | _ { i, j }$$

with

$$c _ { 1 } = \max _ { i, j } | u _ { i, j } ^ { 1 } | \text{ \ and \ } c _ { 2 } = \max _ { i, j } | u _ { i, j } ^ { 2 } |.$$

We remark that the scheme (4.38) with this global c m can be very dissipative for some HamiltonJacobi equations. There is another way of choosing c m which is by taking the max value on the local stencil.

The approximations A 3 ± x | i,j and A 3 ± y | i,j to the derivatives of functions A x ( x, y ) and A y ( x, y ) at ( x , y i j ), respectively calculated directly using one dimensional formulation of the scheme, e.g., when computing A 3 ± x , we fix y and apply 1D scheme in x -direction. For example, when A x is periodic in x -direction, we get

$$A _ { x } ^ { - } \approx \alpha _ { 1 } \mathcal { D } _ { L } ^ { p } [ A ( \cdot, y ), \alpha _ { 1 } ] ( x ) + \alpha _ { 1 } \sum _ { p = 2 } ^ { k } \sigma _ { i, L } ^ { p - 2 } \mathcal { D } _ { L } ^ { p } [ A ( \cdot, y ), \alpha _ { 1 } ] ( x ),$$

$$A _ { x } ^ { + } \approx - \alpha _ { 1 } \mathcal { D } _ { R } [ A ( \cdot, y ), \alpha _ { 1 } ] ( x ) - \alpha _ { 1 } \sum _ { p = 2 } ^ { k } \sigma _ { i, R } ^ { p - 2 } \mathcal { D } _ { R } ^ { p } [ A ( \cdot, y ), \alpha _ { 1 } ] ( x ).$$

$$( x ), \\ \\ 1 ] ( x ).$$

Here, we choose α 1 = β/ c ( 1 ∆ ). Similarly, to approximate t A ± y , we fix x and obtain

$$A _ { y } ^ { - } \approx \alpha _ { 2 } \mathcal { D } _ { L } [ A ( x, \cdot ), \alpha _ { 2 } ] ( y ) + \alpha _ { 2 } \sum _ { p = 2 } ^ { k } \sigma _ { i, L } ^ { p - 2 } \mathcal { D } _ { L } ^ { p } [ A ( x, \cdot ), \alpha _ { 2 } ] ( y ),$$

$$A _ { y } ^ { + } \approx - \alpha _ { 2 } \mathcal { D } _ { R } [ A ( x, \cdot ), \alpha _ { 2 } ] ( y ) - \alpha _ { 2 } \sum _ { p = 2 } ^ { k } \sigma _ { i, R } ^ { p - 2 } \mathcal { D } _ { R } ^ { p } [ A ( x, \cdot ), \alpha _ { 2 } ] ( y ),$$

with α 2 = β/ c ( 2 ∆ ). t Note that in the 2D case, we need to choose β max as half of that for the 1D case to ensure the unconditional stability of the scheme.

If it is non periodic boundary condition, we still use extrapolation with suitable order of accuracy for the derivative values at an outflow boundary, as in the 1D formulation. For the details see [18].

## 4.6 3D magnetic potential equation

Although the evolution equation for the 3D magnetic potential (2.12) is significantly different from the evolution equation for the 2D scalar magnetic potential (2.17), we can still directly apply the

scheme presented in section (4.1) to 3D case. In our previous paper [21], we needed an artificial resistivity term in 3D case. However, with the kernel based approach we no longer need those artificial terms in our equations, since our method is fully implicit.

Writing out the magnetic vector potential equation derived in section (2.2) in component form we have:

$$\partial _ { t } A ^ { 1 } = u ^ { 2 } ( \partial _ { x } A ^ { 2 } ) + u ^ { 3 } ( \partial _ { x } A ^ { 3 } ) - u ^ { 2 } ( \partial _ { y } A ^ { 1 } ) - u ^ { 3 } ( \partial _ { z } A ^ { 1 } )$$

$$\partial _ { t } A ^ { 2 } = - u ^ { 1 } ( \partial _ { x } A ^ { 2 } ) + u ^ { 1 } ( \partial _ { y } A ^ { 1 } ) + u ^ { 3 } ( \partial _ { y } A ^ { 3 } ) - u ^ { 3 } ( \partial _ { z } A ^ { 2 } )$$

$$\partial _ { t } A ^ { 3 } = - u ^ { 1 } ( \partial _ { x } A ^ { 3 } ) - u ^ { 2 } ( \partial _ { y } A ^ { 3 } ) + u ^ { 1 } ( \partial _ { z } A ^ { 1 } ) + u ^ { 2 } ( \partial _ { z } A ^ { 2 } )$$

with A = ( A , A , A 1 2 3 ) and u = ( u , u 1 2 , u 3 ). While A in 3D is not strictly a H-J equation, with the new implicit approach we can simply apply the ideas from the 2D case. Then we can obtain the following equations using the Lax-Friedrichs flux splitting:

$$\frac { d A ^ { 1 } ( t ) } { d t } & = \ u ^ { 2 } \frac { ( A _ { x } ^ { 2 + } + A _ { x } ^ { 2 - } ) } { 2 } - c _ { 2 } \frac { ( A _ { x } ^ { 2 - } - A _ { x } ^ { 2 + } ) } { 2 } + u ^ { 3 } \frac { ( A _ { x } ^ { 3 + } + A _ { x } ^ { 3 - } ) } { 2 } - c _ { 3 } \frac { ( A _ { x } ^ { 2 - } - A _ { x } ^ { 3 + } ) } { 2 } } \\ & \quad - u ^ { 2 } \frac { ( A _ { y } ^ { 1 + } + A _ { y } ^ { 1 - } ) } { 2 } - c _ { 2 } \frac { ( A _ { y } ^ { 1 - } - A _ { y } ^ { 1 + } ) } { 2 } - u ^ { 3 } \frac { ( A _ { z } ^ { 1 + } + A _ { z } ^ { 1 - } ) } { 2 } - c _ { 3 } \frac { ( A _ { z } ^ { 1 - } - A _ { z } ^ { 1 + } ) } { 2 }, \quad ( 4. 4 2 a ) \\ d A ^ { 2 } ( t ) & = - u ^ { 1 } \frac { ( A _ { z } ^ { 2 + } + A _ { x } ^ { 2 - } ) } { 2 } - c _ { 1 } \frac { ( A _ { z } ^ { 2 - } - A _ { z } ^ { 2 + } ) } { 2 } + u ^ { 1 } \frac { ( A _ { y } ^ { 1 + } + A _ { y } ^ { 1 - } ) } { 2 } - c _ { 1 } \frac { ( A _ { y } ^ { 1 - } - A _ { y } ^ { 1 + } ) } { 2 } \\ & \quad + u ^ { 3 } \frac { ( A _ { y } ^ { 3 + } + A _ { y } ^ { 3 - } ) } { 2 } - c _ { 3 } \frac { ( A _ { y } ^ { 3 - } - A _ { y } ^ { 3 + } ) } { 2 } - u ^ { 3 } \frac { ( A _ { z } ^ { 2 + } + A _ { z } ^ { 2 - } ) } { 2 } - c _ { 3 } \frac { ( A _ { z } ^ { 2 - } - A _ { z } ^ { 2 + } ) } { 2 }, \quad ( 4. 4 2 b ) \\ d A ^ { 3 } ( t ) & = - u ^ { 1 } \frac { ( A _ { x } ^ { 3 + } + A _ { x } ^ { 3 - } ) } { 2 } - c _ { 1 } \frac { ( A _ { x } ^ { 3 - } - A _ { x } ^ { 3 + } ) } { 2 } - u ^ { 2 } \frac { ( A _ { y } ^ { 3 + } + A _ { y } ^ { 3 - } ) } { 2 } - c _ { 2 } \frac { ( A _ { y } ^ { 3 - } - A _ { y } ^ { 3 + } ) } { 2 } } \\ & \quad + u ^ { 1 } \frac { ( A _ { z } ^ { 1 + } + A _ { z } ^ { 1 - } ) } { 2 } - c _ { 1 } \frac { ( A _ { z } ^ { 1 - } - A _ { z } ^ { 1 + } ) } { 2 } + u ^ { 2 } \frac { ( A _ { z } ^ { 2 + } + A _ { z } ^ { 2 - } ) } { 2 } - c _ { 2 } \frac { ( A _ { z } ^ { 2 - } - A _ { z } ^ { 2 + } ) } { 2 }, \quad ( 4. 4 2 c ) \\ \text{at} \ ( x _ { i }, y _ { j }, z _ { k } ), \text{where}$$

at ( x , y i j , z k ), where

$$c _ { 1 } = \max _ { i, j, k } | u _ { i, j, k } ^ { 1 } |, \ \ c _ { 2 } = \max _ { i, j, k } | u _ { i, j, k } ^ { 2 } |, \ \text{ and } \ c _ { 3 } = \max _ { i, j, k } | u _ { i, j, k } ^ { 3 } |.$$

Similarly, we use 1D kernel-based formulation to approximate the derivatives of the magnetic vector potential components A ∗± x , A ∗± y , and A ∗± z , where ∗ denotes 1 2 and 3. , In addition, since the density or pressure may become negative in some problems such as the blast wave problem, we use the positivity preserving limiter idea developed in our previous work [21]. For more details on positivity limiter, see [19, 20, 45]

## 5 Central finite difference discretization of ∇× A

According to the Constrained Transport formulation, we applied a discrete curl operator to the magnetic potential at each stage of CT framework for the purpose of obtaining a divergence free magnetic field. Here we will present the strategy we used to approximate the curl operator.

## 5.1 Curl in 2D

We apply a discrete version of 2D curl to the equation (2.16) as

$$B ^ { 1 } _ { i, j } \coloneqq D _ { y } | _ { i, j } A ^ { 3 } \quad \text{and} \quad B ^ { 2 } _ { i, j } \coloneqq - D _ { x } | _ { i, j } A ^ { 3 },$$

where D x and D y are the discrete versions of the partial derivatives ∂ x and ∂ y . We use these version of curls to satisfy the discrete divergence free constraint as follows

$$\nabla \cdot B _ { i, j } \coloneqq D _ { x } | _ { i, j } B ^ { 1 } + D _ { y } | _ { i, j } B ^ { 2 } = D _ { x } D _ { y } | _ { i, j } A ^ { 3 } - D _ { y } D _ { x } | _ { i, j } A ^ { 3 } = 0.$$

Here, we choose 4th-order central finite differences to get high order accuracy

$$D _ { x } | _ { i, j } A ^ { 3 } \coloneqq \frac { 1 } { 1 2 \Delta x } ( A _ { i - 2, j } ^ { 3 } - 8 A _ { i - 1, j } ^ { 3 } + 8 A _ { i + 1, j } ^ { 3 } - A _ { i + 2, j } ^ { 3 } ),$$

$$D _ { y } | _ { i, j } A ^ { 3 } \coloneqq \frac { \dashv \dashv _ { 1 } } { 1 2 \Delta y } ( A ^ { 3 } _ { i, j - 2 } - 8 A ^ { 3 } _ { i, j - 1 } + 8 A ^ { 3 } _ { i, j + 1 } - A ^ { 3 } _ { i, j + 2 } ),$$

where A = A 3 is only third component of magnetic potential.

## 5.2 Curl in 3D

We use the following discrete version of the 3D curl

$$B ^ { 1 } _ { i, j, k } \coloneqq D _ { y } | _ { i, j, k } A ^ { 3 } - D _ { z } | _ { i, j, k } A ^ { 2 },$$

$$B _ { i, j, k } ^ { 2 } \coloneqq D _ { z } | _ { i, j, k } A ^ { 1 } - D _ { x } | _ { i, j, k } A ^ { 3 },$$

$$B _ { i, j, k } ^ { 3 } \coloneqq D _ { x } | _ { i, j, k } A ^ { 2 } - D _ { y } | _ { i, j, k } A ^ { 1 },$$

where D x , D y , and D z are notations for the discrete versions of partial derivatives ∂ x , ∂ y and ∂ z , respectively. With these versions of curls, we get the following divergence free condition satisfied

$$\text{respective} \text{, } \text{with these versions of curs, we get the following divergence tree contain assume} \\ \nabla \cdot B _ { i, j, k } & \colon = D _ { x | i, j, k } B ^ { 1 } + D _ { y | i, j, k } B ^ { 2 } + D _ { z | i, j, k } B ^ { 3 } \\ & = D _ { x } D _ { y | i, j, k } A ^ { 3 } - D _ { x } D _ { z | i, j, k } A ^ { 2 } + D _ { y } D _ { z | i, j, k } A ^ { 1 } - D _ { y } D _ { x | i, j, k } A ^ { 3 } \\ & + D _ { z } D _ { x | i, j, k } A ^ { 2 } - D _ { z } D _ { y | i, j, k } A ^ { 1 } \\ & = 0$$

As in 2D case, we use 4th-order central finite differences to get high order accuracy

$$D _ { x } | _ { i, j, k } A ^ { * } \coloneqq \frac { 1 } { 1 2 \Delta x } ( A ^ { * } _ { i - 2, j, k } - 8 A ^ { * } _ { i - 1, j, k } + 8 A ^ { * } _ { i + 1, j, k } - A ^ { * } _ { i + 2, j, k } ),$$

$$D _ { y } | _ { i, j, k } A ^ { * } \coloneqq \frac { \frac { 1 } { 1 } } { _ { 2 \Delta y } } ( A ^ { * } _ { i, j - 2, k } - 8 A ^ { * } _ { i, j - 1, k } + 8 A ^ { * } _ { i, j + 1, k } - A ^ { * } _ { i, j + 2, k } ),$$

$$D _ { z } | _ { i, j, k } A ^ { * } \coloneqq \frac { 1 } { 1 2 \Delta z } ( A ^ { * } _ { i, j, k - 2 } - 8 A ^ { * } _ { i, j, k - 1 } + 8 A ^ { * } _ { i, j, k + 1 } - A ^ { * } _ { i, j, k + 2 } ).$$

## 6 Numerical Results

In this section, we present the numerical results to demonstrate the accuracy and efficiency of the new method. We use third order SSP RK method for time discretization. Time step is chosen as

$$\Delta t & = \frac { \text{CFL} } { e m * c d }$$

where, the CFL number is 0.5, em = max( λ , λ 5 6 ), cd = 1 ∆ x + 1 ∆ y for 2D and cd = 1 ∆ x + 1 ∆ y + 1 ∆ z for 3D.

## 6.1 Smooth vortex test in MHD

We first test the smooth vortex problem in 2D with non zero magnetic field to show the accuracy of the method within the constrained transport formulation.The initial conditions are

$$( \rho, u ^ { 1 }, u ^ { 3 }, u ^ { 3 }, p, B ^ { 1 }, B ^ { 2 }, B ^ { 3 } ) = ( 1, 1, 1, 0, 1, 0, 0, 0 )$$

with perturbations on u , u 1 2 , B 1 , B 2 and p as :

$$\varphi _ { 1 } & \quad \varphi _ { 2 } \\ & ( \delta u ^ { 1 }, \delta u ^ { 2 } ) = \frac { \mu } { 2 \pi } e ^ { 0. 5 ( 1 - r ^ { 2 } ) } ( - y, x ). \\ & ( \delta B ^ { 1 }, \delta B ^ { 2 } ) = \frac { \kappa } { 2 \pi } e ^ { 0. 5 ( 1 - r ^ { 2 } ) } ( - y, x ). \\ & \delta p = \frac { \mu ^ { y } ( 1 - r ^ { 2 } ) - \kappa ^ { 2 } } { 8 \pi ^ { 2 } } e ^ { ( 1 - r ^ { 2 } ) }.$$

$$8 \pi ^ { \colon }$$

And the initial condition for magnetic potential is

$$A ^ { 3 } ( 0, x, y ) = \frac { \mu } { 2 \pi } e ^ { 0. 5 ( 1 - r ^ { 2 } ) }$$

where r 2 = x 2 + y 2 . The vortex strength is taken as µ = 5 389489439 and . κ = √ 2 µ such that the lowest pressure is around 5 3 . × 10 -12 which happens in the center of the vortex. The domain is [ -10 10] , × -[ 10 10] and periodic boundary condition is used on all four boundaries. , In Table 6.2, we present the errors of ρ at t = 0 05 with the mesh size 160 . × 160, demonstrating that the scheme is third order as designed.

Table 6.2: Smooth vortex problem. Errors of ρ and orders of accuracy.

| N x × N y   |   L 1 error | order   |   L ∞ error | order   |
|-------------|-------------|---------|-------------|---------|
| 20 × 20     |   0.002827  | -       |   0.1479    | -       |
| 40 × 40     |   0.0002982 | 3.245   |   0.01839   | 3.007   |
| 80 × 80     |   1.861e-05 | 4.003   |   0.001327  | 3.793   |
| 160 × 160   |   1.102e-06 | 4.078   |   0.0001126 | 3.559   |
| 320 × 320   |   7.26e-08  | 3.924   |   1.17e-05  | 3.267   |

## 6.2 2D Orszag-Tang Vortex

The Orszang-Tang vortex problem is a standard model problem for testing ∇· B = 0 condition, since the solution of the problem is sensitive to divergence errors at late times. The initial conditions are given as

$$( \rho, u ^ { 1 }, u ^ { 2 }, u ^ { 3 }, p, B ^ { 1 }, B ^ { 2 }, B ^ { 3 } ) = ( \gamma ^ { 2 }, - \sin ( y ), \sin ( x ), 0, \gamma, - \sin ( y ), \sin ( 2 x ), 0 ),$$

where γ = 5 3 is the ideal gas constant and the initial magnetic potential: /

$$A ^ { 3 } = 0. 5 \cos ( 2 x ) + \cos ( y )$$

The computational domain is [0 2 , π ] × [0 , 2 π ] and periodic boundary conditions are used everywhere. We test the schemes with 192 × 192 grid points. Although the problem has smooth initial condition, as solution progresses, several shock waves and a vortex shape appear in the center of the domain. In Figure 6.2, we show density ρ at time t = 3, and compare the results with our previous method [21]. We can see that they are in good agreement.

## 6.3 Cloud Shock

In this section we consider the 2D cloud-shock interaction problem, which models a strong shock passing through a dense stationary bubble. The initial conditions include

$$( \rho, u ^ { 1 }, u ^ { 2 }, u ^ { 3 }, p, B ^ { 1 }, B ^ { 2 }, B ^ { 3 } ) = \begin{cases} \ ( 3. 8 6 8 5 9, 1 1. 2 5 3 6, 0, 0, 1 6 7. 3 4 5, 0, 2. 1 8 2 6 1 8 2, - 2. 1 8 2 6 1 8 2 ) & x < 0. 0 5, \\ ( 1, 0, 0, 0, 1, 0, 0. 5 6 4 1 8 9 5 8, 0. 5 6 4 1 8 9 5 8 ) & x > 0. 0 5, \end{cases}$$

and a circular cloud of density ρ = 10 and radius r = 0 15 centered at ( . x, y ) = (0 25 0 5). . , . The computational domain is [0 1] , × [0 , 1] with mesh 512 × 512. We use inflow boundary condition at left boundary and outflow boundary condition elsewhere. The initial condition for magnetic potential:

$$A ^ { 3 } = \begin{cases} \ - 2. 1 8 2 6 1 8 2 ( x - 0. 0 5 ) & x \leq 0. 0 5, \\ - 0. 5 6 4 1 8 9 5 8 ( x - 0. 0 0 5 ) & x \geq 0. 0 5, \end{cases}$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-61f852e3-3ae1-4216-af81-a5e52780f65b.jpg)

(a) Kernal-based method.

Figure 6.2: Orszag-Tang vortex problem. Contour plots of density at t = 3 with 192 × 192 grid points.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-a122635f-57ea-4095-9b6a-fe7d33e7b180.jpg)

Figure 6.3 presents Schlieren plots of || B || at t = 0 06. . The new method matches well with our numerical results of our previous method [21].

## 6.4 2D Blast wave

Here we investigate the 2D blast wave test problem, which has a strong shock causing negative density or pressure in simulations. To avoid negative density or pressure, we use the positivity preserving limiter we developed in reference [21]. The initial conditions are

$$( \rho, u ^ { 1 }, u ^ { 2 }, u ^ { 3 }, B ^ { 1 }, B ^ { 2 }, B ^ { 3 } ) = ( 1, 0, 0, 0, 5 0 / \sqrt { 2 \pi }, 5 0 / \sqrt { 2 \pi }, 0 )$$

with a spherical pressure pulse

$$p = \begin{cases} \ 1 0 0 0 & r \leq 0. 1, \\ 0. 1 & \text{otherwise}. \end{cases}$$

The initial condition for magnetic potential can be obtained as

$$A ^ { 3 } = ( 0, 0, 5 0 / \sqrt { 2 \pi } ( y - x ) ).$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-8489186a-5552-42a3-9fd6-b42af1ec4326.jpg)

(a) Kernal-based method.

Figure 6.3: Cloud shock problem. Contour plots of ‖ B ‖ at t = 0 06 . with 512 × 512 grid points.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-03ecec73-e775-43ce-9d5a-3752070ecd1c.jpg)

We use a domain [ -0 5 . , 0 5] . × -[ 0 5 . , 0 5] with 256 . × 256 mesh. Outflow boundary conditions are applied everywhere. The results are shown in Figure 6.4 and they have good agreement with the results that we got in our previous paper.

## 6.5 2D Field Loop

The field loop problem is a strenuous test problem that moves a steady state. Here we test our recent method on this strenuous test case. The initial conditions are

$$( \rho, u ^ { 1 }, u ^ { 2 }, p ) = ( 1, \sqrt { 5 } \cos ( \theta ), \sqrt { 5 } \sin ( \theta ), 1 )$$

with the angle θ = arctan(0 5). The initial conditions for magnetic field . B are determined by taking the curl of the magnetic potential, which is initialized with

$$A ^ { 3 } = \begin{cases} \begin{array} { c c } 0. 0 0 1 ( R - r ), & r \leq R, \\ 0, & \text{otherwise}, \end{array} \end{cases}$$

where r = √ x 2 + y 2 and R = 0 3. . We use a domain [ -0 5 . , 0 5] . × -[ 0 5 . , 0 5] with 128 . × 128 mesh. Periodic boundary conditions are applied to all sides. The loop is started in the middle of the domain and advected around the grid once, returned to initial location at time t = 1 as shown in Figure 6.5. We observe that the solution maintains the circular symmetry of the initial condition.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-957a1619-8618-4bfd-8744-a77d5d8776f0.jpg)

(a) density ρ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-b804bc52-647d-4167-9269-c991f742f007.jpg)

(b) pressure p

Figure 6.4: 2D Blast wave problem. Contour plots at time t = 0 01 . with 256 × 256 grid points. (a) density; (b) pressure; (c) the norm of u ; (d) magnetic pressure ‖ B ‖ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-e4840a96-ab87-432f-a279-eba976fd9397.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-e57fe683-d62c-4e31-85e7-5550f7e41c19.jpg)

Figure 6.5: 2D Field loop. Contour plots at time t = 1 with 128 × 128 grid points. a) Magnetic Pressure ‖ B ‖ ; b) Magnetic Potential A 3 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-9604ac02-2a41-465a-aa2e-244bdcce5e4f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-6d8449d5-71be-4c56-94b0-fafd75d5d318.jpg)

## 6.6 3D Field Loop

We tested an advecting field loop which moved diagonally across the boundary with an arbitrary initial angle. The initial conditions are

$$( \rho, u ^ { 1 }, u ^ { 2 }, u ^ { 3 }, p ) = ( 1, 2 / \sqrt { 6 }, 1 / \sqrt { 6 }, 1 / \sqrt { 6 }, 1 )$$

The initial conditions for magnetic field are determined by taking the curl of the magnetic potential, which is given as magnetic potential :

$$A ^ { 3 } = \begin{cases} \begin{array} { c c } 0. 0 0 1 ( R - r ) & r \leq R \\ 0 & \text{otherwise}, \end{array} \end{cases}$$

where A 1 = 0, A 2 = 0, and r = √ x 2 + y 2 and R = 0 3. . We use a domain size [ -0 5 . , 0 5] . × [ -0 5 . , 0 5] . × -[ 0 5 . , 0 5] . with 128 × 128 × 128 mesh. Periodic boundary conditions are applied to all sides. We observe that the field loop integrity is maintained, after advecting diagonally around the domain, until the final time, t = 1. The results shown in Figure 6.6 are a 2D slice of the 3D solution taken at z = 0.

Figure 6.6: 3D Field loop. Contour plots at time t = 1 with 128 × 128 grid points. The loop has been advected around the grid once. a) Magnetic Pressure; b) Magnetic Potential.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-da90ed0c-52cf-4a60-8ffd-4622cd113a6a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-d37b6a52-7240-4898-8a35-fe3a23ef0777.jpg)

## 6.7 3D Blast wave

In this section we investigate the 3D version of the blast wave problem to show the strength of the new method, which eliminating the need for the diffusion limiter. The initial conditions are

$$( \rho, u ^ { 1 }, u ^ { 2 }, u ^ { 3 }, B ^ { 1 }, B ^ { 2 }, B ^ { 3 } ) = ( 1, 0, 0, 0, 5 0 / \sqrt { 2 \pi }, 5 0 / \sqrt { 2 \pi }, 0 )$$

with a spherical pressure pulse

$$p = \begin{cases} \ 1 0 0 0 & r \leq 0. 1 \\ 0. 1 & \text{otherwise}. \end{cases}$$

where r = √ x 2 + y 2 + z 2 . The initial condition for magnetic potential :

$$\mathbf A ( 0, x, y, z ) = ( 0, 0, 5 0 / \sqrt { 2 \pi } y - 5 0 / \sqrt { 2 \pi } x ).$$

We use a domain size [ -0 5 . , 0 5] . × -[ 0 5 . , 0 5] . × -[ 0 5 . , 0 5] . with 150 × 150 × 150 mesh. Outflow boundary conditions are applied everywhere.The results presented in Figure 6.7 and Figure 6.8 are the solutions cut at z = 0. The sub plots highlight the differences between the two methods. The key differences are the new method has a maximum value that achieves a peak of 320, which is

23% higher than the old method, and is far less isotropic around the peek than the old method. As demonstrated in a latter test, the isotropic behavior and lower the peek in the strong bast wave is due to the diffusion limiter that was needed with the old explicit constrained transport method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-5227757e-a023-4318-bbd2-7e8d4a0b615e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-ec6da766-db37-4722-aad9-5d0236841c2a.jpg)

- (a) density with kernel-based method

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-c9b1865a-5077-46ff-a1c6-f6ab9f45f02d.jpg)

(b) density with previous method [21].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-2635e035-45e4-4ac7-bfd8-0f5954011066.jpg)

- (c) pressure with kernel-based method

(d) pressure with previous method [21]

Figure 6.7: 3D Blast wave problem. Contour plots at time t = 0 01 . with 150 × 150 grid points

## 6.8 Kernel-based method with Diffusion Terms

In this section we demonstrate the advantage of the new approach in the context of strong shocks. In the previous approach [21], we needed to add in a diffusion limiter to stabilize the update for A in the vicinity of a sharp change in the solution. This had no impact on smooth solutions, and minimal impact on problems like the cloud shock. But in the context of the blast wave problem, while the diffusion limiter stabilized the solution, comparing ‖ u ‖ for both schemes, we have the max value of the previous method is 261, while it is 320 for the kernel-based method. This involved decreasing the maximum values, making the solution more isotropic, and changing the structure of the contours away from the blast. In Figure 6.9, we show the solution generated by the new method with the diffusion limiter used in our previous code. In this case, the max value of ‖ u ‖ is 265. These results confirm that the limiter is what is causing a change in the solution structure for problems with strong shocks. Adding the diffusion limiter caused the solution to revert from the results of the new method to those obtained by the old method.

## 7 Conclusion

In this work we developed a kernel based constrained transpose scheme based on the magnetic vector potential equations in 2D and 3D for ensuring that the solution of the magnetic filed in ideal MHD is divergence free. The development of the method relies on a kernel based formulation of the spatial derivatives. The framework of the current method is derived from the method of lines transpose methodology and the key idea of successive convolution. The method relies on the idea of replacing a local operator with a global operator that is as efficient as an explicit method, but because it is global it provides an unconditional stable method when coupled with explicit time stepping. For time integration, we coupled the method with the high order explicit strong stability preserving Runge-Kutta scheme. The most important conclusion of this work is that the newly proposed method offers an approach to constrained transport that is mesh aligned for AMR and does not rely on diffusion limiter for stability in 3D, which we needed in our previous work [21]. Eliminating the need of diffusion limiter improves solutions where there are strong shocks, and as demonstrated in numerical simulations, the approach also works well for smooth problems. In addition, the new method is unconditionally stable and achieves high order accuracy. The method is robust and has been tested on a range of 2D and 3D test problems, such as the field loop and

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-7851d5d2-2090-4ca0-9f15-080aa74027ed.jpg)

(a) ‖ u ‖ with kernel-based method

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-3b4f63d7-756b-493b-a793-f7c12fe36c74.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-30cd4c5b-5cb1-4ea8-bc8d-b662d53e54c0.jpg)

(b) ‖ u ‖ with previous method [21]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-9d4963ea-4aab-45e8-89a4-ba29ab72c89b.jpg)

(c)

‖

B

‖

with kernel-based method

(d)

‖

B

‖

with previous method [21]

Figure 6.8: 3D Blast wave problem. Contour plots at time t = 0 01 . with 150 × 150 grid points. While the max ‖ u ‖ value for previous method is 261, it is 320 for the kernel-based method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-048729a7-2b9e-4b38-9d6c-8a6963075130.jpg)

(a) density

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-ca89362b-dce1-4fa1-9a3a-390faee30b77.jpg)

(b) pressure

Figure 6.9: 3D Blast wave with diffusion. The max ‖ u ‖ value is 265.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-4d92f01f-6d49-4e45-b425-41996b6d36bf.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01023v1-figure-69f8a263-c9e7-4a1e-98f8-25b8841973f4.jpg)

blast wave problems.

## References

- [1] D. S. Balsara. Divergence-free adaptive mesh refinement for magnetohydrodynamics. Journal of Computational Physics , 174(2):614-648, 2001.
- [2] D. S. Balsara. Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction. The Astrophysical Journal Supplement Series , 151(1):149, 2004.
- [3] D. S. Balsara. Divergence-free reconstruction of magnetic fields and weno schemes for magnetohydrodynamics. Journal of Computational Physics , 228(14):5040-5056, 2009.
- [4] D. S. Balsara and J. Kim. A comparison between divergence-cleaning and staggered-mesh formulations for numerical magnetohydrodynamics. The Astrophysical Journal , 602(2):1079, 2004.
- [5] D. S. Balsara, C. Meyer, M. Dumbser, H. Du, and Z. Xu. Efficient implementation of ader schemes for euler and magnetohydrodynamical flows on structured meshes-speed comparisons with runge-kutta methods. Journal of Computational Physics , 235:934-969, 2013.
- [6] D. S. Balsara, T. Rumpf, M. Dumbser, and C.-D. Munz. Efficient, high accuracy ader-weno schemes for hydrodynamics and divergence-free magnetohydrodynamics. Journal of Computational Physics , 228(7):2480-2516, 2009.
- [7] D. S. Balsara and D. S. Spicer. A staggered mesh algorithm using high order godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations. Journal of Computational Physics , 149(2):270-292, 1999.
- [8] R. Borges, M. Carmona, B. Costa, and W. S. Don. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of Computational Physics , 227(6):3191-3211, 2008.
- [9] J. U. Brackbill and D. C. Barnes. The effect of nonzero? · b on the numerical solution of the magnetohydrodynamic equations. Journal of Computational Physics , 35(3):426-430, 1980.
- [10] M. Causley, H. Cho, and A. Christlieb. Method of lines transpose: Energy gradient flows using direct operator inversion for phase-field models. SIAM Journal on Scientific Computing , 39(5):B968-B992, 2017.

- [11] M. Causley, A. Christlieb, B. Ong, and L. Van Groningen. Method of lines transpose: An implicit solution to the wave equation. Mathematics of Computation , 83(290):2763-2786, 2014.
- [12] M. F. Causley, H. Cho, A. J. Christlieb, and D. C. Seal. Method of lines transpose: High order L-stable O ( N ) schemes for parabolic equations using successive convolution. SIAM Journal on Numerical Analysis , 54(3):1635-1652, 2016.
- [13] M. F. Causley and A. J. Christlieb. Higher order a-stable schemes for the wave equation using a successive convolution approach. SIAM Journal on Numerical Analysis , 52(1):220-235, 2014.
- [14] M. F. Causley, A. J. Christlieb, Y. Guclu, and E. Wolf. Method of lines transpose: A fast implicit wave propagator. arXiv preprint arXiv:1306.6902 , 2013.
- [15] Y. Cheng, F. Li, J. Qiu, and L. Xu. Positivity-preserving dg and central dg methods for ideal mhd equations. Journal of Computational Physics , 238:255-280, 2013.
- [16] A. Christlieb, W. Guo, and Y. Jiang. A WENO-based Method of Lines Transpose approach for Vlasov simulations. Journal of Computational Physics , 327:337-367, 2016.
- [17] A. Christlieb, W. Guo, and Y. Jiang. Kernel based high order 'explicit' unconditionallystable scheme for nonlinear degenerate advection-diffusion equations. arXiv preprint arXiv:1707.09294 , 2017.
- [18] A. Christlieb, W. Guo, and Y. Jiang. A kernel based high order 'explicit' unconditionally stable scheme for time dependent hamilton-jacobi equations. Journal of Computational Physics , 379:214-236, 2019.
- [19] A. J. Christlieb, X. Feng, D. C. Seal, and Q. Tang. A high-order positivity-preserving singlestage single-step method for the ideal magnetohydrodynamic equations. Journal of Computational Physics , 316:218-242, 2016.
- [20] A. J. Christlieb, Y. Liu, Q. Tang, and Z. Xu. Positivity-preserving finite difference weighted eno schemes with constrained transport for ideal magnetohydrodynamic equations. SIAM Journal on Scientific Computing , 37(4):A1825-A1845, 2015.

- [21] A. J. Christlieb, J. A. Rossmanith, and Q. Tang. Finite difference weighted essentially nonoscillatory schemes with constrained transport for ideal magnetohydrodynamics. Journal of Computational Physics , 268:302-325, 2014.
- [22] W. Dai and P. R. Woodward. A simple finite difference scheme for multidimensional magnetohydrodynamical equations. Journal of Computational Physics , 142(2):331-369, 1998.
- [23] H. De Sterck. Multi-dimensional upwind constrained transport on unstructured grids for 'shallow water' magnetohydrodynamics. In AIAA Computational Fluid Dynamics Conference, 15 th, Anaheim, CA , 2001.
- [24] A. Dedner, F. Kemm, D. Kr¨ oner, C.-D. Munz, T. Schnitzer, and M. Wesenberg. Hyperbolic divergence cleaning for the mhd equations. Journal of Computational Physics , 175(2):645-673, 2002.
- [25] C. R. DeVore. Flux-corrected transport techniques for multidimensional compressible magnetohydrodynamics. Journal of Computational Physics , 92(1):142-160, 1991.
- [26] E. Dorfi. Numerical methods for astrophysical plasmas. Computer Physics Communications , 43(1):1-15, 1986.
- [27] C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows-a constrained transport method. The Astrophysical Journal , 332:659-677, 1988.
- [28] M. Fey and M. Torrilhon. A constrained transport upwind scheme for divergence-free advection. In Hyperbolic problems: theory, numerics, applications , pages 529-538. Springer, 2003.
- [29] P. Fu, F. Li, and Y. Xu. Globally divergence-free discontinuous galerkin methods for ideal magnetohydrodynamic equations. Journal of Scientific Computing , 77(3):1621-1659, 2018.
- [30] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. SIAM review , 43(1):89-112, 2001.
- [31] C. Helzel, J. A. Rossmanith, and B. Taetz. An unstaggered constrained transport method for the 3d ideal magnetohydrodynamic equations. Journal of Computational Physics , 230(10):3803-3829, 2011.

- [32] C. Helzel, J. A. Rossmanith, and B. Taetz. A high-order unstaggered constrained-transport method for the three-dimensional ideal magnetohydrodynamic equations based on the method of lines. SIAM Journal on Scientific Computing , 35(2):A623-A651, 2013.
- [33] G.-S. Jiang and D. Peng. Weighted ENO schemes for Hamilton-Jacobi equations. SIAM Journal on Scientific computing , 21(6):2126-2143, 2000.
- [34] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes. Journal of computational physics , 126(1):202-228, 1996.
- [35] S. Kawai. Divergence-free-preserving high-order schemes for magnetohydrodynamics: An artificial magnetic resistivity method. Journal of Computational Physics , 251:292-318, 2013.
- [36] R. J. LeVeque. Wave propagation algorithms for multidimensional hyperbolic systems. Journal of Computational Physics , 131(2):327-353, 1997.
- [37] F. Li, L. Xu, and S. Yakovlev. Central discontinuous galerkin methods for ideal mhd equations with the exactly divergence-free magnetic field. Journal of Computational Physics , 230(12):4828-4847, 2011.
- [38] P. Londrillo and L. Del Zanna. High-order upwind schemes for multidimensional magnetohydrodynamics. The Astrophysical Journal , 530(1):508, 2000.
- [39] P. Londrillo and L. Del Zanna. On the divergence-free condition in godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method. Journal of Computational Physics , 195(1):17-48, 2004.
- [40] G. K. Parks. Physics of space plasmas-an introduction. Redwood City, CA, Addison-Wesley Publishing Co., 1991, 547 p. , 1991.
- [41] K. G. Powell. An approximate riemann solver for magnetohydrodynamics (that works more than one dimension). Technical report, INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING HAMPTON VA, 1994.
- [42] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw. A solutionadaptive upwind scheme for ideal magnetohydrodynamics. Journal of Computational Physics , 154(2):284-309, 1999.

- [43] J. A. Rossmanith. An unstaggered, high-resolution constrained transport method for magnetohydrodynamic flows. SIAM Journal on Scientific Computing , 28(5):1766-1797, 2006.
- [44] D. Ryu, F. Miniati, T. Jones, and A. Frank. A divergence-free upwind code for multidimensional magnetohydrodynamic flows. The Astrophysical Journal , 509(1):244, 1998.
- [45] D. C. Seal, Q. Tang, Z. Xu, and A. J. Christlieb. An explicit high-order single-stage singlestep positivity-preserving finite difference weno method for the compressible euler equations. Journal of Scientific Computing , 68(1):171-190, 2016.
- [46] C.-W. Shu. High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM review , 51(1):82-126, 2009.
- [47] M. Torrilhon. Locally divergence-preserving upwind finite volume schemes for magnetohydrodynamic equations. SIAM Journal on Scientific Computing , 26(4):1166-1191, 2005.
- [48] G. T´ oth. The ∇· B = 0 constraint in shock-capturing magnetohydrodynamics codes. Journal of Computational Physics , 161(2):605-652, 2000.
- [49] J. R. Wilson. Some magnetic effects in stellar collapse and accretion fn1. Annals of the New York Academy of Sciences , 262(1):123-132, 1975.
- [50] K. Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. IEEE Transactions on antennas and propagation , 14(3):302-307, 1966.
- [51] A. L. Zachary, A. Malagoli, and P. Colella. A higher-order godunov method for multidimensional ideal magnetohydrodynamics. SIAM Journal on Scientific Computing , 15(2):263-284, 1994.
Firat Cakir Andrew Christlieb Yan Jiang
A Kernel Based High Order "Explicit" Unconditionally Stable Constrained Transport Method for Ideal Magnetohydrodynamics
The ideal Magnetohydrodynamics (MHD) equations are challenging because one
needs to maintain the divergence free condition, $\nabla \cdot \Bv = 0$. Many
numerical methods have been developed to enforce this condition. In this work,
we further our work on mesh aligned constrained transport by developing a new
kernel based approach for the vector potential in 2D and 3D. The approach for
solving the vector potential is based on the method of lines transpose and is
A-stable, eliminating the need for diffusion limiters needed in our previous
work in 3D. The work presented here is an improvement over the previous method
in the context of problems with strong shocks due to the fact that we could
eliminate the diffusion limiter that was needed in our previous version of
constrained transport. The method is robust and has been tested on the 2D and
3D cloud shock, blast wave and field loop problems.