A KERNEL BASED HIGH ORDER EXPLICIT UNCONDITIONALLY STABLE

CONSTRAINED TRANSPORT METHOD FOR IDEAL AND NON-IDEAL

MAGNETOHYDRODYNAMICS

By

Firat Cakir

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Mathematics — Doctor of Philosophy

2020

ABSTRACT

A KERNEL BASED HIGH ORDER EXPLICIT UNCONDITIONALLY STABLE

CONSTRAINED TRANSPORT METHOD FOR IDEAL AND NON-IDEAL

MAGNETOHYDRODYNAMICS

By

Firat Cakir

The Magnetohydrodynamics (MHD) model is an important model to describe the be-

havior of the plasma phenomenon. The 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 thesis we propose a hybrid MHD solver via the

constrained transport method which is based on kernel based approximations to correct the

solution of magnetic field equation so that it satisfies the divergence free condition. The main

focus of this work is on the constrained transport methodology that intends to eliminate the

need of artificial diffusion limiters (introduced in the previous work [29]) for the magnetic

vector potential equation for ideal MHD and extend the proposed scheme to resistive MHD,

which includes a term that takes the form of diffusion.

In the first part of the thesis, we focus on ideal MHD model and we further our work on

mesh aligned constrained transport [29] 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.

In the second part of the thesis, we turn our focus to resistive MHD model, which allows

magnetic diffusion terms in the magnetic field induction equation. We derive magnetic

potential equation with diffusion terms and develop our kernel based method for the second

order derivative terms appearing in the magnetic potential equation. We update the magnetic

field solutions using the constrained transport framework based on our novel kernel based

scheme. Our novel scheme is robust and unconditionally stable. Numerical examples in 2D

are provided to verify the order of accuracy and to test the performance of the proposed

scheme.

ACKNOWLEDGMENTS

First and foremost, I would like to express my deepest and sincere gratitude to my advisor,

Dr. Andrew Christlieb, for his continuous support and encouragement throughout my Ph.D

study and related research. Andrew provided me tremendous knowledge, great insight,

motivation, patient guidance and pleasant research environment. His guidance helped me in

all the time of research and writing this thesis. I could not have imagined having a better

advisor and mentor for my Ph.D study. It was also truly helpful to have weekly meetings

with my advisor and other group members to discuss over the works with them. Besides my

advisor, I would like to thank to Yan Jiang, Hyoseon Yang for their close collaborations and

to William Sand for the help he offered me. I would also like to thank to my dissertation

committee: Dr. Yingda Cheng, Dr. Qian Jianliang, and Dr. Brian O’Shea for their insightful

comments and encouragement.

I am eternally grateful to my wife, Hayriye, for her support, encouragement, patience

and endless love she has for me. I wouldn’t have made it without her support, help and love.

I always consider her as my life mentor. I love you my beautiful wife! I also want to express

my deepest love to my lovely kids, Miran Aras and Nehir Lavin. Their existence made my

life more meaningful. I also want to thank to the rest of my family: my parents, my parents

in law, my siblings and my siblings in law, for their endless love, support and prayers. They

have always been there for me and I have always felt their love in my heart. I am truly

thankful to my mom for her continuous praying and love to me. Moreover, I would like to

thank my special friends Ayse and Metin for their great friendship and support.

Last but not least, I am beyond words in my gratitude to Allah.

iv

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Plasma Phenomenon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Theoretical Overview of Plasma
. . . . . . . . . . . . . . . . . . . . . . . .
1.4 Review of Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 Numerical Methods for Hyperbolic Conversation Laws
. . . . . . . .
1.4.2 Numerical Methods for ideal MHD . . . . . . . . . . . . . . . . . . .
1.4.3 Method of Lines Transpose . . . . . . . . . . . . . . . . . . . . . . . .
1.4.4 Numerical Methods for non-ideal MHD . . . . . . . . . . . . . . . . .

Chapter 2 The Ideal MHD Equations . . . . . . . . . . . . . . . . . . . . . .
2.1 The Ideal MHD Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Hyperbolicity of the governing equations . . . . . . . . . . . . . . . .
2.2 Overview of the Base Scheme WENO-HCL . . . . . . . . . . . . . . . . . . .
2.3 Magnetic potential in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Magnetic potential in 2D . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 The advantage of MOLT based current method over WENO-HJ method . .
2.5 Outline of the constrained transport methodology . . . . . . . . . . . . . . .

Chapter 3 Kernel Based Approximation with Constrained Transport for
Ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Hamilton-Jacobi equations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1
1D Hamilton-Jacobi equations . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Approximation of the first order derivative ∂x . . . . . . . . . . . . .
3.1.2.1 Periodic boundary conditions . . . . . . . . . . . . . . . . .
3.1.2.2 Non periodic boundary conditions
. . . . . . . . . . . . . .
3.1.2.3
Space Discretization and WENO-based Quadrature . . . . .
2D magnetic potential equation . . . . . . . . . . . . . . . . . . . . .
3D magnetic potential equation . . . . . . . . . . . . . . . . . . . . .
3.2 Central finite difference discretization of ∇ × A . . . . . . . . . . . . . . . .
3.2.1 Curl in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Curl in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Positivity Limiter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Smooth vortex test in MHD . . . . . . . . . . . . . . . . . . . . . . .
3.4.1
3.4.2
2D Orszag-Tang Vortex . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Cloud Shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.1.3
3.1.4

1
1
3
4
12
12
13
17
18

21
21
22
23
27
30
31
33

36
37
37
38
41
44
46
51
54
55
56
56
57
61
61
62
64

v

2D Blast wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.4
2D Field Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.5
3D Field Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.6
3.4.7
3D Blast wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.8 Kernel-based method with Diffusion Terms . . . . . . . . . . . . . . .

Chapter 4 MOLT for PEC for ideal MHD . . . . . . . . . . . . . . . . . . .
4.1 MOLT for PEC for ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Derivation of the Boundary Conditions for MOLT for Dirichlet BC .
4.2 Hartmann Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Numerical Results for ideal MHD with PEC Boundary Conditions . .

Chapter 5 Non-Ideal MHD Equations . . . . . . . . . . . . . . . . . . . . . .
5.1 NON- IDEAL MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Magnetic Vector Potential Formulation for Resistive MHD . . . . . . . . . .
5.3 Outline for Constrained Transport Framework for Non-ideal MHD . . . . . .

65
68
68
70
71

76
76
78
81
84

87
87
89
93

Chapter 6 Kernel Based Approximation with Constrained Transport for
Resistive MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1 Hamilton Jacobi with Diffusion terms . . . . . . . . . . . . . . . . . . . . . .
6.2 Approximation of the second order derivative ∂xx . . . . . . . . . . . . . . .

96
97
97
6.2.0.1 Periodic boundary conditions . . . . . . . . . . . . . . . . . 101
6.3 Numerical Results for Resistive MHD for Periodic Boundary Condition . . . 101
6.3.1 Refinement Study for Smooth Vortex Problem . . . . . . . . . . . . . 102
6.3.2
Smooth Vortex Problem for Resistive MHD . . . . . . . . . . . . . . 103
6.4 Orszang Tang Vortex problem for Resistive MHD . . . . . . . . . . . . . . . 106
6.4.1 Refinement Study for Orszag Tang Vortex problem for Resistive MHD 106
6.4.2 Orszag Tang Vortex problem with diffusion for Resistive MHD . . . . 107

Chapter 7 MOLT for PEC for non-ideal MHD Equations . . . . . . . . . . 110
7.1 MOLT for PEC for non-ideal MHD . . . . . . . . . . . . . . . . . . . . . . . 110
7.2 Numerical results for Hartmann Flow problem with diffusion . . . . . . . . . 111
7.3 Magnetic reconnection: GEM challenge problem . . . . . . . . . . . . . . . . 112

Chapter 8 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . 119
8.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

vi

LIST OF TABLES

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

. . . . . . . . . . . . . . . . . . .

43

Table 3.2: Smooth vortex problem. Errors of B and orders of accuracy.

. . . . . . .

63

Table 4.1: Hartmann Flow for ideal MHD. Errors of B and orders of accuracy.

. . .

84

Table 6.1: Smooth vortex problem. Errors of A and orders of accuracy.

. . . . . . . 103

Table 6.2: Errors of A and orders of accuracy.

. . . . . . . . . . . . . . . . . . . . . 107

Table 6.3: Errors of A and orders of temporal accuracy.

. . . . . . . . . . . . . . . . 107

Table 7.1: Hartmann flow problem. Errors of A and orders of accuracy.

. . . . . . . 112

vii

LIST OF FIGURES

Figure 3.1: The structure of the stencils in WENO integration.
. . . . . . . . . . .
Figure 3.2: Smooth vortex problem. L2 norms of u and B at t = 20 with 160 × 160
. . . . . . . . . . .

grid points. Domain for the plots is [−5, 5] × [−5, 5].

48

64

grid points.

192 × 192 grid points.

Magnetic Pressure (cid:107)B(cid:107); b) Magnetic Potential A3.

Figure 3.3: Orszag-Tang vortex problem. Contour plots of density at t = 3 with
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 3.4: Cloud shock problem. Contour plots of (cid:107)B(cid:107) at t = 0.06 with 512 × 512
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 3.5: 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 (cid:107)B(cid:107).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 3.6: 2D Field loop. Contour plots at time t = 1 with 128× 128 grid points. a)
. . . . . . . . . . . .
Figure 3.7: 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.
Figure 3.8: 3D Blast wave problem. Contour plots at time t = 0.01 with 150 × 150
grid points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 3.9: 3D Blast wave problem. Contour plots at time t = 0.01 with 150 × 150
grid points. While the max (cid:107)u(cid:107) value for previous method is 261, it is
320 for the kernel-based method. Note that this is 23% higher than the
old method, and is far less isotropic around the peak than the old method 73

67

69

65

66

70

72

Figure 3.10: 3D Blast wave from the kernel based method with having added the dif-
fusion limiter. The max (cid:107)u(cid:107) value is 265. We note that these results
are nearly identical to the old WENO-HJ method, which required the
diffusion limiter.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

Figure 4.1: Sketch of a Hartmann type of flow [58].
83
Figure 4.2: 1D cross-sections of Hartmann problem for Ha = 10 at t = 1 with 128×128. 85

. . . . . . . . . . . . . . . . . .

Figure 4.3: Velocity profile for the Hartmann problem for Ha = 10 in 2D at t = 1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

with 128 × 128.

85

viii

Figure 4.4:

Induced magnetic field for the Hartmann problem for Ha = 10 at t = 1
with 128 × 128.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

Figure 4.5: Magnetic potential A for the Hartmann problem at Ha = 10.

. . . . . .

86

Figure 6.1:

Figure 6.2:

Figure 6.3:

Figure 6.4:

(FD) Contour plots of (cid:107)B(cid:107) at time t = 20 with ∆x = ∆y = 0.2.
(Kernel based method) Contour plots of (cid:107)B(cid:107) at time t = 20 with ∆x =
∆y = 0.2.
(FD) Contour plots of ρ at time t = 3 with 192 × 192.
. . . . . . . . . . 108
(MOLT) Contour plots of ρ at time t = 20 with 192 × 192. . . . . . . . . 109

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

. . . . . 104

Figure 7.1: Velocity profile for the Hartmann problem at Ha = 0, 2, 5, 10.

. . . . . . 113

Figure 7.2: Exact and Numerical velocity solutions for the Hartmann problem at

Ha = 0, 2, 5, 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 7.3:

Induced magnetic field for the Hartmann problem at Ha = 0, 2, 5, 10.

. . 115

Figure 7.4: Exact and Numerical induced magnetic field solutions for the Hartmann

problem at Ha = 0, 2, 5, 10.

. . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure 7.5: Magnetic Reconnection in Space. Credits: NASA . . . . . . . . . . . 117

Figure 7.6: Breaking and rejoining of magnetic field lines due to the presence of lo-

calized diffusion region (shaded). ( Priest, E. R [63])

. . . . . . . . . . . 118

Figure 7.7: GEM Magnetic Reconnection Problem for diffusion η = 0.001.

. . . . . . 118

ix

Chapter 1

Introduction

1.1 Outline of the Thesis

In this thesis, our first intention is to develop a hybrid Magnetohydrodynamics (MHD) solver

to solve ideal MHD so that the solutions satisfy the divergence free constraint. Our second

goal is to extend a similar approximation to the non-ideal MHD models, in particular to the

resistive MHD model. In this section we will explain how this thesis is organized.

In Chapter 2, we briefly review the ideal MHD equations and the evolution of the mag-

netic vector equations with respect to magnetic potential equations. We propose the steps of

our version of the constrained transpose framework, which is based on the our novel kernel

based approximation scheme. Our proposed schemes are coupled with base-scheme weighted

essentially non-oscillatory (WENO) with the kernel based scheme for the constrained trans-

port (CT) methodology. Also, we give a short overview of WENO scheme which is designed

to solve Hyperbolic Conservation Laws [48], which we refer as WENO-HCL. Note that this

thesis is not about WENO, however WENO-HCL is a main tool we use in developing our

solver.

In Chapter 3, we represent our novel kernel based numerical scheme for 1D Hamilton-

Jacobi equations and present how we transform the first spatial derivative terms via con-

volution integrals. We develop a novel MOLT-HJ solver to address the Hamilton Jacobi

1

equations as a part of our MHD solver in 2D and 3D as well as the WENO-based quadrature

formulation. The main advantage of the proposed scheme is to be able to avoid the artificial

diffusion terms, which were needed in the previous work [29]. Moreover, the stability analysis

is shown for the linear problems. Several test problems in 2D and 3D are presented in the

numerical results section to verify the order of accuracy and to show the robustness of the

scheme. We compare the new method with the old method to validate that the new method

works well. We observe that the limiter needed in the 3D vector potential in the old method

is causing a change in the structure of the solutions in the context of the strong shocks,

while the new method conserve the solution structure for the problems with strong shocks,

in particular for the blast wave problem.

In Chapter 4, we introduce our proposed scheme with respect to PEC boundary condi-

tions and show it is equivalent to enforce Dirichlet boundary conditions. In fact, we develop

the Dirichlet boundary conditions for magnetic potential in MOLT-HJ. Then, we go over

the derivation of an important benchmark problem, Hartmann Flow problem [58]. Hart-

mann Flow problem is a well known test problem to check PEC boundary conditions. Our

proposed kernel based scheme in the context of PEC boundary conditions are tested via

the Hartmann Flow problem. We demonstrate that we could get high order accuracy and

compare the numerical solutions with the analytical solutions.

In Chapter 5, we extend the MHD model to include diffusion terms in the magnetic field

equation. Then we derive a set of magnetic potential equations for the new case of having

diffusion terms in the magnetic field equation. Later, the constrained transport strategy for

the case of resistive MHD is presented.

In Chapter 6, we extend the MOLT-HJ solver to include magnetic diffusion terms. We

present kernel based representations for the second order spatial derivative terms in the

2

magnetic potential equation. Such representations are based on a successive convolution.

Then we test the scheme with 2D benchmark problems to test our proposed method with

respect to periodic boundary conditions.

In Chapter 7, we then present kernel based method for PEC boundary condition in

the case of resistive MHD. We then work on several test problems such as Hartmann Flow

and GEM Magnetic Reconnection problem. And in Chapter 8, we conclude with a brief

discussion about the conclusions.

In the next section, we will explain the plasma phenomenon and then present some

important models, which describe the behavior of plasma. We especially derive the equations

for the MHD model. In the rest of this chapter we will present a literature review.

1.2 Plasma Phenomenon

One of the most common but least understood phenomenon in the universe is plasma. A

plasma is a collection of charged particles, containing many interacting free electrons and

ionized atoms or molecules which are obtained by heating a gas until it breaks down into

ints constitutive elements. A plasma is the fourth state of the matter: solid, liquid, gas, and

plasma. The primary differences among solids, liquids, gases and plasmas are based on the

molecular relationships between the particles in the states. The strength of the bonds that

hold their constituent particles together is strong in a solid state, moderate in liquid state,

and weak in the gaseous state. The thermal energy (random kinetic energy) level of the

particles and the strength of the binding forces determines the state of a given substance.

Applying a high enough constant temperature to a solid or liquid substance generates thermal

kinetic energy to the molecules such that they exceed the binding potential energy and this

3

leads to red a phase transition. However, in the sense of thermodynamics, the transition

from a gas to a plasma is not considered as a phase transition since a plasma occurs by

heating a substance with gradually increasing temperature until a reasonable high fractional

ionization is produced.

The particle interactions determine the properties of a plasma. For example, a basic

feature of a plasma is that charged particles simultaneously interact each other over long

distances, which is a distinct behavior of a plasma from ordinary fluids and solids. Since the

particles in a plasma are interacting with electromagnetic fields, this provides many novel

phenomena that do not exist in fluids or solids. A significant feature of a plasma is the

ability to have a various wave phenomena. Therefore, the studies of waves in plasma present

important information on plasma properties.

The plasma state commonly exists in the Universe. For example, common instances

of matter in the plasma state are stars and galaxies. There are a range of man made

plasma applications including thermonuclear fusion reactions, which is a possible source of

energy; the magnetohydrodynamic generator, which converts kinetic energy of a plasma into

electrical energy (MHD); and plasma propulsion systems for space propulsion.

1.3 Theoretical Overview of Plasma

The plasma dynamics is defined by the interaction between the plasma charged particles

with internal and external fields. The motion of the charged particles can generate electric

fields and magnetic fields. The dynamic behavior of the particles in plasma is described by

the classical (non-quantum) mechanics laws. The Lorentz F orce gives the interaction of

charged particles with electromagnetic fields. For a particle with charge q, the motion of the

4

equation is

= q(E + v × B),

dp
dt

(1.3.1)

where m is mass, v is the velocity of the moving particle, E is the electric field and B is the

magnetic field, and p = mv denotes the particle momentum. The dynamics of a plasma are

described by solving the equation of motion for the each charged particle under the combined

effect of the internally generated fields by all the other particles and external fields. For N

number particles, N nonlinear coupled motion equations needs to be solved simultaneously.

The electromagnetic fields are governed by the M axwell equations

∇ × E = −∂B
∂t

,

∇ × B = µ0(J + 
0

∂E
∂t

),

,

∇ · E =

ρ

0
∇ · B = 0,

(1.3.2)

(1.3.3)

(1.3.4)

(1.3.5)

where ρ is the total charge density, J is the total electric current density, 
0 is the electric

permittivity, and µ0 is the magnetic permeability of free space.

There are several different choices of approximations for the theoretical description of

plasma phenomenon. Each of these models applies to different situations of plasma.

The first choice is to use a kinetic model of the plasma which is based on the distribution

function for the system of particles at each point. The kinetic equations that describe the

evolution of the distribution function in phase space are 6D plus time. One example of

differential kinetic equations to describe plasma is the Boltzmann equation. It governs the

5

evolution of a probability density function fα(t, x, v) of plasma species α in phase space (x, v)

at a given time t. The function fα(t, x, v) is also called the distribution function in phase

space and represents the possibility of density of representative points of particles in phase

space.

Another approach is two-fluid modeling, which is based on a local distribution function

for each species. If there exists a very frequent collisions between the plasma particles, then

each species (electrons and ions) can be considered as a fluid described by a local variables

such as density, velocity and temperature. In this two fluid model, there exist conservation of

mass, conservation of momentum, and conservation of energy for each species in the plasma.

A third model is magnetohydrodynamic (MHD), which is single fluid model, which treats

two-fluid models as a single conducting fluid using unified macroscopic variables and their

related conservation equations. Thus, the single fluid model is a set of coupled equations of

two-fluid model equations, which is simpler to study. Single fluid model is less complete, but

it is still gives important information about the behavior of a plasma.

The MHD model is applicable to the studies of highly conducting fluid plasmas engaged

with magnetic fields. It’s a model to understand the solar and stellar winds and transport

of plasmas inside magnetospheres.

The derivation of fluid models from the Boltzmann equation requires the introduction

of a ”closure” model. This ”closure” model is a key weakness in a fluid model because

it introduces assumptions about the behavior of the underlying micro scale physics, but is

necessary to arrive at a finite set of moments. Now we briefly show the general strategy to

derive the fluid equations from kinetic equations, [12, 60].

6

=

qα
mα

(E + v × B),

aα =

Fα
mα

Ä ∂fα
ä
Ä ∂fα

∂t

ä

where mα is the mass of species α and qα is the charge of α. The term at the right hand

side of the Boltzmann equation

coll describes the changes of plasma due to particle

collisions and is called a general collision operator. Also note that in the case where the

collisions in particles are neglected,

∂t

coll = 0, the Boltzmann equation (1.3.6) is called

the Vlasov equation.

The Boltzmann equation is multiplied by a weight function χ(v) and then the resulting

equation is integrated with respect to the velocity space

Ç∂fα

∂t

(cid:90)

χ

å

(cid:90)

χ

dv =

Ç∂fα

å

Mathematically, the Boltzmann equation is as follows

∂fα
∂t

+ v · ∇fα + aα · ∇vfα =

with the acceleration

Ç∂fα

∂t

å

,

coll

(1.3.6)

(1.3.7)

+ v · ∇fα + aα · ∇vfα

dv.

(1.3.8)

∂t

coll

These weight functions χ(v) represent some physical features of particles in plasma. By

properly choosing the function χ(v), we end up with continuity equation, the momentum

equation, and the energy equation. More precisely, if χ(v) = mα, we get the continuity

equation, if χ(v) = mαv, we get the momentum equation, and if χ(v) = 1

2 mαv2, we

get the energy equation. Consequently, there is a system of transport equations for each

α species since a plasma consists of multiple particle species such as electrons, ions, and

7

neutral particles. If we denote < χ >α as the average over v-space as

(cid:90)

χfαdv,

< χ >α=

then we obtain the system of transport equations as follows

Ç∂ρα

å

(1.3.9)

(1.3.10)

(1.3.11)

dv,

(1.3.12)

,

v

Ç∂fα
å
Ç∂fα

∂t

å

(cid:90)
(cid:90)

dv,

coll

∂t

coll

î ∂uα

∂t

ρα

+ (uα · ∇)uα

ó

∂ρα
∂t

+ ∇ · (ραuα) =

∂t

coll

+ ∇ · (Pα)− < Fα >α= mα

+ ∇ · (
αuα) + ∇ · (Pα · uα) − nαqαuα · E =

∂
α
∂t

1
2

mα

||v||2

where the right hand side of the above equations are collision terms. The fluid unknown

quantities are the mass density ρα = nαmα = mα < 1 >α, the momentum density ραuα =

ρα < vα >α, and the energy density


α =

pα
γ − 1

+

1
2

ρα||uα||2,

(1.3.13)

where γ is the ideal gas constant, and Pα = mα < (v − uα) ⊗ (v − uα) >α is the pressure

tensor. If there exist Lorentz force, then the term < Fα >α becomes

< Fα >α= qα(E + uα × B),

(1.3.14)

where E is electric field and B is magnetic field.

In the case of collisonless plasma, the

8

conservation equations for each species α reduces to

+ ∇ · (ραuα) = 0,

î ∂uα

∂t

ρα

∂ρα
∂t

+ (uα · ∇)uα

ó

= nαqα(E + uα × B) − ∇ · (Pα),

+ ∇ · (
αuα) + ∇ · (Pα · uα) − nαqαuα · E = 0.

∂
α
∂t

(1.3.15)

(1.3.16)

(1.3.17)

What we come up with is a fluid model that treats each species (ions and electrons) in a

plasma as a different fluid, [12]. Therefore we end up with a model for two species so called

two fluid model where each individual plasma species lead to a different system of transport

equations. However, there is a simpler description of plasma dynamics as the evolution of

one single fluid, where the plasma is considered a conducting fluid. It is a system of transport

equations which treats a plasma as a whole by combining each species variables. In other

words, the quantities that are relative to each species, we can actually deduce to the one-fluid

quantities that are evolved by the single-fluid model. For example, some of the unknown

variables such as the mass density ρ(t, x), the momentum density ρu(t, x), and the electric

current density J are given, respectively, by the sum of the electron contribution and the ion

contribution as follows

ρ(t, x) =

nαmα,

ρu(t, x) =

nαmαuα,

J =

nαqαuα.

α

(cid:88)
(cid:88)
(cid:88)

α

α

9

(1.3.18)

(1.3.19)

(1.3.20)

Then the conservation equations for the plasma as a whole can be given as

+ ∇ · (ρu) = 0,

î ∂u

∂t

ρ

∂ρ
∂t

+ (u · ∇)u

ó

= ρcE + J × B − ∇ · (Pα),

+ ∇ · (
u) + ∇ · (P · u) − ρcu · E = 0,

∂

∂t

where ρc is charge density so that ρc =(cid:80) nαqα.

(1.3.21)

(1.3.22)

(1.3.23)

Furthermore, for the type of collisionless case of the transport equations with a closure

condition on the pressure tensor Pα and some additional assumptions, the system leads to

the three conservation laws, which become the Magnetohydrodynamics (MHD) model, when

an additional approximation is introduced for E. The fourth equation of the MHD system

is the magnetic field equation comes from the ”modified” Maxwell equations,

∇ × E = −∂B
∂t

,

∇ · B = 0,

∇ × H = J,

∇ × D = 0,

(1.3.24)

(1.3.25)

(1.3.26)

(1.3.27)

where B = µ0H, D = 
0E and the electric field E is obtained by the Ohm’s law, [12, 60].

The generalized Ohm’s law is

u × B

c

E +

= ηJ +

J × B
enec

− ∇ · Pe

neec

+

me
nee2

dJ
dt

.

(1.3.28)

10

There are multiple versions of the MHD models. Depending on which term we keep at

the right hand side of the generalized Ohm’s law, we end up with a specific MHD model.

For example, keeping the all terms at the right hand side as zero is called the ideal Ohm’s

law

E + u × B = 0,

(1.3.29)

and the resulting system becomes the ideal MHD model by decoupling the system with

Maxwell’s equations by doing the following assumptions

the current density:

J = e(niui − neue),

where ni = ne = n since the plasma assumed to be quasi-neutral because of Gauss’s law,

the charge density:

ρc = e(ni − ne) = 0.

So, the charge density is ignored by the MHD formulation. Also, another assumption is

to neglect the electromagnetic waves which in practice means neglecting the displacement

current 
0

∂E
∂t in the Maxwell’s equations.

Moreover, keeping just the ηJ term in the Ohm’s law gives the resistive MHD. Ohm’s

law then becomes

J = σ(E + u × B).

Using this ohm’s law in the Faraday’s equation (1.3.29) to eliminate E

∇ ×Ä J

σ

− u × B

= −∂B
∂t

.

ä

11

(1.3.30)

(1.3.31)

Then using ∇ × H = J and B = µ0H relations, we get

∇ ×Ä∇ × B

ä − ∇ × (u × B) = −∂B

µ0σ

.

(1.3.32)

∂t

Then applying the vector identity ∇ × (∇ × B) = −∇2B + ∇(∇ · B) and ∇ · B = 0, we

obtain the magnetic field evolution equation

∂B
∂t

=

1
µ0σ

∇2B + ∇ × (u × B),

(1.3.33)

where is σ conductivity and 1

µ0σ is the magnetic viscosity. Also,

1

µ0σ∇2B describes the

effects of diffusion in the resistive MHD model.

Another MHD model is the Hall MHD model, which is derived by keeping ηJ + J×B
enec

terms in the Ohm’s law. In this work, our focus is on both the ideal MHD model and the

resistive MHD model and we intend to work on the Hall MHD system in the future. For

more details, see [12, 60].

1.4 Review of Previous Work

1.4.1 Numerical Methods for Hyperbolic Conversation Laws

As stated, the ideal and resistive MHD models are mathematically hyperbolic conservation

laws. This is the base system of the MHD model before applying the constrained transport

correction methodology. There have been many numerical methods applied to MHD system

to solve it in the literature. For example, the discontinuous Galerkin method [30], the finite

volume / finite difference essentially non-oscillatory (ENO) methods [41, 68, 69], and finite

12

volume / finite difference Weighted essentially non-oscillatory (WENO) methods [45, 53].

The WENO schemes are the most robust and efficient schemes for the problems that contains

discontinuities comparing to other methods. Thus, we use the WENO scheme in this thesis

as a base-scheme for the ideal and resistive MHD equations. In section (2.2), we will provide

the WENO scheme described for Hyperbolic Conservation Laws [48], which we will refer to

this method as the WENO-HCL scheme.

For the conservation laws type equations, the non-oscillatory property is required for

the time discretization. Therefore, We use strong stability preserving (SSP) Runge Kutta

methods for temporal evolution.

In particular, we use the optimal third order SSP-RK

scheme (SSP-RK3) in our works which can be given by [68]. For the following equation

∂tq = L(q),

(1.4.1)

the SSP-RK3 can be expressed as

q(1) = qn + ∆tL(qn),

q(2) =

qn +

3
4

1
4

q(1) +

1
4

∆tL(q(1)),

qn+1 =

qn +

1
3

2
3

q(2) +

2
3

∆tL(q(2)),

(1.4.2)

where ∆t denotes time step and the supper scripts mean stage of SSP-RK.

1.4.2 Numerical Methods for ideal MHD

The ideal Magnetohydrodynamics (MHD) equations are one of the most important clas-

sical 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

13

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 ap-

proaches overcoming the difficult in the literature: 8-wave formulation [61, 62]; projection

methods [4, 71]; hyperbolic divergence cleaning methods [33]; and constrained transport

methods [2, 8, 31, 37, 38, 42, 54, 55, 64, 65, 32, 71, 70, 29].

The first methodology used for divergence free condition is developed in [14] 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 [79] 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 [61]. 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 [71] 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 [33]

is the hyperbolic divergence-cleaning method. This method is a similar to the projection

14

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 [37]. The original CT method is considered to be a modification

of the Yee method [78] from electromagnetics for the ideal MHD equations. In the origi-

nal 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 [34] 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 [8], Dai and Woodard [31], and Ryu et al. [65]. Londrillo

and Zanna [54, 55] constructed one of the first high order upwind scheme based on the work

of Evans and Hawley. De Sterck [32] 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 [71] compared several of these schemes maintaining divergence free

15

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 [38] developed an unstaggered upwind

method satisfying the divergence free constraint. Rossmanith [64] designed an ustaggered

wave propagation scheme for MHD flows based on the algorithms in [51] using a constrained

transport method to keep the divergence free magnetic field. Helzel et al. [42, 43] 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 equa-

tions 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] in-

troduced a high accuracy ADER-WENO schemes for divergence free magnetohydrodynamics

on structured meshes, again using a staggered magnetic field. Li et al. [52, 39] and Cheng

et al. [20] presented high order central discontinuous Galerkin schemes satisfying the diver-

gence free constraint globally for ideal MHD simulations on two overlapping meshes, called

primal and dual meshes, by utilizing different discretization for magnetic induction equa-

tions. Kawai [49] 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 [29], 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

16

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 appear in HJ equation.

x and A+

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.

1.4.3 Method of Lines Transpose

In this current work, we further our work [29] 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

[22, 23], which is derived from the Method of Lines Transpose (MOLT ) [19, 16, 18, 17, 21, 15].

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 [29] since the kernel based

method is a fully implicit method. The corrected magnetic field is computed by 4th-order

17

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 [47] as the base scheme for ideal MHD equations. Third order explicit

strong-stability-preserving (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.

1.4.4 Numerical Methods for non-ideal MHD

Magnetic reconnection is a significant phenomenon occurring in our solar system. When

oppositely directed magnetic field lines in plasma come together, a strong current sheet

is established. Even a small resistivity in this region becomes important, allowing plasma

diffusion and thus, magnetic reconnection happens. In an ideal MHD model, there is a frozen-

in flux condition, therefore there is no reconnection happening. However, Biskamp showed

in his work [10] that for the resistive MHD model, the frozen-in constraint is broken in the

dissipation region and reconnection occurs. The dissipation region establishes Sweet-Parker

layer for a small resistivity value and the rate of reconnection is very low. In [50, 74], the

inertia of electrons breaks the frozen-in constraint and the collision rate is very small in the

magnetosphere. In the work about hybrid simulation with resistivity [57], in the work about

two fluid simulations [56, 11] and in the work about the particle simulations [67, 44] it is

18

found that the rate of reconnection is insensitive to the strength of dissipation. Furthermore,

the narrow current layers remain microscopic in length opposed to macroscopic Sweet-Parker

current sheet. In [11], it is shown that the macroscopic current reappears, when the Hall

term is removed from the simulations, thus the reconnection rate is decreased and it becomes

sensitive to the electron inertia. Therefore, the Hall effect in the MHD model is an important

term in determining collisionless magnetic reconnection rates. In [9] they investigated the

magnetic reconnection in 2D in a simple Harris sheet configuration for resistive MHD and

hall MHD. Our goal is to test our kernel based approximation on some magnetic reconnection

benchmark problems.

Moreover, the choice of time stepping is also important for solving the non-ideal MHD

equations, which have diffusion terms. For choice of explicit time stepping, there is very strict
constraint for ∆t due to the stability condition (i.e., ∆t ≤ (∆x)2

), where b is the diffusivity

2b

parameter. This stability condition of explicit schemes gives the size of the time step and

forces ∆t to be very small, which leading to many time steps and taking quite computational

time in the simulations. However, there are other ways of solving diffusion equations like

using implicit solver for time stepping. Implicit schemes reduce the stability condition and

allow larger time steps, but then there is a resulting a system of linear equations needed to

be solved which may also take time. In the end, using implicit schemes may still save some

computational time. In [72] Udrea develops an implicit algorithm for solution of the resistive

MHD. The implicit formulation employed by their solver (WARP3) is a two level scheme

and they use symmetric Gauss-Seidel method to solve the linear system resulting from the

implicit scheme. Their solver addresses to solve steady state problems and their code also

has an explicit mode capability. They tested that for time dependent simulations the implicit

scheme in their solver WARP3 is faster than the explicit scheme. In a similar work, in [36]

19

Dubuc used a generalized conjugate gradient scheme to solve the linear system that results

from the implicit discretization of the Euler equations. For the linear system resulting from

the implicit discretization there are other direct methods (such as Gaussian elimination),

which requires larger computer memory and thus not desired. Vanden and Whitfield [73]

investigated the use of a direct method for a linear system that results from the discretization

of Euler equations and they showed that their implementation of the direct method required

approximately 3.5 more memory and 14.5 more CPU time than a Gauss-Seidel iterative

scheme.

In this work, we choose explicit time stepping to solve resistive MHD and aim to develop

an unconditionally stable kernel based approximation to solve Hamilton Jacobi equations

with diffusion terms, which will allow us to take larger time steps. We will test our solver

with our previous work [29], which is originally designed to solve ideal MHD. However, we

will add second order approximation for diffusion terms in previous finite difference weighted

essentially non-oscillatory (FD-WENO) solver.

20

Chapter 2

The Ideal MHD Equations

2.1 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

∂t

ρ





ρu
E

B



ρu

+ ∇ ·

ρu ⊗ u + ptotI − B ⊗ B
u(E + ptot) − B(u · B)

u ⊗ B − B ⊗ u

∇ · B = 0,



= 0,

(2.1.1)

(2.1.2)

(2.1.3)

with the equation of state as

E =

p

γ − 1

+

ρ(cid:107)u(cid:107)2

2

(cid:107)B(cid:107)2
2

.

+

where the total mass ρ, the momentum ρu = (ρu1, ρu2, ρu3)T , the energy densities E of the

plasma system, and the magnetic field vector B = (B1, B2, B3)T are all conserved variables.
2(cid:107)B(cid:107)2, together with the hydrodynamic

The velocity u and the total pressure ptot = p + 1

21

pressure, which is given by the ideal gas law as

p = (γ − 1)(E − 1
2

(cid:107)B(cid:107)2 − 1
2

ρ(cid:107)u(cid:107)2),

(2.1.4)

are derived quantities. γ = 5/3 is the ideal gas constant and the notation (cid:107) · (cid:107) is used for

the purpose of the Euclidean vector norm. The MHD equations (2.1.1)-(5.1.2) are derived

and discussed in many standard plasma textbooks (e.g., [59]).

2.1.1 Hyperbolicity of the governing equations

The system (2.1.1) along with the equation of state (5.1) comprise a system of hyperbolic

conservation laws

qt + ∇ · F(q) = 0.

(2.1.5)

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 ((cid:107)n(cid:107) = 1), they can be written as

λ1,8 = u · n ∓ cf
λ2,7 = u · n ∓ ca
λ3,6 = u · n ∓ cs
λ4 = u · n
λ5 = u · n

fast magnetosonic waves,

(2.1.6a)

Alv´en waves,

(2.1.6b)

slow magnetosonic waves,

(2.1.6c)

entropy waves,

(2.1.6d)

divergence waves,

(2.1.6e)

22

where

a ≡

ρ
(B · n)2

  γp
Ã
1
a2 +
1
a2 +

ρ

2

2

ca ≡

cf ≡

cs ≡

ÃÇ
ÃÇ

a2 +

a2 +

å2 − 4a2 (B · n)2
å2 − 4a2 (B · n)2

ρ

ρ

(2.1.7a)

(2.1.7b)

(2.1.7c)

(2.1.7d)

1/2



1/2

(cid:107)B(cid:107)2
ρ
(cid:107)B(cid:107)2
ρ

(cid:107)B(cid:107)2
ρ
(cid:107)B(cid:107)2
ρ

+

−

The eight eigenvalues are well ordered as

λ1 ≤ λ2 ≤ λ3 ≤ λ4 ≤ λ5 ≤ λ6 ≤ λ7 ≤ λ8

(2.1.8)

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

linearly degenerate.

2.2 Overview of the Base Scheme WENO-HCL

Assuming that B is always divergence free, we can simply apply the standard WENO

methodology. Here we review this methodology before talking about of form of constrained

transport. We present the 5th-order Weighted Essentially Non-Oscillatory method developed

for Hyperbolic Conservation Laws (WENO-HCL) [48] for the MHD system that comprises

the base scheme in the CT procedure. WENO-HCL method is described by Jiang and Wu

23

in [48]. We will briefly explain the scheme for a one dimensional MHD model

∂q
∂t

+

∂f (q)

∂x

= 0

q = (ρ, ρu1, ρu2, ρu3,E, B1, B2, B3)T ,

ρu1, ρu1u1 + p +

1
2

||B||2 − B1B1, ρu1u2 − B1B2, ρu1u3 − B1B3,

||B||2ä − B1(u · B), 0, u1B2 − u2B1, u1B3 − u3B1

åT

where

f (q) =

Ç
u1ÄE + p +

1
2

(2.2.1)

(2.2.2)

(2.2.3)

Also, we denote the primitive variables as

w = (ρ, u1, u2, u3, p, B1, B2, B3)T

(2.2.4)

The Jacobian of the flux ∂f

∂q of the above equation can be decomposed based on eigen-

vectors since the MHD system is hyperbolic

∂f
∂q

= RΛL

(2.2.5)

where Λ is the diagonal matrix of real eigenvalues, R is the right eigenvectors , and L = R−1

is the left eigenvectors.

We discretize the space into uniform intervals with ∆x size. Denoting qi(t) as the nu-

merical solution of the MHD model at the point x = xi, the WENO-HCL method for the

24

MHD system in the conservative form

dqi(t)

dt

=

1
∆x

Ç

ˆF

− ˆF

i− 1
2

i+ 1
2

å

(2.2.6)

The following WENO procedure is applied to get the numerical flux ˆF

for the above semi

i+ 1
2

discrete form:

1. Compute the physical flux and the solution qi for all i

fi = f (qi).

(2.2.7)

2. Do the following steps at each fixed x

i+ 1
2

(a) Compute the average state w

using Roe average

i+ 1
2

=

w

i+ 1
2

1
2

(wi + wi+1)

(2.2.8)

(b) Compute the eigenvalues Λ

of the Jacobian of flux, ∂f

i+ 1
2

Also, compute the right eigenvector R

and left eigenvector L

Ç

å

i+ 1
2

Ç

å

∂q , at the point x

.

i+ 1
2

Ç

i+ 1
2

å

Λ

i+ 1
2

= Λ

w

,

i+ 1
2

R

i+ 1
2

= R

w

,

i+ 1
2

L

i+ 1
2

= L

w

i+ 1
2

.

(2.2.9)

(c) Project the solution and the physical flux to local characteristic field (right eigen-

vector space)

Vj = L

i+ 1
2

qj

and Gj = L

fj

i+ 1
2

(2.2.10)

where j = i − 2, i − 1, i, i + 1, i + 2, i + 3 for the 5th-order WENO.

25

(d) Perform the Lax Friedrichs flux splitting procedure for the each component of the

characteristic variables, to obtain the corresponding component of the Vj and

Gj. We can assume the m-th component of the characteristic variables of the Vj

and Gj as vj and gj, respectively and compute

ä

g±
j =

1
2

gj ± α(m)vj

,

(2.2.11)

where the maximum wave speed of the m-th component of characteristic variables

is

α(m) = max

k

|λ(m)(qk)|

(2.2.12)

(e) Compute the numerical flux by applying a WENO reconstruction on the computed

flux components g±

j ’ Let’s WENO5 denotes the 5-th order WENO reconstruction,

then the numerical flux is obtained as

å
å

,

.

ˆg+
i+ 1
2

ˆg−
i+ 1
2

= W EN O5

g+
i−2, g+

i−1, g+

i , g+

i+1, g+
i+2

= W EN O5

g−
i+3, g−

i+2, g−

i+1, g−

i , g−
i−1

Ä

Ç
Ç

Then we compute the m-th component of the numerical flux ˆG

i+ 1
2

as

ˆg
i+ 1
2

= ˆg+

i+ 1
2

+ ˆg−
i+ 1
2

(f) Project the numerical flux back into physical space

ˆF

i+ 1
2

= ˆR

ˆG

i+ 1
2

i+ 1
2

26

(2.2.13)

(2.2.14)

(2.2.15)

(2.2.16)

The scheme described above is the WENO-HCL scheme that serves as the base scheme in

our Constrained Transport framework. Also, the high order SSP-RK methods are used for

the time discretization for the MHD model.

As mentioned earlier, the direct application of the WENO-HCL scheme to the MHD

model leads to divergence error in the magnetic field, and we update the magnetic field

solution through magnetic potential equation by the CT method. In the next section, we

will derive the magnetic potential equation from the magnetic field equation for the ideal

MHD model.

2.3 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 the magnetic vector potential equation

from the magnetic field equation given in (2.1.1)-(5.1.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 = ∇ × A.

(2.3.1)

We show in the section (3.2) that if we construct the curl using 4th order central differences,

we locally preserve the divergence free property. The key step of the constrained transport

scheme is to solve the magnetic potential for correcting the magnetic field. The evolution of

27

the magnetic field in (2.1.1) can be written in the following form

Bt + ∇ × (B × u) = 0,

(2.3.2)

using the relation

∇ · (u ⊗ B − B ⊗ u) = ∇ × (B × u).

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

∇ ׶

©

At + (∇ × A) × u

= 0.

This implies that there exists a scalar function ψ such that

At + (∇ × A) × u = −∇ψ.

There are various choices of the gauge conditions depending on how we chose the ψ. Helzel
[42] showed that using the Weyl gauge, i.e., setting ψ ≡ 0, one can achieve stable

et al.

solutions. This condition results in the evolution equation for the magnetic vector potential

as

which can be written as:

∂tA + (∇ × A) × u = 0,

(2.3.3)

∂tA + N1Ax + N2Ay + N3Az = 0,

(2.3.4)

28

with

N1 =



0

0

0 −u2 −u3

u1

0

0

u1

 , N2 =



0

u2
0
−u1 0 −u3

0

0

u2

 , N3 =



u3

0

0

u3

0
0
−u1 −u2 0

 .

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 = (n1, n2, n3) to show the weakly hyperbolicity:



n2u2 + n3u3

−n2u1
−n3u1

−n1u2

n1u1 + n3u3

−n1u3
−n2u3

−n3u2

n1u1 + n2u2

 , (2.3.5)

M (n, u) := n1N1 + n2N2 + n3N3 =

which has real eigenvalues for all (cid:107)n(cid:107) = 1 as

λ1 = 0,

λ2 = λ3 = n · u,

(2.3.6)

and the right eigenvectors matrix is

r(1)

(cid:12)(cid:12)(cid:12)(cid:12)r(2)

(cid:12)(cid:12)(cid:12)(cid:12)r(3)

 =

R =



n1 n2u3 − n3u2 u1(u · n) − n1(cid:107)n(cid:107)2
n2 n3u1 − n1u3 u2(u · n) − n2(cid:107)n(cid:107)2
n3 n1u2 − n2u1 u3(u · n) − n3(cid:107)n(cid:107)2

 .

(2.3.7)

If we assume that (cid:107)u(cid:107) (cid:54)= 0 and (cid:107)n(cid:107) = 1, then the determinant of the matrix of the right

29

eigenvector R can be found as

det(R) = −(cid:107)u(cid:107)3 cos(α) sin2(α)

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.3.4) is only weak hyperbolic.

2.3.1 Magnetic potential in 2D

In the 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

∇ · B = B1

x + B2
y .

There is no point to define B3 since any value of B3 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, B1 and B2, can be found as

B1 = A3
y

and B2 = −A3
x.

(2.3.8)

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

A3
t + u1A3

x + u2A3

y = 0

(2.3.9)

30

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

2.4 The advantage of MOLT based current method

over WENO-HJ method

As we will explain in the next chapter about the Constrained Transport Methodology, the

magnetic vector field is updated so that it satisfies the divergence free condition by solving

the magnetic vector potential equations (2.3.9) in 2D and (2.3.4) in 3D. Although the CT

method that we use is the same idea that we used in our previous work [29], in this section

we will explain the challenges of the previous work and present the advantages of the current

work.

The magnetic potential equation is a Hamilton-Jacobi equation since the components of

the velocity are given from the previous time step such that

A3
t + H(x, y, A3

x + u2A3

y) = 0

with Hamiltonian

H(x, y, A3

x + u2A3

y) = u1A3

x + u2A3
y

(2.4.1)

(2.4.2)

In the previous work [29], the WENO-HJ scheme is directly applied to two-dimensional

magnetic scalar potential equation (2.3.9). However, the scheme cannot be directly applied

to the three-dimensional magnetic vector potential equation due the problem being weakly

hyperbolic of the 3D system (2.3.4). Because of the problem is weakly hyperbolic, we ap-

plied a finite difference weighted essentially non-oscillatory (WENO-FD) discretization using

31

arithmetic averages to define derivatives at grid points in the old work [29]. But this type

of discretization approach leads a lack of insufficient numerical resistivity in order to control

unphysical oscillations in the 3D magnetic field equation. Therefore, artificial resistivity

terms into the magnetic vector potential equation were introduced for solving 3D system.

For brevity, we give only the equation of the first component of the magnetic vector potential

system with artificial resistivity term:

t − u2A2
A1

x − u3A3

z = ε1 ∂2A1
∂x2

x + u2A1

y + u3A1

(2.4.3)

The above equation is the equation of the first component of magnetic vector potential

A1, but this equation also contains other two components A2 and A3. With explicit time

stepping, it tends to be unstable. Therefore, we added the the term that is on the right

hand side of the (2.4.3) to stabilize the equation. The additional term on the right hand

side of the equation (2.4.3) is an artificial resistivity in the missing direction, which in the

above equation for Ax, it has an artificial resistivity term in the x-direction. The artificial

resistivity is defined as

ε1 = 2νγ1 ∆x2
δ + ∆t

,

(2.4.4)

where δ is a very small parameter to ensure that ε1 remains bounded as ∆t → 0+, γ1 is the

smoothness indicator of Ax, and ν is a constant to control the magnitude of the artificial

resistivity. Note that the limiter is off in smooth regions.

In the current work, we will circumvent the stabilization challenges by introducing a new
O(N) implicit method that solves magnetic potential equation with an unconditionally stable

solver using the Method of Lines Transpose idea. We no longer need artificial resistivity terms

for solving the magnetic vector potential equations in this novel method. The rest of the

32

scheme is still using the same WENO architecture that we’ve used in the previous work. We

are targeting a better solver that handle things like strong shocks since we do not need to

introduce a limiter for A. We will show the examples of this in the results section.

In the next section we will present the overview of the Method of Lines Transpose which

have been presented in the previous work [23].

2.5 Outline of the constrained transport methodology

The major challenge when dealing with the numerical solution of the system (2.1.1)-(5.1) is

to satisfy the divergence-free condition (5.1.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 [75], as well as Dorfi [35]. However, modern shock capturing strategies were not

used in those works and hence caused strong numerical diffusion. Londrillo and Del Zanna

[54] used the magnetic potential solutions in the context of shock capturing methods, along

with De Sterck [32], Londrillo and Del Zanna [55], and Rossmanith [64]. Helzet et al. [43]

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 [29], 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

33

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.1)

q(cid:48)
mhd(t) = L(qmhd(t))

(2.5.1)

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

And also consider the update for the magnetic potential equation ((2.3.9) for 2D case

and (2.3.4) for 3D case) on the mesh,

q(cid:48)
A(t) = H(qA(t), u(t)),

(2.5.2)

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

The key steps advancing the solution from its current time step t = tn (or the initial

condition at t0) to its new time step tn+1 are listed below:

• step 0: Start with the given current time step qn

mhd = (ρn, ρun,E n, Bn)T and qn
A.

• step 1: Obtain q∗

mhd and qn+1

A

q∗
mhd =

separately, where

ρn+1, ρn+1un+1,E∗, B∗ä
Ä

.

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

34

method before the end of the time step.

• step 2: Replace B∗ to Bn+1 by a discrete curl of qn+1
A .

Bn+1 = ∇ × qn+1
A .

• 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

E n+1 = E∗.

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

E n+1 = E∗ +

(cid:13)(cid:13)(cid:13)Bn+1(cid:13)(cid:13)(cid:13)2 − (cid:107)B∗(cid:107)2)

(

1
2

(Second option sometimes helps to prevent negative pressure).

In Section 3.2, we show that our approach leads to a divergence free solution. We now

describe how we constraint the update for qn+1

A

using our kernel based approach.

35

Chapter 3

Kernel Based Approximation with

Constrained Transport for Ideal MHD

In this chapter we develop a kernel based approximation for magnetic scalar/vector potential

equations, which are considered as Hamilton Jacobi equations since the velocity variables

are given in the previous time step. We also present a WENO quadrature formulation to

approximate the convolution integrals, which appears in our kernel based formulation. We

will refer this new kernel based approximation coupled with WENO quadrature as MOLT-

HJ scheme [23]. This proposed scheme is a key part of the constrained transport method

used to update the solutions such that the magnetic field solution satisfies the divergence

free condition. Our goal is to eliminate the need of artificial terms such as the diffusion

limiter (introduced in the previous work [29]) and to understand what that means for strong

shocks. We develop our solver for both 2D and 3D. We later include the artificial diffusion

terms into the MOLT-HJ scheme in 3D to be able to compare the the new scheme with the

old method. We demonstrate the robustness of the new method on several test problems

in 2D and 3D and validate that the new method works well for the blast wave and strong

shocks by comparing with the old method. Further, we validate the negative impact of the

diffusion limiter by adding the artificial limiters to MOLT-HJ scheme. MOLT with diffusion

limiter obtains very similar results with the previous method.

36

3.1 Hamilton-Jacobi equations

In this section we introduce the main ideas in our kernel-based method [22, 23], which is

derived from the Method of Lines Transpose (MOLT ) [19, 16, 18, 17, 21, 15]. 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.

3.1.1 1D Hamilton-Jacobi equations

We start with a 1D Hamilton-Jacobi equation

At + H(Ax) = 0,

A(x, 0) = A0(x),

x ∈ [a, b],

(3.1.1)

where H is the Hamiltonian flux. For the time discretization purpose to evolve the solution

from time tn to tn+1, we use the classical explicit SSP RK schemes [40]. In this work, we

propose to use the following SSP RK schemes such as the first order forward Euler scheme

An+1 = An − ∆t ˆH(An,−

x , An,+
x );

(3.1.2)

the second order SSP RK scheme

A(1) = An − ∆t ˆH(An,−

x , An,+
x ),

37

Å

A(1) − ∆t ˆH(A

(1),−
x

ã

)

;

(1),+
, A
x

(3.1.3)

An+1 =

An +

1
2

1
2

and the third order SSP RK scheme

Å
Å

A(1) = un − ∆t ˆH(An,−

A(2) =

An +

3
4

1
4

An+1 =

An +

1
3

2
3

x , An,+
x ),
A(1) − ∆t ˆH(A

A(2) − ∆t ˆH(A

(1),−
x

)

(1),+
, A
x
(2),−
x

(2),+
, A
x

ã
ã

,

)

.

(3.1.4)

Here, ∆t denotes the time step and ˆH is the numerical Hamiltonian, e.g., the local Lax-

Friedrichs Hamiltonian flux [46]:

ˆH(φ−, φ+) = H(

φ− + φ+

2

) − c(φ−, φ+)

(φ+ − φ−)

2

,

(3.1.5)

with c(φ−, φ+) = maxφ∈[min(φ−,φ+),max(φ−,φ+)] |H(cid:48)(φ)|. A−
tives with left-biased and right-biased methods, respectively, to approximate Ax. Finding

x are one-side deriva-

x and A+

these derivatives has the dominant role in our work and we will present the details of the

construction of them in the following subsection.

3.1.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 [22]. Let’s define operators LL

38

and LR on the closed interval [a, b], and their inverse operators

LL = I +
1
α
LR = I − 1
α

∂x ⇒ L−1
∂x ⇒ L−1

L [v, α](x) = α

R [v, α](x) = α

(cid:90) x
(cid:90) b

a

x

e−α(x−y)v(y)dy + ALe−α(x−a),
e−α(y−x)v(y)dy + BRe−α(b−x),

(3.1.6a)

(3.1.6b)

where I is the identity operator and α is a positive constant. We use IL and IR to denote

the convolution integral as

(cid:90) x

a

IL = α

e−α(x−y)v(y)dy,

IR = α

(cid:90) b

x

e−α(y−x)v(y)dy.

(3.1.7)

We can see that IL depends on the function values of v from left end point a to x, as well as

IR is for right end point b to x. Also note that AL and BR are determined by the boundary

conditions. For instance, if we assume periodic boundary conditions, i.e.,

L−1
L [v, α](a) = L−1

L [v, α](b),

then we obtain

and, L−1

R [v, α](a) = L−1

R [v, α](b).

AL =

IL[v, α](b)

1 − µ

,

and BR =

IR[v, α](a)

1 − µ

,

(3.1.8)

where µ = e−α(b−a). If we require

L−1
L [v, α](a) = Ca,

and L−1

R [v, α](b) = Cb

with given numbers Ca and Cb, then, the boundary coefficients are obtained as

AL = Ca,

and BR = Cb.

(3.1.9)

39

Now, let’s define new operators DL and DR,

DL = I − L−1
L ,

and DR = I − L−1
R .

Then 1

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

∂x = LL − I = LL(I − L−1

L ) = DL(I − DL)−1 =

∞(cid:88)

p=1

∂x = I − LR = LR(L−1

R − I) = −DR(I − DR)−1 = −

L

Dp
∞(cid:88)

p=1

1
α

1
α

(3.1.10a)

(3.1.10b)

Dp

R

where, Dp is successively defined as Dp = D[Dp−1]. In the kernel-based method, we will

truncate the series (3.1.10) and only compute the corresponding partial sum to approximate
the A±
x in (3.1.2)-(3.1.4). Moreover, an additional term D0 is introduced for reasons related

to stability of our high order method

D0[v, a] = v(a) − α
2

(cid:90) b

a

e−α|x−y|v(y)dy − A0e−α(x−a) − B0e−α(b−x),

(3.1.11)

where A0 and B0 are determined by boundary conditions. For example, if D0 is a periodic

function such that

D0[v, α](a) = D0[v, α](b),

and ∂xD0[v, α](a) = ∂xD0[v, α](b),

then we get

A0 =

I0[v, α](b)

1 − µ

, B0 =

I0[v, α](a)

1 − µ

,

(3.1.12)

40

with I0[v, α](x) = α
2

(cid:82) b
a e−α|x−y|v(y)dy. And if we require

D0[v, α](a) = Ca,

and D0[v, α](b) = Cb,

Ä
Ä

Ä
Ä

ä −Ä
ä −Ä

ää
ää

with given numbers Ca and Cb, then the boundary coefficients are obtained as

A0 =

B0 =

1

1 − µ2

1

1 − µ2

I0[v, α](b) − v(b) + Cb
I0[v, α](a) − v(a) + Ca

µ

µ

I0[v, α](a) − v(a) + Ca

I0[v, α](b) − v(b) + Cb

,

.

(3.1.13a)

(3.1.13b)

In fact, −α2D0[v, α] is an approximation to vxx, which we present the details in the section

(6.2).

In the following, we will study the properties of the partial sums with different boundary

conditions. Additionally, the parameter α would be chose carefully to make this scheme high

order of accuracy and unconditionally stable properties.

3.1.2.1 Periodic boundary conditions

If A is a periodic function, we define the approximations of the first derivatives A±

x using

partial sums as

x ≈ P−
A−

k [A, α](x) =



α(cid:80)k
p=1 Dp
α(cid:80)k
p=1 Dp

L[A, α](x),

k = 1, 2,

(3.1.14a)

L[A, α](x) − αD0 ∗ D2

L[A, α](x), k = 3,

41

and

x ≈ P+
A+

k [A, α](x) =



−α(cid:80)k
−α(cid:80)k

p=1 Dp

R[A, α](x),

k = 1, 2,

(3.1.14b)

p=1 Dp

R[A, α](x) + αD0 ∗ D2

R[A, α](x), k = 3.

Note that we have an extra term for the case of k = 3. This term is necessary to ensure

unconditional stability and show up in all MOLT approximations with k > 2. The scheme
for k > 2 is A(α) − stable, however with the modification shown above for k = 3 scheme
A − stable for properly chosen of the constant β [22, 23]. Then, we have the following

theorem, which gives an error estimate for the partial sums approximation and has been

proven in [22].

Theorem 3.1.1. Suppose v(x) ∈ Ck+1[a, b] is a periodic smooth function. Consider the
operators D∗ with the boundary treatment D∗(a) = D∗(b). Here, ∗ can be L, R, and 0.

Then,

k [A, α](cid:107)∞ = O(cid:16)

(cid:107)Ax − P−

1/αk(cid:17)

,

k [A, α](cid:107)∞ = O(cid:16)

(cid:107)Ax − P+

1/αk(cid:17)

.

(3.1.15)

Here, we take α = β/(c∆t), where c is the maximum wave speed and β is a constant

independent of the time step ∆t. Therefore, the accuracy for the approximation (3.1.14) to
the ∂xA is O(∆tk). Note that, to achieve kth order accuracy in time, k = 1, 2, 3, we should
employ the kth order SSP RK method (3.1.2)-(3.1.4) as well as P±
k in (3.1.14). Moreover, if
β is chosen appropriately, the semi-discrete scheme (without numerical approach to Dp) 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 [22].

42

Theorem 3.1.2. For the linear equation φt + cφx = 0, (i.e. the Hamiltonian is linear) with

periodic boundary conditions, we consider the kth order SSP RK method as well as the kth
partial sum in (3.1.14), with α = β/(|c|∆t). Then there exists a constant βk,max > 0 for
k = 1, 2, 3, such that the scheme is A-stable provided 0 < β ≤ βk,max. The constants βk,max

for k = 1, 2, 3 are summarized in Table 3.1.

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

k
βk,max

1
2

2
1

3
1.243

Proof. Let’s review the key points of proof for the case of k = 1. A von Neumann analysis

is applied to the advection equation to obtain the amplification factor λ for the given ansatz
φn = ˆφeikx. According to the von Neumann analysis if |λ| ≤ 1 for any k and time step ∆t,

then the scheme is unconditionally stable. Now we apply the Fourier transform in space to
the operators DL and LL

LL = I +

1
αL

∂x ⇒ ˆLL = 1 +

ik
αL

,

DL = I − L−1

L ⇒ ˆDL = 1 − 1
ˆLL

=

ik/αL

1 + ik/αL

.

Now, once we apply the Forward Euler scheme and integral approximation

where αL = β

(c∆t). Thus,

φn+1 = φn − ∆tαLDL[cφn, αL],

λ =

φn+1
φn = 1 − β

ikc∆t/β

1 + ikc∆t/β

.

43

(3.1.16)

(3.1.17)

(3.1.18)

(3.1.19)

For β ≤ 2, the amplification factor satisfies |λ| ≤ 1, which means the scheme is A - stable.

In other words, the scheme is unconditionally stable for appropriately chosen constant value

β from the table (3.1). The same argument holds for k > 1. This is discussed in [23, 24].

3.1.2.2 Non periodic boundary conditions

Next, we introduce a modification of the partial sums (3.1.14) to achieve a higher order

accuracy for the non-periodic case. Assume that we have some derivative values at boundary,
x A(b), m ≥ 1. Using integration by parts one can derive the following

i.e., ∂m
modified partial sums for k ≤ 3 to deal with the non-periodic boundary conditions

x A(a) and ∂m

k [A, α](x) =

k(cid:80)
k(cid:80)

p=1

p=1

α

α

DL[A1,p, α](x),
DL[A1,p, α](x) − αD0[A1,3, α](x), k = 3,

k = 1, 2,

(3.1.20a)

A+

k [A, α](x) =

k(cid:80)
k(cid:80)

p=1

p=1

DR[A2,p, α](x),
DR[A2,p, α](x) + αD0[A2,3, α](x), k = 3.

k = 1, 2,

(3.1.20b)

−α

−α




A−

x ≈‹P−
x ≈‹P+


And A1,p and A2,p are given as

A1,1 = A,

A1,2 = DL[A1,1, α] − k(cid:88)
k(cid:88)
A1,3 = DL[A1,2, α] +

m=2

åm

x A(a)e−α(x−a),
∂m
− 1
α

x A(a)e−α(x−a),
∂m

(3.1.21a)

Ç

− 1
α

åm
Ç

(m − 1)

m=2

44



Ç 1

åm

α

(m − 1)

Ç 1

åm

α

x A(b)e−α(b−x),
∂m

x A(b)e−α(b−x),
∂m

A2,1 = A,

A2,2 = DR[A2,1, α] − k(cid:88)
k(cid:88)
A2,3 = DR[A2,2, α] +

m=2

m=2

(3.1.21b)

where the boundary conditions for the operators are imposed as

αDL[A1,1, α](a) = Ax(a), αDR[A2,1, α](b) = −Ax(b),
αDL[A1,p, α](a) = αDR[A2,p, α](b) = 0,
αD0[A∗,3, α](a) = αD0[A∗,3, α](b) = 0,

for p ≥ 2,
∗ could be 1 or 2.

Then, the modified partial sum (3.1.20) agrees with the derivative values at the boundary

‹P−
k [A, α](a) = Ax(a), ‹P+

k [A, α](b) = Ax(b).

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

Theorem 3.1.3. Suppose A ∈ Ck+1[a, b]. Then the modified partial sums (3.1.20) satisfy

(cid:107)Ax −‹P−

(cid:107)Ax −‹P+

k [A, α](cid:107)∞ = O(1/αk),

k [A, α](cid:107)∞ = O(1/αk),

k = 1, 2, 3.

(3.1.22)

Again, we take α = β/(c∆t) in application, where c is the maximum wave speed and β

is taken as the same as periodic boundary condition, see Table 3.1.

As we have seen in (3.1.21), 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

45

are no physical boundary given). The details of the general non-periodic conditions can be

found in our previous work [23].

3.1.2.3 Space Discretization and WENO-based Quadrature

According to the kernel based method, the spatial derivatives can be approximated by the

partial sums based on an integral formulation. A fully discrete numerical solution is then
obtained by discretizing DL and DR operators in space.
In this subsection we will give
the details of the spatial discretization of the DL and DR operators and the WENO-based
quadrature formulation to approximate the convolution integrals appearing in the DL and
DR operators.

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

∆x = (b − a)/N,

xi = a + i∆x,

i = 0, 1, . . . , N.

The convolution integrals IL

i = IL[v, α](xi) and IR

i = IR[v, α](xi) satisfy a recursive relation

i = e−α∆xiIL
IL
i = e−α∆xi+1IR
IR

i−1 + J L
i ,

i+1 + J R
i ,

i = 1, . . . , N,

IL
0 = 0,
i = 0, . . . , N − 1,

IR
N = 0,

(3.1.23a)

(3.1.23b)

where,

(cid:90) xi

xi−1

J L
i = α

v(y)e−α(xi−y)dy,

J R
i = α

(cid:90) xi+1

xi

v(y)e−α(y−xi)dy.

(3.1.24)

Since there is a unique polynomial p(x) of degree at most k that interpolates v(x) at the

46

nodes in the interpolation stencil S(i) = {xi−r, . . . , xi−r+k}, which contains xi−1 and xi,
we can approximate J L

i with kth order accuracy by

(cid:90) xi

xi−1

i ≈ α
J L

p(y)e−α(xi−y)dy.

(3.1.25)

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

i can

be approximated by

(cid:90) xi+1

xi

i ≈ α
J R

p(y)e−α(y−xi)dy,

(3.1.26)

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

xi and xi+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 will use WENO-based quadra-

ture formula and the nonlinear filter to control oscillations and capture the correct solution,

which is proprosed in [23]. Here we will provide a brief description of methodology.

To summarize, if Ax is periodic function, then we will be using the following modified

sums framework instead of (3.1.14) for approximation of A±

x at xi:

A−
x,i = αDL[A, α](xi) + α

L[A, α](xi),

x,i = −αDR[A, α](xi) − α
A+

R[A, α](xi);

k(cid:88)
σp−2
i,L Dp
k(cid:88)

p=2

σp−2
i,R Dp

p=2

(3.1.27a)

(3.1.27b)

and if Ax is non periodic function, then we will use the following formulation instead of

47

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

(3.1.20)

A−
x,i = αDL[A1,1, α](xi) + α

k(cid:88)
σp−2
i,L DL[A1,p, α](xi),
k(cid:88)

σp−2
i,R DR[A2,p, α](xi).

(3.1.28a)

(3.1.28b)

x,i = −αDR[A2,1, α](xi) − α
A+

p=2

p=2

Moreover, we only use the WENO formulation when p = 1 while we use cheap high order

linear formulation (3.1.25) and (3.1.26) for the case p > 1. The filters σi,L and σi,R are

obtained based on the smoothness indicators from the WENO quadrature.

The fifth order 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) = {xi−3, . . . , xi+2} and the
three small stencils as Sr(i) = {xi−3+r, . . . , xi+r}, r = 0, 1, 2.

1. We approximate the integrals on each small stencils Sr(i) as follows

(cid:90) xi

xi−1

J L
i,r = α

e−α(xi−y)pr(y)dx, =

3(cid:88)

j=0

c

(r)
j vi−3+r+j,

(3.1.29)

where pr(x) is the polynomial interpolating v(x) on nodes Sr(i), and the coefficients

c

(r)
j depend on α and the cell size ∆x, but not on v.

48

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

(cid:90) xi

xi−1

J L
i = α

e−α(xi−y)p(y)dx =

with the linear weights dr satisfying (cid:80)2

r=0 dr = 1.

2(cid:88)

r=0

drJ L
i,r,

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

ωr = ˜ωr/

with

2(cid:88)

s=0

Ç

˜ωs, r = 0, 1, 2,

å

˜ωr = dr

1 +

τ5


 + βr

.

(3.1.30)

(3.1.31)

We take 
 = 10−6 as a small positive number, 
 > 0, in our numerical test problems

to avoid zero at the denominator. The smoothness indicator βr is determined as

3(cid:88)

(cid:90) xi

xi−1

l=2

βr =

∆x2l−3

i

(cid:32) ∂lpr(x)

(cid:33)2

∂xl

dx,

(3.1.32)

which is used to measure the relative smoothness of the function v(x) in the stencil
Sr(i). Here, τ5 = |β0 − β2|. Furthermore, we introduce a parameter ξi as

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

(3.1.33)

å2

.

Ç

τ5


 + max(β0, β2)

ξi =

βmin
βmax

,

å2

,

βmin = 1 +

49

Ç

βmax = 1 +

τ5


 + min(β0, β2)

Note that we have developed the nonlinear weights using the idea of the WENO-Z

method proposed in [13], which has less dissipation and higher resolution comparing

to the original WENO method.

4. Lastly, we obtain the approximation

2(cid:88)

r=0

J L
i =

ωrJ L
i,r.

(3.1.35)

The filter σi,L is determined as

σi,L = min(ξi−1, ξi).

(3.1.36)

All coefficients are given in Appendix. 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 xi.

To close this section, we summarize the proposed method for approximating one-dimensional

problem (3.1.1) in the following algorithm flowchart.

Given the function un, the approximation order k ≤ 3, the mesh size ∆x and the time step ∆t.

1. Choose β from the Table 3.1. Compute c = max|H(cid:48)(φ)| at time tn.

2. On each inner stage of the kth order SSP RK scheme:

(a) Compute c, and then obtain parameters α = β/(c∆t). For non periodic boundary

conditions, compute the derivatives ∂m

x A(a) and ∂m

x A(b), m = 1, . . . , k.

50

(b) Apply DL and DR on A, respectively. Use the WENO quadrature to calculate J L,R
and at the same time obtain the nonlinear filter σL,R. Compute IL,R via (6.2.6) and

then calculate the parameter AL and BR based on the boundary condition. Combine
IL,R, AL and BR to construct DL[A] and DR[A].

(c) For k > 1, further construct Dp

R [DR[A]], or
the modified functions DL,R[A1,2,k] in (3.1.21), by a similar procedure for 1 < p ≤ k.

R[A] = Dp−1

L[A] = Dp−1

L [DL[A]] and Dp

WENO quadrature is not needed for construction of these high order terms.

(d) Construct A±

x based on the the partial sum approximations (3.1.14) or (3.1.20).

(e) Substitute into the Lax-Friedrichs Hamiltonian flux (3.1.5) and update the solution

accordingly.

3.1.3 2D magnetic potential equation

In this section we develop our proposed MOLT-HJ scheme in 2D, which still use the WENO-

quadrature formulation presented in the previous section to approximate the convolution

integrals. According to the Constrained Transport formulation described in Section 2.5, we

must update the solution of the magnetic potential equation by solving a discrete version of

the following equation

A3
t + u1(x, y)A3

x + u2(x, y)A3

y = 0

(3.1.37)

where the velocity components u1 and u2 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

51

(3.1.37) as a Hamilton-Jacobi equation

A3
t + H(A3

x, A3

y) = 0

with Hamiltonian flux

H(A3

x, A3

y) = u1(x, y)A3

x + u2(x, y)A3
y.

(3.1.38)

(3.1.39)

We can directly apply a two dimensional version of the kernel-based framework presented in

Section 3.1.1. The 2D semi-discrete scheme can be written as

dA3
i,j(t)
dt

= − ˆH(A3−

x |i,j, A3+

x |i,j, A3−

y |i,j, A3+

y |i,j)

(3.1.40)

on each point (xi, yj), where ˆH is a Lipschitz continuous Hamiltonian flux. We use the

Lax-Friedrichs flux and the equation (3.1.40) becomes

Å A3−
ÅA3+

x + A3+
x
x − A3−

2

x

2

ã|i,j − u2
ÅA3−
ã|i,j + c2
ÅA3+

i,j

2

y + A3+
y

2

y − A3−

y

ã|i,j
ã|i,j

dA3
i,j(t)
dt

= − u1

i,j

+ c1

with

(3.1.41)

c1 = max
i,j

i,j| and c2 = max
|u1

i,j

i,j|.
|u2

We remark that the scheme (3.1.41) with this global cm can be very dissipative for some

Hamilton-Jacobi equations. There is another way of choosing cm which is by taking the max

52

value on the local stencil.

The approximations A3±

x |i,j and A3±

y |i,j to the derivatives of functions Ax(x, y) and

Ay(x, y) at (xi, yj), respectively calculated directly using one dimensional formulation of
the scheme, e.g., when computing A3±

x , we fix y and apply 1D scheme in x-direction. For

example, when Ax is periodic in x-direction, we get

A−
x ≈ α1Dp

L[A(·, y), α1](x) + α1

L[A(·, y), α1](x),

x ≈ −α1DR[A(·, y), α1](x) − α1
A+

σp−2
i,R Dp

R[A(·, y), α1](x).

Here, we choose α1 = β/(c1∆t). Similarly, to approximate A±

y , we fix x and obtain

k(cid:88)
σp−2
i,L Dp
k(cid:88)

p=2

p=2

k(cid:88)
σp−2
i,L Dp
k(cid:88)

p=2

p=2

A−
y ≈ α2DL[A(x,·), α2](y) + α2

L[A(x,·), α2](y),

y ≈ −α2DR[A(x,·), α2](y) − α2
A+

σp−2
i,R Dp

R[A(x,·), α2](y),

with α2 = β/(c2∆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.

At an outflow boundary condition, we still use extrapolation with suitable order of ac-

curacy for the derivative values, as in the 1D formulation. In the case of outflow boundary

conditions, in the sense of characteristic propagation, it is natural to apply high order ex-

trapolation to get the boundary values.

53

3.1.4 3D magnetic potential equation

In this section we present MOLT-HJ scheme in 3D, which is a key scheme in our proposed

3D MHD solver. The main accomplishment in 3D MOLT-HJ scheme is that we avoid the

artificial diffusion limiters, which were needed in the previous work [29]. Note that we can

include those artificial limiters to our 3D MOLT-HJ solver to compare the new method with

the old method. Although the evolution equation for the 3D magnetic potential (2.3.4)

is significantly different from the evolution equation for the 2D scalar magnetic potential

(2.3.9), we can still directly apply the scheme presented in section (3.1.1) to the 3D case.

In our previous paper [29], we needed an artificial resistivity term in the 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.3) in component

form we have:

∂tA1 = u2(∂xA2) + u3(∂xA3) − u2(∂yA1) − u3(∂zA1)

∂tA2 = −u1(∂xA2) + u1(∂yA1) + u3(∂yA3) − u3(∂zA2)

∂tA3 = −u1(∂xA3) − u2(∂yA3) + u1(∂zA1) + u2(∂zA2)

(3.1.42)

(3.1.43)

(3.1.44)

with A = (A1, A2, A3) and u = (u1, u2, u3). 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:

dA1(t)

dt

= u2 (A2+

x + A2−
x )

2

x + A3−
x )

+ u3 (A3+

2

− c3

(A3−

x − A3+
x )

2

− c2

(A2−

x − A2+
x )

2

54

y + A1−
y )

− u2 (A1+

2

dA2(t)

dt

dA3(t)

dt

= − u1 (A2+
+ u3 (A3+

x + A2−
x )
y + A3−
y )

2

2

= − u1 (A3+
+ u1 (A1+

x + A3−
x )
z + A1−
z )

2

2

− c2

− c1
− c3

− c1
− c1

(A1−

y − A1+
y )

2

(A2−
(A3−

x − A2+
x )
y − A3+
y )

2

2

(A3−
(A1−

x − A3+
x )
z − A1+
z )

2

2

− u3 (A1+

z + A1−
z )

2

+ u1 (A1+
− u3 (A2+

y + A1−
y )
z + A2−
z )

2

2

y + A3−
y )
z + A2−
z )

− u2 (A3+
+ u2 (A2+

2

2

− c3

− c1
− c3

− c2
− c2

(A1−

z − A1+
z )

2

,

(3.1.45a)

(A1−

y − A1+
y )
z − A2+
z )

2

(A2−

2

,

(3.1.45b)

(A3−
(A2−

y − A3+
y )
z − A2+
z )

2

2

,

(3.1.45c)

at (xi, yj, zk), where

c1 = max
i,j,k

i,j,k|,
|u1

c2 = max
i,j,k

i,j,k|,
|u2

and c3 = max
i,j,k

i,j,k|.
|u3

Similarly, we use 1D kernel-based formulation to approximate the derivatives of the mag-

netic 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 [29]. For more details on positivity limiter, see [26, 28, 66]

3.2 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

55

curl operator.

3.2.1 Curl in 2D

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

i,j := Dy|i,jA3
B1

and

i,j := −Dx|i,jA3,
B2

(3.2.1)

where Dx and Dy 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

∇ · Bi,j := Dx|i,jB1 + Dy|i,jB2 = DxDy|i,jA3 − DyDx|i,jA3 = 0.

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

Dx|i,jA3 :=
Dy|i,jA3 :=

1

12∆x

1

12∆y

(A3

(A3

i−2,j − 8A3
i,j−2 − 8A3

i−1,j + 8A3

i,j−1 + 8A3

i+1,j − A3
i,j+1 − A3

i+2,j),

i,j+2),

(3.2.2a)

(3.2.2b)

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

3.2.2 Curl in 3D

We use the following discrete version of the 3D curl

i,j,k := Dy|i,j,kA3 − Dz|i,j,kA2,
B1
i,j,k := Dz|i,j,kA1 − Dx|i,j,kA3,
B2

(3.2.3a)

(3.2.3b)

56

i,j,k := Dx|i,j,kA2 − Dy|i,j,kA1,
B3

(3.2.3c)

where Dx, Dy, and Dz 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

∇ · Bi,j,k : = Dx|i,j,kB1 + Dy|i,j,kB2 + Dz|i,j,kB3

= DxDy|i,j,kA3 − DxDz|i,j,kA2 + DyDz|i,j,kA1 − DyDx|i,j,kA3
+ DzDx|i,j,kA2 − DzDy|i,j,kA1

(3.2.4)

= 0

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

Dx|i,j,kA∗ :=
Dy|i,j,kA∗ :=
Dz|i,j,kA∗ :=

1

12∆x

1

12∆y

1

12∆z

i−2,j,k − 8A∗
(A∗
i,j−2,k − 8A∗
(A∗
i,j,k−2 − 8A∗
(A∗

i−1,j,k + 8A∗
i,j−1,k + 8A∗
i,j,k−1 + 8A∗

i+1,j,k − A∗
i,j+1,k − A∗
i,j,k+1 − A∗

i+2,j,k),

i,j+2,k),

i,j,k+2).

(3.2.5a)

(3.2.5b)

(3.2.5c)

3.3 Positivity Limiter

In this section we will give an explanation of the positivity limiter strategy that is used in

several test problems in this work. Some test problems that has strong shock sometimes

cause negative density or pressure in simulations. These negative quantities can create

complex wave speeds which violates the hyperbolicity of the MHD system and ruin the

numerical simulations. In the literature there have been many works addressing this problem

57

[7, 80, 77, 76, 27, 66]. A positivity limiting strategy was developed in the previous work [29]

is also being utilized in this work to avoid negative density and pressure.

We will present a brief overview of the positivity limiter on one-dimensional MHD system,

∂q
∂t

+

∂f (q)

∂x

= 0

(3.3.1)

The base scheme WENO-HCL is applied to the MHD equations to update the conserved

quantities with a conservative form as

dqi(t)

dt

+

1
∆x

Ç

ˆF

i+ 1
2

− ˆF

i− 1
2

= 0

(3.3.2)

å

å

where ˆF

i+ 1
2

is a numerical flux obtained by WENO-HCL scheme. We can apply time

discretization to the semi-discrete equation (3.3.2) with third-order SSP-RK method

L(qn

i ) = − 1
∆x

( ˆFn

i+ 1
2

− ˆFn

i− 1
2

.

(3.3.3)

The numerical fluxes reconstructions based on qn, q(1), and q(2) are ˆFn

, ˆF

(2)
i+ 1
2
a linear combination of numerical fluxes from different

, and ˆF

(1)
i+ 1
2

i+ 1
2

,

respectively. Thus, considering ˆFrk
i+ 1
2

stages, the final stage of RK scheme can be written as

qn+1
i

= qn

i − λ( ˆFrk
i+ 1
2

− ˆFrk
i− 1
2

)

+ 4 ˆF

(2)
i+ 1
2

+ ˆF

(1)
i+ 1
2

å

.

(3.3.4)

(3.3.5)

where

λ =

∆t
∆x

,

ˆFrk
i+ 1
2

=

1
6

ˆFn

i+ 1
2

Ç

58

The idea of the positivity limiter is to modifying the numerical flux ˆFrk
i+ 1
2

by a positivity

preserving flux to guarantee the positivity of the high order solutions by the WENO method

as follows

˜F

i+ 1
2

= θ

( ˆFrk
i+ 1
2

i+ 1
2

− ˆf

i+ 1
2

) + ˆf

i+ 1
2

(3.3.6)

with the limiting parameter θ

is the Lax-Friedrichs flux and

∈ [0, 1]. Here, the ˆf

i+ 1
2

i+ 1
2

Ä

ä

formulated as

ˆf
i+ 1
2

=

1
2

f (qn

j+1) + f (qn

i ) − α(qn

j+1 − qn
i )

where the α is the maximal wave speed

α = max

i

Ä|ux| + cf

ä

(3.3.7)

(3.3.8)

with cf as the fast speed of the MHD system. There are two steps for the positivity preserving

limiting procedure. The first step is to guarantee the positivity of the computed density.

We will give the following notations: ˆf ρ is the first order flux of the density in ˆf ; f ρ and ˜f ρ

are the flux components in Frk and ˜F, respectively. We seek upper bounds Λρ

± 1

2 ,Ii

limiting parameters θ

i± 1
2

at each cell Ii so that

(θ

i− 1
2

, θ

i+ 1
2

) ∈ Sρ,Ii

= [0, Λρ

− 1

2 ,Ii

] × [0, Λρ
+ 1

2 ,Ii

]

of the

(3.3.9)

and the computed density is being positive.

ρn+1
i

(θ

i− 1
2

, θ

i+ 1
2

) ≥ 
n+1

p

.

(3.3.10)

59

The second step is to preserve positive pressure. We describe a subset of the Sρ,Ii

as

Sp,Ii

= {(θ

i− 1
2

, θ

i+ 1
2

) ∈ [0, Λρ
− 1

2 ,Ii

] × [0, Λρ
+ 1

2 ,Ii

]}

(3.3.11)

and then we describe a rectangle inside of Sp,Ii

as

Rρ,p,Ii

= [0, Λ− 1

2 ,Ii

] × [0, Λ

].

+ 1

2 ,Ii

such that the pressure function pn+1

i

(θ

i− 1
2

, θ

i+ 1
2

) is positive,

pn+1
i

(θ

i− 1
2

, θ

i+ 1
2

) ≥ 
n+1

p

.

At the end of the this procedure for all j, we get

(3.3.12)

(3.3.13)

θ

i+ 1
2

= min(Λ− 1

2 ,Ii

, Λ

)

+ 1

2 ,Ii

(3.3.14)

The positivity preserving limiter strategy is used to guarantee the positivity of the solution

at the final RK stage. So, we preserve numerical solutions with positive density and pressure

at the end of the step 1 in the Constrained Transport procedure. Also, we use option 2 in

the CT framework to preserve positive pressure after magnetic field is updated.

The details are omitted here. For the details, see [29].

60

3.4 Numerical Results

In this section, we present the numerical results to demonstrate the accuracy and efficiency

of the new method. We demonstrate the robustness of the proposed method on the field

loop test problem in 2D and 3D. We validate that the new method works well for the blast

wave and strong shocks by comparing with the old method [29] and then adding diffusion

limiter into the new method. The main accomplishment is that the new method achieves a

peak of 320, which is 23% higher than the old method, and thus is far less isotropic around

the peak than the old method. We use third order SSP RK method for time discretization.

Time step is chosen as

∆t =

CFL
em ∗ cd

(3.4.1)

where the β = 0.8, 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. In all examples, unless stated, the results being presented are from

kernel based-CT. We denote the new method as kernel based-CT and our WENO-HJ-CT

method as our previous method.

3.4.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

(ρ, u1, u3, u3, p, B1, B2, B3) = (1, 1, 1, 0, 1, 0, 0, 0)

61

with perturbations on u1, u2, B1, B2 and p as :

(δu1, δu2) =

(δB1, δB2) =

e0.5(1−r2)(−y, x).
µ
2π
e0.5(1−r2)(−y, x).
κ
2π

µy(1 − r2) − κ2

e(1−r2).

δp =

8π2

And the initial condition for magnetic potential is

A3(0, x, y) =

e0.5(1−r2)

µ
2π

where r2 = x2 + y2. 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 3.2, we present the errors of B at t = 0.05 with the mesh size
160 × 160. We observe fourth order accuracy in space for B. Note that we use fifth order

WENO scheme for the solution of A, but since we take derivative of the magnetic potential

equation to get magnetic field equation, we loose an order of accuracy for B, thus we get

fourth order in space.

Additionally, in Figure 3.2 we show the results at the final time t = 20 with 160 × 160

mesh size. Those plots show the solution after one period at t = 20.

3.4.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.

62

Table 3.2: Smooth vortex problem. Errors of B and orders of accuracy.

Nx × Ny
20 × 20
40 × 40
80 × 80
160 × 160
320 × 320

L1 error
2.827E-03
2.982E-04
1.861E-05
1.102E-06
7.260E-08

order

–

3.245
4.003
4.078
3.924

L∞ error
1.479E-01
1.839E-02
1.327E-03
1.126E-04
1.170E-05

order

–

3.007
3.793
3.559
3.267

The initial conditions are given as

(ρ, u1, u2, u3, p, B1, B2, B3) = (γ2,− sin(y), sin(x), 0, γ,− sin(y), sin(2x), 0),

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

A3 = 0.5 cos(2x) + 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 3.3, we show density ρ at time t = 3, and

compare the results with our previous method [29]. We can see that they are in good

agreement.

63

(a) (cid:107)u(cid:107)2

(b) (cid:107)B(cid:107)2

Figure 3.2: Smooth vortex problem. L2 norms of u and B at t = 20 with 160 × 160 grid
points. Domain for the plots is [−5, 5] × [−5, 5].

3.4.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

(ρ, u1, u2, u3, p, B1, B2, B3) =

 (3.86859, 11.2536, 0, 0, 167.345, 0, 2.1826182,−2.1826182) x < 0.05,

(1, 0, 0, 0, 1, 0, 0.56418958, 0.56418958)

x > 0.05,

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

64

(a) Kernal-based method.

(b) Previous method [29].

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

for magnetic potential:



A3 =

x ≤ 0.05,
−2.1826182(x − 0.05)
−0.56418958(x − 0.005) x ≥ 0.05,

Figure 3.4 presents Schlieren plots of ||B|| at t = 0.06. The new method matches well with

our numerical results of our previous method [29].

3.4.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

65

(a) Kernal-based method.

(b) Previous method [29].

Figure 3.4: Cloud shock problem. Contour plots of (cid:107)B(cid:107) at t = 0.06 with 512 × 512 grid
points.

the positivity preserving limiter we developed in reference [29]. The initial conditions are

√
(ρ, u1, u2, u3, B1, B2, B3) = (1, 0, 0, 0, 50/

√
2π, 50/

2π, 0)

with a spherical pressure pulse

p =

 1000 r ≤ 0.1,

0.1

otherwise.

The initial condition for magnetic potential can be obtained as

√
A3 = (0, 0, 50/

2π(y − x)).

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 3.5 and they have good agreement

with the results that we got in our previous paper.

66

(a) density ρ

(b) pressure p

(c) (cid:107)u(cid:107)

(d) (cid:107)B(cid:107)

Figure 3.5: 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 (cid:107)B(cid:107).

67

 0.001(R − r), r ≤ R,

0,

A3 =

otherwise,

»

where r =

3.4.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

(ρ, u1, u2, p) = (1,

√

5 cos(θ),

√

5 sin(θ), 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

x2 + y2 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 3.6. We observe that the solution maintains the circular symmetry

of the initial condition.

3.4.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

√
(ρ, u1, u2, u3, p) = (1, 2/

√
6, 1/

√
6, 1/

6, 1)

68

(a) (cid:107)B(cid:107)

(b) A

Figure 3.6: 2D Field loop. Contour plots at time t = 1 with 128 × 128 grid points. a)
Magnetic Pressure (cid:107)B(cid:107); b) Magnetic Potential A3.

The initial conditions for magnetic field are determined by taking the curl of the magnetic

potential, which is given as magnetic potential :

 0.001(R − r) r ≤ R
»

0

A3 =

otherwise,

x2 + y2 and R = 0.3. We use a domain size [−0.5, 0.5] ×
where A1 = 0, A2 = 0, and r =
[−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 3.7 are a 2D

slice of the 3D solution taken at z = 0.

69

(a) (cid:107)B(cid:107)

(b) A

Figure 3.7: 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.

3.4.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

√
(ρ, u1, u2, u3, B1, B2, B3) = (1, 0, 0, 0, 50/

√
2π, 50/

2π, 0)

with a spherical pressure pulse

p =

 1000 r ≤ 0.1

0.1

otherwise.

70

»

where r =

x2 + y2 + z2. The initial condition for magnetic potential :

√
A(0, x, y, z) = (0, 0, 50/

√
2πy − 50/
2π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 3.8 and Figure

3.9 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

peak than the old method. As demonstrated in a latter test, the isotropic behavior and

lower the peak in the strong bast wave is due to the diffusion limiter that was needed with

the old explicit constrained transport method.

3.4.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 [29], 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 (cid:107)u(cid:107)

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

3.10, 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 (cid:107)u(cid:107) is 265. These results confirm that the

71

(a) density with kernel-based method

(b) density with previous method [29].

(c) pressure with kernel-based method

(d) pressure with previous method [29]

Figure 3.8: 3D Blast wave problem. Contour plots at time t = 0.01 with 150 × 150 grid
points

72

(a) (cid:107)u(cid:107) with kernel-based method

(b) (cid:107)u(cid:107) with previous method [29]

(c) (cid:107)B(cid:107) with kernel-based method

(d) (cid:107)B(cid:107) with previous method [29]

Figure 3.9: 3D Blast wave problem. Contour plots at time t = 0.01 with 150 × 150 grid
points. While the max (cid:107)u(cid:107) value for previous method is 261, it is 320 for the kernel-based
method. Note that this is 23% higher than the old method, and is far less isotropic around
the peak than the old method

.

73

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.

74

(a) density

(b) pressure

(c) (cid:107)u(cid:107)

(d) (cid:107)B(cid:107)

Figure 3.10: 3D Blast wave from the kernel based method with having added the diffusion
limiter. The max (cid:107)u(cid:107) value is 265. We note that these results are nearly identical to the old
WENO-HJ method, which required the diffusion limiter.

75

Chapter 4

MOLT for PEC for ideal MHD

In this chapter we expand our proposed scheme to work for perfectly electrically conducting

(PEC) boundary conditions. We develop PEC boundary conditions for magnetic potential

in MOLT-HJ scheme. Then, we go over the derivation of a key test problem called the

Hartmann Flow problem and demonstrate the new method on this test problem.

4.1 MOLT for PEC for ideal MHD

We have seen in the previous chapters that our kernel based approximation is working

well for the periodic and outflow boundary conditions.

In this section, we will explain

the numerical boundary condition for bounded domains in ideal MHD. We will define the

perfectly electrically conducting (PEC) boundary conditions for the conserved variables and

the magnetic potential. In particular, we will derive these boundary conditions in the context

of MOLT and show that it is equivalent to enforcing Dirichlet conditions.

Following the ideas in [25], we can define the PEC boundary condition. Let’s assume

that the PEC boundary exists at ξ = ξ0. Then on a PEC boundary, we have the following

conditions,

u1(ξ0, η, t) = 0,

B1(ξ0, η, t) = 0.

76

(4.1.1)

(4.1.2)

According to the local characteristics analysis, forward fast magnetosonic waves is moving

into the domain, while backward fast magnetosonic waves is exiting the domain. Remaining

characteristics are shown to propagate along the boundary. Note that the PEC boundary

condition on the magnetic field (4.1.2) guarantees that the remaining four characteristics,

Alfv´en and slow magnetosonic waves, are moving along the boundary. Thus, these PEC

boundary conditions are consistent with the analysis of local characteristic. A numerical

implementation of these boundary conditions requires a careful treatment, as many charac-

teristics are moving parallel to the boundary.

We will describe the analytic boundary condition for the magnetic potential A in the

constrained transport method before discussing the numerical boundary conditions. The

PEC boundary condition for magnetic field (4.1.2) at the boundary ξ = ξ0 satisfies the

following relation

B1 = ∂xξ∂yA − ∂yξ∂xA,

(4.1.3)

= (∂xξ∂yη − ∂xη∂yξ)∂ηA,

= 0,

which gives that

∂ηA = 0.

(4.1.4)

Hence, the magnetic potential equation at the boundary satisfies

∂tA + (u1∂xξ + u2∂yξ)∂ξA = 0.

(4.1.5)

Also, due to the velocity boundary condition (4.1.1), along with (u1∂xξ +u2∂yξ)∂ξ = ∇xξ·u,

77

we obtain the analytical boundary condition for the magnetic potential as

A(ξ0, η, t) = A0,

(4.1.6)

where A0 is some constant number that is given in the initial condition.

A Dirichlet boundary condition is applied to the magnetic potential

Ai = A0,

on i = 0.

(4.1.7)

Note that values at the ghost points are approximated, here, by using extrapolation of

sufficiently high order. More details can be found in [25].

4.1.1 Derivation of the Boundary Conditions for MOLT for Dirich-

let BC

In the previous section, we showed that the PEC boundary condition is equivalent to en-

forcing a Dirichlet boundary condition on the magnetic vector potential. Our next task is

to connect this to the MOLT scheme using a simple example.

we define the operator LR as

∂tv − c∂xv = 0,

LR = I − 1
α

∂x,

78

(4.1.8)

(4.1.9)

where α = β

c∆t. For the Backward Euler time discretization, we obtain

vn+1 − vn

∆t

− c∂xvn+1 = 0,

(4.1.10)

where tn is the current time level and tn+1 is the new time step. Then,

(vn+1 − vn) − c∆t∂xvn+1 = 0,

(4.1.11)

then if we rearrange terms, we obtain the semi-discrete equation

(1 − c∆t∂x)vn+1 = vn,

which can be compactly stated as

LRvn+1 = vn.

(4.1.12)

(4.1.13)

We can then solve this equation by applying the integrating factor method. With the aid of

equation (4.1.13), we find that

∂x(e−αxvn+1(x)) = −αe−αxvn+1(x) + e−αx∂xvn+1(x),

(4.1.14)

= −αe−αx(vn+1(x) − 1
α

∂xvn+1(x)),

= −αe−αxvn(x).

We integrate the relation (4.1.14) from x to b, which result in

(cid:90) b

x

∂x(e−αyvn+1(y))dy = −α

(cid:90) b

x

79

e−αyvn(y)dy,

(4.1.15)

after which, we see that

e−αbvn+1(b) − e−αxvn+1(x) = −α

(cid:90) b

x

e−αyvn(y)dy.

(4.1.16)

If we rewrite this equation to solve for vn+1, we obtain the formula

(cid:90) b

x

vn+1(x) = α

e−α(y−x)vn(y)dy + e−α(b−x)vn+1(b).

(4.1.17)

In the above, we use the definitions

IR[vn](x) = α

(cid:90) b

x

e−α(y−x)vn(y)dy

and

BR = vn+1(b).

(4.1.18)

The equation (4.1.17) is actually vn+1(x) = L−1
series of DR, which is defined as DR = I − L−1
of DR to approximate the first order derivative term as follow in short

R [vn](x) and we need L−1
R . Remember that we use the infinite series

R in our infinite

(cid:88)Dp

R.

∂x = −α

(4.1.19)

Now we show how to approximate the boundary term BR using a Taylor expansion. For

non-linear problems, we have observed that using linear operators results in inconsistencies

at the boundaries. Using a Taylor expansion around the previous time level’s data allows us

to correct this issue. Applying Taylor expansion on vn+1 gives

BR := vn+1(b) = vn(b) + ∆tvn

t (b) + ...

.

(4.1.20)

80

Employing Cauchy-Kowalevski theorem and equation (4.1.8), we have

vn
t = cvn

x + O(∆x2),

(4.1.21)

hence equation (4.1.20) becomes

BR = vn(b) + c∆tvn

x (b) + O(∆x2∆t),

(4.1.22)

= vn(b) +

1
α

vn
x (b).

In the next section, we describe a well-known benchmark problem, called Hartmann

Flow, to test our formulation of PEC boundary conditions in the MOLT using the above

correction.

4.2 Hartmann Flow

In general, the non-linearity of the MHD equations makes it difficult to write analytical

solutions. However, there is a restricted version of MHD model where non-linear terms

of MHD equations are excluded. For this special case of the MHD model, with the aid

of stationary conditions, it is possible to obtain a steady-state analytical solution. One

such example of this is the so-called Hartmann Flow problem, in which the flow is taken

to be incompressible. The importance of this test problem is that it allows us to assess

the capabilities of our code for both restrictive and ideal MHD models with PEC boundary

conditions. Even though the model for Hartmann Flow includes diffusion terms, we can

ignore them to test the PEC boundaries in ideal MHD. Here, we provide a brief derivation

81

of the Hartmann Flow problem, which includes the diffusion terms.

We will start with the MHD equations and end up with the Hartmann Flow. Since the

model for MHD is taken to be incompressible,

∇ · v = 0.

(4.2.1)

Further, in Hartmann flow, ρ = 1 in the domain as well as p = 1. Additionally, the

momentum and magnetic field equations are given as

ñ ∂v

∂t

ô

ρ

+ (v · ∇)v

= −∇p + f + N (j × B) +

∇2v

1
R

∂B
∂t

+ (v · ∇)B = (B · ∇)v +

∇2B

1

Rm

(4.2.2)

(4.2.3)

where the last term at the right hand side of the equation (4.2.3) represents magnetic dif-

fusion. With the assumptions of linearity in equations (4.2.2) and (4.2.3), we have that

v = vx(y)ˆx and B = ˆy + Rm

Ha Bx(y)ˆx, where the subscript x denotes the x-component of

the velocity and magnetic field, respectively. Note that, these components are functions of

a single variable y only.

The geometry of the problem consists of two infinite parallel plates and there is flow

in the gap between the plates. The externally applied magnetic field is perpendicular to

the parallel plates, and there is no boundary at the x direction. Because of the stationary

condition and the assumptions above, left hand side of equation (4.2.3) is identically zero.

Thus, the induction equation simplifies to

∂2Bx
∂y2 + Ha

∂vx
∂y

= 0.

82

(4.2.4)

Figure 4.1: Sketch of a Hartmann type of flow [58].

One can also verify that, under these assumptions, the left hand side of the momentum

equation (4.2.2) vanishes. Terms on the right hand side can be multiplied by R so that the

magnitude of the gradient of the uniform pressure becomes unity. Also, the Ampere’s law

gives a relation between the induced magnetic field and the current density

j =

1
Ha

∇ × Bxˆx =

−∂Bx
∂y

1
Ha

ˆz,

(4.2.5)

so that the Lorenz force j × B can be given with respect to Bx. Therefore, the momentum

equation simplifies to

∂2vx
∂y2 + Ha

∂Bx
∂y

= −1.

Thus, we obtain the Hartman flow problem,

∂2vx
∂y2 + Ha

∂Bx
∂y

= −1,

∂2Bx
∂y2 + Ha

∂vx
∂y

= 0, −1 < y < 1,

and

vx = 0, ±∂Bx
∂y

= 0,

at

y = ±1,

83

(4.2.6)

(4.2.7)

(4.2.8)

(4.2.9)

Ç

å

where the walls are perfectly conducting. The analytic solutions of the equations above can

be found as

vx(y) =

−1

Ha2cosh(Ha)

cosh(Ha y) − cosh(Ha)

Bx(y) =

−y
Ha

+

1

Ha2

sinh(Ha y)
cosh(Ha)

(4.2.10)

(4.2.11)

More details can be found in [58]. Moreover, we can derive the magnetic potential equation

in 2D as

−y2
2Ha

A =

+

1

Ha3

cosh(Ha y)
cosh(Ha)

.

(4.2.12)

4.2.1 Numerical Results for ideal MHD with PEC Boundary Con-

ditions

A refinement study is done for the Hartmann Flow problem for the ideal MHD, where the

diffusion effect is zero, η = 0 and the Hartmann number is Ha = 10. We use a CFL number

of 0.5 and β = 0.8. In Table (4.1), we obtained high order accuracy for magnetic field, in

the case of ideal MHD with PEC boundary conditions.

Table 4.1: Hartmann Flow for ideal MHD. Errors of B and orders of accuracy.
order

L∞ error

L1 error

order

9.6462E-007
7.2391E-008
4.9457E-009
3.2161E-010

–

3.74
3.87
3.94

5.3007E-006
7.9204E-007
1.0718E-007
1.1407E-008

–

2.74
2.89
3.23

Nx × Ny
64 × 64
128 × 128
256 × 256
512 × 512

In Figure 4.2, we have 1D cross-sections of the magnetic potential A and the induced
magnetic field Bx for the resolution 128 × 128 and the final time t = 1. We see that exact

solutions overlap the numerical results.

84

(a) A

(b) Bx

Figure 4.2: 1D cross-sections of Hartmann problem for Ha = 10 at t = 1 with 128 × 128.

In Figures 4.3, 4.4, 4.5, we have 2D contour plots of exact and numerical results of ux, Bx,

and A, respectively. The numerical results have good agreement with the exact solutions.

(a) ux exact.

(b) ux obtained with the MOLT scheme

Figure 4.3: Velocity profile for the Hartmann problem for Ha = 10 in 2D at t = 1 with
128 × 128.

85

(a) Bx exact solution.

(b) Bx obtained with the MOLT scheme

Figure 4.4: Induced magnetic field for the Hartmann problem for Ha = 10 at t = 1 with
128 × 128.

(a) A exact solution.

(b) A obtained with the MOLT scheme

Figure 4.5: Magnetic potential A for the Hartmann problem at Ha = 10.

86

Chapter 5

Non-Ideal MHD Equations

5.1 NON- IDEAL MHD

In this chapter we will describe the non-ideal MHD equations. The structure of this system

is similar to that of the ideal case, with the difference being the introduction of certain terms

in the equation for the magnetic field. These terms, which come from a generalized form of

Ohm’s law, result in the collection of equations given by





ρ

ρu
E

B

∂t

+ ∇ ·







ρu

ρu ⊗ u + ptotI − B ⊗ B
u(E + ptot) − B(u · B)

=

u ⊗ B − B ⊗ u



0

0

0
∇ × (ηJ + 1

ne J × B)

,

(5.1.1)

∇ · B = 0,

(5.1.2)

with the equation of state as

E =

p

γ − 1

+

ρ(cid:107)u(cid:107)2

2

(cid:107)B(cid:107)2
2

.

+

87

Here, ηJ is commonly referred to as the resistive term, while 1

neJ × B is known as the Hall

term. In this work, we focus on the resistive form of MHD, leaving the Hall term as part of

future work.

Compared with plasmas modeled with ideal MHD, there is no ”frozen-in” condition for

the resistive MHD model. This condition can be explained as follows: if some particles of

the plasma are attached by a magnetic field line at one time, they will stay attached by a

magnetic field line at future times and the magnetic topology of plasma is preserved. In

other words, the magnetic field and the plasma move together once we have this condition.

The frozen-in condition is valid when the following conditions are satisfied:

The ideal Ohm’s law :

E +

u × B

c

= 0,

Faraday’s law :

= −∇ × E.

∂B
∂t

However, the resistive Ohm’s law is

u × B

c

E +

= ηJ,

(5.1.3)

(5.1.4)

(5.1.5)

As will be derived in the next section, under equation 5.1.5, the filed equation becomes which

doesn’t satisfy the frozen-in condition.

= ∇ × (u × B) + cη∇2B,

∂B
∂t

(5.1.6)

where cη∇2B is the diffusion term. Thus the magnetic topology of the resistive MHD is

changing with the effect of resistivity since there is no ”frozen-in” condition anymore and it

88

leads to a ”reconnection” of different magnetic field lines, which is an important effect in the

macroscopic processes of MHD plasmas. Magnetic field lines reconnection happens in a very

thin layer in which the resistivity effect cannot be neglected since the parallel electric field

locally vanishes. The main goal of this work will be analyzing this important phenomenon.

Our approach to solve the resistive MHD system is through the magnetic potential equa-

tion derived from the magnetic induction equation for resistive MHD. The derivation of the

magnetic potential equation is given in the next section.

5.2 Magnetic Vector Potential Formulation for Resis-

tive MHD

In this section we will derive the magnetic potential equation for the resistive MHD model.

The resulting equations will be the dominant part for our constrained transport method for

resistive MHD.

Let’s write down Generalized Ohm’s Law for MHD [12, 60],

u × B

c

E +

= ηJ +

J × B
enec

− ∇ · Pe

neec

+

me
nee2

dJ
dt

,

(5.2.1)

where η is the collision frequency. In deriving the resistive model, we shall keep only the

first term on the right hand side of the equation (5.2.1), i.e.,

u × B

c

E +

= ηJ.

(5.2.2)

89

Let’s recall the Maxwell’s equations

with involutions

We also know that

1

c2 ∂tE = −η0J + ∇ × B,

∂tB = −∇ × E,

∇ · B = 0,

∇ · E = − ρ

0

.

B = ∇ × A,

1
c2 ∂tE = 0

in MHD. Then, using (5.2.4) and (5.2.7) we obtain

∇ × ∂tA = −∇ × E

∇ ×î

⇒

∂tA = −E

ó

.

The relation (5.2.10) implies the existence of a scalar function ψ such that

∂tA = −E + ∇ψ.

(5.2.3)

(5.2.4)

(5.2.5)

(5.2.6)

(5.2.7)

(5.2.8)

(5.2.9)

(5.2.10)

(5.2.11)

Moreover, if we apply the curl operator to the final term of the (5.2.11), we will get ∇×∇ψ =

0. Note that the function ψ is a gauge function, which is a degree of the freedom. Here, we

90

make use of the Weyl gauge, ψ = 0 [42]. Hence (5.2.11) becomes

∂tA = −E.

Using the equation (5.2.2) into (5.2.13), we get

∂tA = −î−u × B

ó

,

+ ηJ

c

Because of the (5.2.8) and the (5.2.3), we get

Using (5.2.7) into (5.2.14), we get

J =

∇ × B.

1
η0

J =

1
η0

∇ × (∇ × A).

(5.2.12)

(5.2.13)

(5.2.14)

(5.2.15)

Using the Grassmann identity, a × (b × c) = (a · c)b − (a · b)c = b(a · c) − (a · b)c, we obtain

î∇(∇ · A) − (∇ · ∇)A

ó

J =

1
η0

.

(5.2.16)

Equation (5.2.16) can be rewritten as



J =

1
η0

∂xyA2 + ∂xzA3 − ∂yyA1 − ∂zzA1
∂xyA1 + ∂yzA3 − ∂xxA2 − ∂zzA2
∂xzA1 + ∂yzA2 − ∂xxA3 − ∂yyA3

 .

(5.2.17)

Now we need to get u × B to complete the equation (5.2.13). So,

91

ï
ï
ï

η0

ò − η
ò − η
ò − η

η0

η0

ò
ò
ò

Ä
Ä

Ä

u × (∇ × A) =

u2(∂xA2 − ∂yA1) − u3(∂zA1 − ∂xA3)
u3(∂yA3 − ∂zA2) − u1(∂xA2 − ∂yA1)
u1(∂zA1 − ∂xA3) − u2(∂yA3 − ∂zA2)

ä
ä
ä



(5.2.18)

Then combining (5.2.17) and (5.2.18) will give the equation (5.2.13) in component-wise

ï
ï
ï

as

∂A1
∂t

=

∂A2
∂t

=

∂A3
∂t

=

1
c

1
c

1
c

u2∂xA2 + u3∂xA3 − u2∂yA1 − u3∂zA1

∂xyA2 + ∂xzA3 − ∂yyA1 − ∂zzA1

u3∂yA3 + u1∂yA1 − u3∂zA2 − u1∂xA2

u1∂zA1 + u2∂zA2 − u1∂xA3 − u2∂yA3

(5.2.19)
∂xyA1 + ∂yzA3 − ∂xxA2 − ∂zzA2

(5.2.20)
∂xzA1 + ∂yzA1 − ∂xxA3 − ∂yyA3

The resulting system for magnetic vector potential equation becomes a system of Hamilton-

Jacobi equations with diffusion terms since the velocity variables are given from the previous

time step. The approximation used for the Hamilton-Jacobi terms was studied in Chapter-3,

in which we applied the MOLT. In the next chapter, we discuss the approximations used for

(5.2.21)

the second order derivative terms. Note that in 2D, A = (0, 0, A3), so the system (5.2.19) -

(5.2.21) reduces to a single Hamilton-Jacobi equation with diffusion equation given by

ï − u1∂xA3 − u2∂yA3

ò − η

ï − ∂xxA3 − ∂yyA3

ò

η0

∂A3
∂t

=

1
c

(5.2.22)

In this work we focus on the 2D resistive MHD system, so we will use the scalar equation

92

for the magnetic potential.

5.3 Outline for Constrained Transport Framework for

Non-ideal MHD

In this section, we provide an outline of the steps used to solve the 2D non-ideal MHD system

using our constrained transport methodology. The basic idea of constrained transport (CT)

for non-ideal MHD is similar to the one used for ideal MHD. The approach begins by updating
the conserved variables q = (ρ, ρu,E, B) using a standard HCL-WENO scheme. Then, the

updated magnetic field is replaced with a discrete curl of the magnetic potential. The

importance of the last step is that it automatically maintains the divergence free condition

on the magnetic field. Here, we provide the important steps of the CT methodology to

advance the solution from current time step t = tn to the new time step tn+1.

Because we are using explicit time stepping, the values in the update only depend on

the old data in each of the stages of an RK method. We start by considering a single

Euler step. The values of B at level tn are all that are needed for the update of the energy

and momentum. Hence, our proposed method is to update the variables using ideal MHD,

followed by correcting B at time tn+1 with the non-ideal solution to A at tn+1. This avoids

the still CFL restrictions of having a diffusion term on the right hand side of B when doing

the update. This has to be done in the stage values in RK.

Consider a semi-discrete system of ordinary differential equations for ideal MHD equa-

tions (2.1.1)

q(cid:48)
mhd(t) = L(qmhd(t))

(5.3.1)

93

where qmhd(t) represents the grid function at time t consisting of all point-wise values of
the conserved quantities in the ideal MHD system qmhd = (ρ, ρu,E, B). The details about
L(qmhd(t)) were presented in [29]. Next, note that the evolution equation for the magnetic

potential in 2D is of the form

q(cid:48)
A(t) = H(qA(t), u(t)) + G(qA(t)),

(5.3.2)

where qA denotes the magnetic scalar potential, A3, in 2D, and G is the diffusion part of
the equation which involves second order derivative terms of the A3. In order to update this

equation, we proceed as follows:

• step 0: Start with the given current time step qn

mhd = (ρn, ρun,E n, Bn)T and qn
A.

• step 1: Apply a standard HCL - WENO method designed for the ideal MHD and

obtain q∗

mhd, where

q∗
mhd =

Ä
ρn+1, ρn+1un+1,E n+1, B∗ä

.

Here, B∗ is given with a ∗ superscript instead of n + 1 to indicate that the predicted B.

Note that since explicit time stepping depends on old data, the only wrong variable is
B∗ and it will be corrected using our CT method at each stage of Runge Kutta method

before the end of the time step.

• step 2: Solve the magnetic potential equation with diffusion qn+1

A

by applying our

proposed MOLT scheme for the Hamilton-Jacobi with diffusion equation (5.3.2).

94

• step 3: Replace B∗ to Bn+1 by a discrete curl of qn+1
A :

Bn+1 = ∇ × qn+1
A .

Note that we do not provide two options for energy density E as we saw for in the ideal

MHD case. Instead, we choose to conserve the total energy of the system. Also, note that as

we explained in the section 3.2, we satisfy the divergence free condition automatically, since

the divergence of the curl is identically zero.

In the next section, we will explain how we apply a MOLT scheme to approximate the

second derivatives appearing in the magnetic potential equation.

95

Chapter 6

Kernel Based Approximation with

Constrained Transport for Resistive

MHD

The magnetic potential equation that we have derived for resistive MHD in (5.2) is a Hamil-

ton Jacobi equation with diffusion terms. In this chapter, we will develop a kernel based

approximation for that HJ with diffusion equation [22], which now involves second order

spatial derivative terms different than the magnetic potential equation for ideal MHD. As in

the MOLT framework for the first spatial derivative term we developed in the chapter (3),

we also transform second order spatial derivative terms into a kernel based representation.

Such representations are also based on successive convolutions. The proposed method in

this chapter is a key part of the constrained transport methodology to update the magnetic

field solutions in our hybrid resistive MHD solver. We test the performance of the proposed

scheme on some benchmark problems in 2D and show some results for different diffusion

constant numbers η.

96

6.1 Hamilton Jacobi with Diffusion terms

The magnetic potential equation for the resistive-MHD system is considered as a Hamilton-

Jacobi equation with a diffusion term

At + H(Ax) = G(A)xx.

(6.1.1)

The main difference with the ideal case is the appearance of second derivative terms in the

equation for the magnetic vector potential. As with the first derivative terms, we transform

the second order spatial derivative terms into a kernel based representation. The scheme we

propose in this section also keeps the unconditional stability property. For the idea we use

in this chapter, we utilize the works [24, 17] and combined it with our HJ solver.

6.2 Approximation of the second order derivative ∂xx

In this subsection, we review the approximations used to form second order derivatives,
following the ideas of [24]. Let’s define L0 operator

L0 = I − 1

α2 ∂xx,

x ∈ [a, b],

(6.2.1)

where the I is the identity operator and α > 0 is a constant. For a given function v(x),

assume w(x) is the solution of the following differential equation

L0[w, α](x) = w(x) − 1

α2 wxx(x) = v(x).

(6.2.2)

97

and w(x) can be found by analytically constructing the inverse of the operator L0 as

w(x) = L−1

0 [v, α](x) = I0[v, α](x) + A0e−α(x−a) + B0e−α(b−x),

(6.2.3)

where

I0[v, α](x) :=

(cid:90) b

a

α
2

e−α|x−y|v(y)dy,

(6.2.4)

and A0 and B0 are constant values calculated by boundary conditions. For instance, for

periodic boundary conditions we have w(a) = w(b) and wx(a) = wx(b), then we obtain

A0 =

I0[v, α](b)

1 − µ

and B0 =

I0[v, α](a)

1 − µ

,

(6.2.5)

with µ = e−α(b−a).

Before we show how we approach second order derivative terms, we present how we

evaluate the convolution integral I0[v, α0](x). Note that, as stated in section (3.1.2.3), the

convolution integrals IL

i = IL[v, α](xi) and IR

i = IR[v, α](xi) satisfy a recursive relation

i = e−α∆xiIL
IL
i = e−α∆xi+1IR
IR

i−1 + J L
i ,

i+1 + J R
i ,

i = 1, . . . , N,

IL
0 = 0,
i = 0, . . . , N − 1,

IR
N = 0,

(6.2.6a)

(6.2.6b)

where,

(cid:90) xi

xi−1

J L
i = α

v(y)e−α(xi−y)dy,

J R
i = α

(cid:90) xi+1

xi

v(y)e−α(y−xi)dy.

(6.2.7)

Then the convolution integral I0[v, α0](x) can be broken into a left-biased integral IL[v, α0](x)

98

and right-biased integral IR[v, α0](x) such that

I0[v, α0](x) =

1
2

(IL[v, α0](x) + IR[v, α0](x)).

(6.2.8)

Hence, I0[v, α0](x) can be calculated based on left-biased and right-biased integrals, which
leads to O(N ) evaluation, see [22, 24].
Now we define a new operator D0 as

D0 = I − L−1
0 ,

(6.2.9)

and find that L0 = (I −D0)−1. By using the definition (6.2.1), we can construct the second
derivative term with respect to infinite series of the D0 as

1

α2 ∂xx = I − L0 = L0(L−1

0 − I) = −D0(I − D−1

0 ) = −

where operators are described successively as

0 = D0[Dp−1
Dp

0

].

∞(cid:88)

p=1

Dp
0,

(6.2.10)

(6.2.11)

Therefore, the term at the right hand side of the equation (6.1.1) can be approximated as

G(A)xx = −α2

∞(cid:88)

p=1

Dp
0[G(A), α](x).

(6.2.12)

Remember that we have already seen the approximation to the ∂x in the chapter 3. With that
in mind, we can now approximate to the −H(Ax) + G(A)xx with the following combination

99

of infinite series,

Ç

ˆH

∞(cid:88)

p=1

− αL

å

∞(cid:88)

p=1

− α2

0

∞(cid:88)

p=1

R[A−, αR]
Dp

Dp
0[G(A), α0]

(6.2.13)

where ˆH is numerical Hamiltonian. Here, α0, αL, and αR are chosen as

Dp
L[A+, αL], + αR

 

α0 =

β
b∆t

,

(cid:48)
|G

(A)|,

b = max

A

αL = αR =

β
c∆t

,

(cid:48)

|H

(A)|.,

c = max

A

where ∆t is time step and β is a prescribed constant. We need to truncate the series
In particular, −H(Ax) + G(A)xx is

and calculate only the corresponding partial sums.

approximated by kth partial sums as

−H(Ax) + G(A)xx ≈ ˆH

β
c∆t

](x), +

β
c∆t

k(cid:88)

p=1

R[A−,
Dp

β
c∆t

k(cid:88)

p=1

Dp
L[A+,

 

− β
c∆t
k(cid:88)

p=1

− β
b∆t

Dp
0[G(A),

β
b∆t

](x)

with accuracy O(∆tk). As explained in [22, 24], we have an extra term for the advection part

of the Hamilton Jacobi equation, which achieves fourth order accuracy for the case b << c,

i.e., convection dominates. Hence,

−H(Ax) + G(A)xx ≈ ˆH

β
c∆t

](x), +

β
c∆t

k(cid:88)

p=1

R[A−,
Dp

β
c∆t

å

å

](x)

(6.2.14)

](x)

(6.2.15)

Ç

Ç

k(cid:88)

p=1

Dp
L[A+,

 

− β
c∆t
k(cid:88)

p=1

− β
b∆t

Dp
0[G(A),

β
b∆t

](x)

100

Ç

ˆH

ñ

D2
L[A+,

+

D0

β
c∆t

],−D2

R[A−,

β
c∆t

β
c∆t

],

β
c∆t

ô

å

(x)

6.2.0.1 Periodic boundary conditions

For periodic boundary conditions, we require

A(a, t) = A(b, t), Ax(a, t) = Ax(b, t),

t ≥ 0,

(6.2.16)

the convolution operators are required to satisfy

Dp
0[G(A), α0](a) = Dp

0[G(A), α0](b),

L[H+(A), αL](a) = Dp
Dp
R[H−(A), αR](a) = Dp
Dp

L[H+(A), αL](b),
R[H−(A), αR](b).

(6.2.17)

For periodic boundary conditions, the coefficients A0 and B0 are prescribed by equation

(6.2.5). Recall that the coefficients AL and BR for periodic boundary conditions were pro-

vided in section (3.1.2).

6.3 Numerical Results for Resistive MHD for Periodic

Boundary Condition

In this section we will show some numerical results using resistive MHD equations for our

proposed scheme. The results shown here are to test our method using periodic boundary

conditions. We use third order SSP RK method for time discretization. Time step is chosen

101

as

∆t =

CFL
em ∗ cd

(6.3.1)

where the CFL number is 1, em = max(λ5, λ6), cd = 1

∆x + 1

∆y for 2D. This is CFL condition

of hyperbolic conservation laws, where the magnetic potential equation A with diffusion

is solved by kernel based method which is unconditionally stable. We choose β = 0.8 for

αL and αR whereas β = 0.4 for α0.

In the numerical experiments here, we have used

kernel based method for k=3. As stated in section (6.2), the approximation to magnetic

potential equation is expected to be fourth order accurate. In this section, we compare our

MHD solver based on MOLT with the WENO-FD solver, which was originally designed for

Hamilton Jacobi equations, but we add a second order explicit in time approximation for

diffusion term. Note that, the WENO-FD scheme is designed based on explicit time stepping

and has a restricted CFL condition for diffusion terms by von Neumann stability analysis.

Thus, for finite difference with explicit diffusion we use ∆t < min(

1

em∗cd , (∆x)2

2b

), where b

is the diffusivity and em is the maximum wave speed coming from the eigenvalues. We are

comparing the CFL restriction of the advection part and the diffusion part and make sure we

take the smallest time step. However, our MOLT solver is capable to take larger time steps

comparing to WENO-FD since our proposed HJ based on MOLT with diffusion behaves

unconditionally stable.

6.3.1 Refinement Study for Smooth Vortex Problem

We performed a successive spatial refinement study for the smooth vortex problem with
diffusion with 100× 100 at final time t = 0.03 for the diffusion η = 0.1. For refinement study

102

we choose CFL=1 and β = 0.6 for advection while β = 0.4 for diffusion part.

In Table (6.1), we provide a table of error results, which demonstrate that the method

achieves high order accuracy (almost 5th order of accuracy) in space since 5th order WENO-

quadrature is used.

Table 6.1: Smooth vortex problem. Errors of A and orders of accuracy.

Nx × Ny
32 × 32
64 × 64
128 × 128
256 × 256

L1 error

1.9162E-04
1.1876E-05
4.7350E-07
1.6577E-08

order

–

4.01
4.65
4.84

L∞ error
6.6734E-02
3.9467E-03
2.6466E-04
1.0801E-05

order

–

4.08
3.90
4.61

6.3.2 Smooth Vortex Problem for Resistive MHD

In Figures 6.1 and 6.2 we show the results at the final time t = 20 using a 100 × 100

mesh and CFL=. At time t = 20, the solution has advanced by one period. The results

are for different diffusion coefficients η = 0, 0.001, 0.1, 0.5, 1. We compare the FD-WENO

scheme with our MOLT scheme. While the FD-WENO scheme can only give results up to

η = 0.5, our method works well, even for η = 1. Moreover, the ability of taking CFL=1

for MOLT is computationally efficient. For example, for the case of η = 0.5 our solver

requires %23 less time than the FD with diffusion method due to its CFL restriction ∆t <

min(

2∗em∗cd , (∆x)2

2b

1

). Note that once we refine the problem in space, (∆x)2

2b will dominate

the CFL selection and we will get more time steps, thus we will get square more efficiency

comparing to explicit scheme. This is a simple benchmark problem and there are better

ways to solve the diffusion terms.

103

(a) η = 0

(b) η = 0.01

(c) η = 0.1

(d) η = 0.5

(e) η = 1

Figure 6.1: (FD) Contour plots of (cid:107)B(cid:107) at time t = 20 with ∆x = ∆y = 0.2.

104

-10-50510-10-505100.10.20.30.40.50.60.70.8-10-50510-10-505100.10.20.30.40.5-10-50510-10-505100.010.020.030.040.050.060.070.08-10-50510-10-505101234567810-3-10-50510-10-50510(a) η = 0

(b) η = 0.01

(c) η = 0.1

(d) η = 0.5

(e) η = 1

Figure 6.2: (Kernel based method) Contour plots of (cid:107)B(cid:107) at time t = 20 with ∆x = ∆y = 0.2.

105

-10-50510-10-505100.10.20.30.40.50.60.70.8-10-50510-10-505100.10.20.30.40.5-10-50510-10-505100.010.020.030.040.050.060.070.08-10-50510-10-505101234567810-3-10-50510-10-505100.511.522.510-36.4 Orszang Tang Vortex problem for Resistive MHD

In this section we conduct a refinement study and show 2D contour plots of the density

for the Orzsag Tang Vortex problem with a resistive MHD model. The Orszag Tang Vortex

problem is a good benchmark for numerical schemes, as the solution is sensitive to divergence

errors in the magnetic field. We have two-way periodic boundary conditions on a domain
[0, 2π] × [0, 2π]. We use 192 × 192 for mesh size, and β = 0.8 for advection part whereas

β = 0.2 for the diffusion part. We are able to take a CFL of 1 for MOLT, while CFL relation

is ∆t < min(

1

em∗cd, (∆x)2

2b

) for WENO-FD due to the stability condition for the diffusion

terms for the explicit time stepping. With this CFL, MOLT requires less time to run the

problem.

6.4.1 Refinement Study for Orszag Tang Vortex problem for Re-

sistive MHD

We show the successive refinement study in this section. The mesh size is 192 × 192 at

time t = 0.03. we choose β = 0.8 for advection part whereas β = 0.4 for the diffusion

part and use CFL=0.5 for refinement study. In table (6.2) we obtain 5th order accuracy in

space for A and in table (6.3) we show the time refinement and observe 4th order in time,

which thereby verifying the analysis presented in section (6.2). Note that we use 5th order

WENO-based quadrature for space discretization and observe 5th order accuracy in space

as expected based on the analysis presented in section (6.2), which confirms the theory.

106

Nx × Ny
16 × 16
32 × 32
64 × 64
128 × 128

L1 error

Table 6.2: Errors of A and orders of accuracy.
L∞ error
9.4445E-03
7.7115E-04
3.1896E-05
8.3855E-07

1.7425E-03
1.0593E-04
3.1380E-06
6.7905E-08

4.04
5.08
5.53

order

–

order

–

3.61
4.60
5.25

Table 6.3: Errors of A and orders of temporal accuracy.

CF L

0.5
0.25

L1 error

2.4474E-03
1.5242E-04

order

–

4.01

L∞ error
1.1966E-03
6.6924E-05

order

–

4.16

6.4.2 Orszag Tang Vortex problem with diffusion for Resistive

MHD

We add diffusion terms for the magnetic potential for this problem and we want to see the

evolution of the problem for the different diffusion constants η = 0, 0.001, 0.1 and 1.

In

Figures 6.3 and 6.4 we compare the solutions obtained from a FD scheme to those of our
proposed MOLT scheme, at the final time t = 3 with a 192 × 192 mesh. We are able to
em∗cd, (∆x)2

use CFL of one for MOLT, while we can use the CFL relation ∆t < min(

2b

)

1

for WENO-FD. With the ability of taking larger time steps, we need less time to run the

problems. For example, for the case of η = 0.01 MOLT code requires %40 less time comparing

to FD-WENO code. Note that as stated in the previous example, once we refine the problem

in space, then the diffusion restriction of CFL will dominate the CFL selection for Finite

Differencing method and cause more time steps, thus we will get square more efficiency

comparing the WENO-FD method.

107

(a) η = 0

(b) η = 0.01

(c) η = 0.1

(d) η = 1

Figure 6.3: (FD) Contour plots of ρ at time t = 3 with 192 × 192.

108

(a) η = 0

(b) η = 0.01

(c) η = 0.1

(d) η = 1

Figure 6.4: (MOLT) Contour plots of ρ at time t = 20 with 192 × 192.

109

Chapter 7

MOLT for PEC for non-ideal MHD

Equations

In this chapter we will expand our resistive MHD solver to work for perfectly electrically

conducting (PEC) boundary condition. In the chapter 4, we have presented the PEC bound-

ary conditions for ideal MHD case, which we only had first order spatial derivative terms in

the magnetic potential equation. We have utilized the Cauchy-Kowalevski theorem there to

derive boundary terms AL and BR. However, it is not that easy to use Cauchy-Kowalevski

theorem for the second order derivative cases. Thus, we are using the method for diffusion

that is consistent with our approach that is only 1st order at the boundary. Note that we

are not correcting at the boundary and leaving that for future work to achieve higher or-

der accuracy. Then we will demonstrate the results on several test problems such as the

Hartmann Flow and the GEM challenge problem.

7.1 MOLT for PEC for non-ideal MHD

In this section we examine the operator D0 for non periodic boundary condition, which we

use for the approximation of the second order derivative terms.

D0[v, α](x) = v(x) − α
2

(cid:90) b

a

e−α|x−y|v(y)dy − A0e−α(x−a) − B0e−α(b−x).

(7.1.1)

110

where A0 and B0 are determined by the boundary conditions. Following the idea in [23], if

we require non periodic boundary conditions, we assume that we are given numbers Ca and

Cb. If we require

D0[v, α](a) = Ca,

and D0[v, α](b) = Cb,

then, the boundary coefficients are obtained as

A0 =

B0 =

1

1 − µ2

1

1 − µ2

I0[v, α](b) − v(b) + Cb
I0[v, α](a) − v(a) + Ca

µ

µ

Ä
Ä

Ä
Ä

ä −Ä
ä −Ä

I0[v, α](a) − v(a) + Ca

I0[v, α](b) − v(b) + Cb

ää
ää

,

.

(7.1.2a)

(7.1.2b)

In the next section we will test our code on several test problems such as the Hartmann

Flow problem and the GEM Magnetic Reconnection problem.

7.2 Numerical results for Hartmann Flow problem with

diffusion

The kernel based approximation based on the MOLT for PEC boundary conditions for the

second order derivative term will be tested on the Hartmann Flow which is derived in the

section (4.2). This benchmark problem is a good test problem to test the code both for ideal

and non-ideal MHD models. Here, we will let the diffusion variable η be nonzero and we will

test our code for the diffusion effect.

We choose β = 0.8 for αL and αR whereas β = 0.6 for α0. Since the Hartmann problem

assumes the MHD model to be incompressible, we use ρ = 1 and p = 1. We show refinement

study in the table (7.1) at time t = 0.1 for diffusion η = 0.1 and we choose CFL=0.5 for

111

refinement study. We get low order of accuracy in space since we only use first order accurate

boundary terms A0 and B0 for PEC boundary conditions for diffusion, while we use higher

order approximations in the domain. In this test problem we observe that we get first order

accuracy for magnetic potential equation as expected by the theory. Note that we intend

to achieve higher order approximations for PEC boundary conditions for diffusion in the
future by correcting boundary terms A0 and B0. For resolution 128 × 128, final time t = 1,

Table 7.1: Hartmann flow problem. Errors of A and orders of accuracy.

Nx × Ny
32 × 32
64 × 64
128 × 128
256 × 256

L1 error

3.1362E-07
1.8958E-07
1.0127E-07
5.2168E-08

order

–

0.73
0.90
0.96

L∞ error
3.0554E-06
2.1523E-06
1.2318E-06
6.5960E-07

order

–

0.51
0.81
0.90

CFL=1 and η = 0.001, the velocity and the induced magnetic field solutions for the different

Hartmann numbers are shown in figures 7.1 and 7.3. In Figures 7.2 and 7.4, we presented

1D cross sections of velocity and induced magnetic field for different Hartmann numbers

Ha = 0, 2, 5, 10. The numerical results overlap on the exact solutions.

7.3 Magnetic reconnection: GEM challenge problem

The main motivation of going beyond the ideal MHD by including additional physics in the

plasma model is to observe the magnetic reconnection phenomenon, which is a significant

cosmic plasma process. This phenomenon has an effect on the dynamics of the magneto-

sphere by facilitating the entry of particles and energy from the Sun’s atmosphere into the

Earth’s magnetosphere at the magnetopause by the solar wind. When the high energy parti-

112

(a) Muller [58]

(b) Molt

Figure 7.1: Velocity profile for the Hartmann problem at Ha = 0, 2, 5, 10.

cles reach Earth and magnetic reconnection happens, then this restricts the performance of

fusion reactors and also affects the geospace weather which can damage power grids, commu-

nication systems, GPS navigation and threaten spacecraft and satellites. One other effect

of reconnection is believed to be helping trigger aurora, also known as the northern and

southern lights. What the magnetic reconnection really is the breaking and reconnecting of

magnetic field lines in a highly conducting plasma. The magnetic field lines get squeezed

together as magnetized plasma flows are going toward each other. The particles resist being

squeezed together and the field lines are broken and the reconnection happens in the tiny

diffusion region.

In plasma physics, it is known that for magnetic reconnection to occur, a mechanism is

required to break the ”frozen-in” flux condition, which could be resistivity, electron inertia or

anisotropic electron pressure, which appear in the generalized Ohm’s law. This phenomenon

does not occur in the ideal MHD plasma (infinitely conductive plasma). In this work, it is

113

(a) Ha = 10

(b) Ha = 5

(c) Ha = 2

(d) Ha = 0

Figure 7.2: Exact and Numerical velocity solutions for the Hartmann problem at Ha =
0, 2, 5, 10.

the resistivity that breaks this condition since we work within the resistive MHD model.

The Geospace Environment Modeling (GEM) Magnetic Reconnection Challenge Problem

[9] is one of the benchmark problems of magnetic reconnection and it is studied in this work.

The study of Magnetic Reconnection is on a simple 2D current-sheet, called the Harris sheet

configuration, with a specified set of initial conditions and magnetic island perturbation to

trigger the dynamics. The Harris sheet configuration is a planar MHD equilibrium with a

horizontal magnetic field given by

B1(y) = B0 tanh(y/λB).

(7.3.1)

114

(a) Muller [58]

(b) Molt

Figure 7.3: Induced magnetic field for the Hartmann problem at Ha = 0, 2, 5, 10.

Assume a uniform temperature T = T0 and a density profile given by

ρ(y) = ρ0 cosh−2(y/λB) + ρ∞.

The initial equilibrium is perturbed by an additional magnetic perturbation

δB1(x, y) = −ψ(π/Ly) cos(2πx/Lx) sin(πy/Ly)

δB2(x, y) = ψ(2π/Lx) sin(2πx/Lx) cos(πy/Ly).

(7.3.2)

(7.3.3)

The standard case takes − 1

2 Lx < x < 1

2 Lx and − 1

2 Ly < y < 1

2 Ly with Lx = 25.6 and

Ly = 12.8. The simulations use periodic boundary conditions horizontally, while perfectly
conducting impenetrable walls are assumed at y = ± 1

2 Ly. A perfectly conducting wall is a

solid wall boundary and can be given as follows

115

(a) Ha = 10

(b) Ha = 5

(c) Ha = 2

(d) Ha = 0

Figure 7.4: Exact and Numerical induced magnetic field solutions for the Hartmann problem
at Ha = 0, 2, 5, 10.

∂yρ = 0,

∂yBx = By = ∂yBz = 0

∂yux = uy = ∂yuz = 0,

Ex = ∂yEy = Ez = 0

(7.3.4)

(7.3.5)

In Figure (7.7), we have presented our kernel based method results for the GEM challenge

problem and compared it with the work done by Birn et. al., [9]. The mesh spacing is
128 × 128 with CFL = 1. We use β = 0.8 for the advection part of the equation while

β = 0.6 for the diffusion part. The plots show the evolution of the magnetic field lines and

current density for the resistive MHD for uniform resistivity η = 0.001.The formation of the

116

Figure 7.5: Magnetic Reconnection in Space. Credits: NASA

magnetic island in the center of the current sheet matches the work done by the Birn et. al.,

[9]. The magnetic island is increasing over the time and then merging with the bigger island

to the left.

117

Figure 7.6: Breaking and rejoining of magnetic field lines due to the presence of localized
diffusion region (shaded). ( Priest, E. R [63])

(a) Birn et. al., [9]

(b) Kernel-based method

Figure 7.7: GEM Magnetic Reconnection Problem for diffusion η = 0.001.

118

-10-50510(a) t=50-6-4-2024600.511.522.53-10-50510(b) t=100-6-4-20246-10-50510(c) t=150-6-4-20246Chapter 8

Conclusion and Future Work

8.1 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 [29]. 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

119

and blast wave problems.

In Chapter 4, we developed the PEC boundary conditions for magnetic potential equation

in MOLT-HJ scheme. In fact, the PEC boundary condition is equivalent enforcing Dirichlet

boundary conditions. We derived the MOLT-HJ scheme for Dirichlet boundary conditions

for magnetic potential equation. Then we go over the derivation of the Hartmann Flow

problem, which is a benchmark problem to test PEC boundary conditions both for ideal and

non ideal MHD equations. We get high order accuracy for ideal Hartmann problem. We also

showed that the exact solution and the numerical solutions overlap each other by showing

1D cross sections and 2D contour plots.

In Chapter 5, we further our work from ideal MHD to resistive Magnetohydrodynamics

model. We derived the magnetic potential equation from the magnetic field equation which

has the diffusion term in it. We obtained Hamilton-Jacobi with diffusion type of equation

for magnetic potential equation, which still satisfies the divergence free condition for the

magnetic field. We had a kernel based approximation for the first derivative terms in the

ideal MHD in Chapter 3, then in the Chapter 6 we developed a similar approximation for

the second order derivative terms existed in the magnetic potential equation for the resistive

MHD. We presented MOLT-HJ method for the periodic boundary conditions. We tested our

code for periodic boundary conditions on the several 2D test problems. We get high order

accuracy on these examples, Smooth Vortex and Orszag Tang Vortex problems. We also

presented 2D plots of these problems both for FD-WENO code and our MOLT-HJ code for

different diffusion constants η. We observe that MOLT-HJ is working well and requires less

time since our kernel based approximation is unconditionally stable and allows taking larger

time steps.

In Chapter 7, we introduced the PEC boundary conditions for the MOLT method for

120

non-ideal case of MHD equations. We tested our work on both Hartmann Flow problem and

GEM Magnetic Reconnection Problems. For the resistive flow example, we get first order

of accuracy for the Hartmann problem and left the study of making it higher order for the

future.

8.2 Future Work

The main objective of resistive MHD work is to solve the Magnetic Reconnection problem

for the PEC boundary condition at the top and bottom and periodic boundary conditions at

the left and right boundaries. The non-ideal evaluation of code gave us low order of accuracy

so far. For such an important problem we will endeavor to obtain a higher order of accuracy

in future work.

We also will seek to further our resistive MHD work to the Hall MHD equations by

including more physics in Ohm’s law. That problem is more sensitive to diffusions and is

challenging to solve. Keeping the Hall terms in the equations will not change the structure

of the magnetic potential equation, thus we will be able to use our Constrained Transpose

methodology for resistive MHD model to the Hall MHD model.

121

APPENDIX

122

Appendix

WENO-Integration Formulation

Here, we show all coefficients in WENO quadrature. We denote γ = α∆x, then we have

6 − 6ν + 2ν2 − (6 − ν2)e−ν

6ν3

1 = − 6 − 8ν + 3ν2 − (6 − 2ν − 2ν2)e−ν

(0)

,

6 − 10ν + 6ν2 − (6 − 4ν − ν2 + 2ν3)e−ν

,

2ν3

2ν3

3 = − 6 − 12ν + 11ν2 − 6ν3 − (6 − 6ν + 2ν2)e−ν

(0)

,

6ν3
6 − ν2 − (6 + 6ν + 2ν2)e−ν

6ν3

1 = − 6 − 2ν − 2ν2 − (6 + 4ν − ν2 − 2ν3)e−nu

(1)

6 − 4ν − ν2 + 2ν3 − (6 + 2ν − 2ν2)e−ν

2ν3

3 = − 6 − 6ν + 2ν2 − (6 − ν2)e−ν

(1)

,

6ν3

2ν3

(2)
0 =

6 + 6ν + 2ν2 − (6 + 12ν + 11ν2 + 6ν3)e−ν
1 = − 6 + 4ν − ν2 − 2ν3 − (6 + 10ν + 6ν2)e−ν

6ν3

(2)

,

2ν3

6 + 2ν − 2ν2 − (6 + 8ν + 3ν2)e−ν

,

,

,

,

(2)
2 =

3 = − 6 − ν2 − (6 + 6ν + 2ν2)e−ν

(2)

ui+2.

2ν3

6ν3

(0)
0 =

(0)
2 =

(1)
0 =

(1)
2 =

c

c

c

c

c

c

c

c

c

c

c

c

,

,

123

And the linear weights are

6 − ν2 − (6 + 6ν + 2ν2)e−ν
3ν((2 − ν) − (2 + ν)e−ν)
60 − 60ν + 15ν2 + 5ν3 − 3ν4 − (60 − 15ν2 + 2ν4)e−ν

10ν2(6 − ν2 − (6 + 6ν + 2ν2)e−ν)

d0 =

d2 =

d1 = 1 − d0 − d2.

The smoothness indicator βr have the expressions as

β0 =

β1 =

β2 =

13
12
13
12
13
12

(−vi−3 + 3vi−2 − 3vi−1 + vi)2 +
(−vi−2 + 3vi−1 − 3vi + vi+1)2 +
(−vi−1 + 3vi − 3vi+1 + vi+2)2 +

(vi−3 − 5vi−2 + 7vi−1 − 3vi)2,
(vi−2 − vi−1 − vi + vi+1)2,
(−3vi−1 + 7vi − 5vi+1 + vi+2)2.

1
4
1
4
1
4

124

BIBLIOGRAPHY

125

BIBLIOGRAPHY

[1] D. S. Balsara. Divergence-free adaptive mesh refinement for magnetohydrodynamics.

Journal of Computational Physics, 174(2):614–648, 2001.

[2] D. S. Balsara.

divergence-free reconstruction.
151(1):149, 2004.

Second-order-accurate schemes for magnetohydrodynamics with
The Astrophysical Journal Supplement Series,

[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. Spicer. Maintaining pressure positivity in magnetohydrodynamic

simulations. Journal of Computational Physics, 148(1):133–148, 1999.

[8] D. S. Balsara and D. S. Spicer. A staggered mesh algorithm using high order godunov
fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations. Jour-
nal of Computational Physics, 149(2):270–292, 1999.

[9] J. Birn and M. Hesse. Geospace environment modeling (gem) magnetic reconnection
challenge: Resistive tearing, anisotropic pressure and hall effects. Journal of Geophysical
Research: Space Physics, 106(A3):3737–3750, 2001.

[10] D. Biskamp. Magnetic reconnection via current sheets. The Physics of fluids,

29(5):1520–1531, 1986.

[11] D. Biskamp, E. Schwarz, and J. F. Drake. Two-fluid theory of collisionless magnetic

reconnection. Physics of Plasmas, 4(4):1002–1009, 1997.

[12] J. A. Bittencourt. Fundamentals of plasma physics. Springer Science & Business Media,

2013.

126

[13] R. Borges, M. Carmona, B. Costa, and W. S. Don. An improved weighted essen-
tially non-oscillatory scheme for hyperbolic conservation laws. Journal of Computational
Physics, 227(6):3191–3211, 2008.

[14] 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.

[15] 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.

[16] 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.

[17] 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.

[18] M. F. Causley and A. J. Christlieb. Higher order a-stable schemes for the wave equa-
tion using a successive convolution approach. SIAM Journal on Numerical Analysis,
52(1):220–235, 2014.

[19] 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.

[20] 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.

[21] 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.

[22] A. Christlieb, W. Guo, and Y. Jiang. Kernel based high order “explicit” unconditionally-
stable scheme for nonlinear degenerate advection-diffusion equations. arXiv preprint
arXiv:1707.09294, 2017.

[23] A. Christlieb, W. Guo, and Y. Jiang. A kernel based high order “explicit” uncondition-
ally stable scheme for time dependent hamilton–jacobi equations. Journal of Computa-
tional Physics, 379:214–236, 2019.

[24] A. Christlieb, W. Guo, Y. Jiang, and H. Yang. Kernel based high order ?explicit?
unconditionally stable scheme for nonlinear degenerate advection-diffusion equations.
Journal of Scientific Computing, 82(3):52, 2020.

127

[25] A. J. Christlieb, X. Feng, Y. Jiang, and Q. Tang. A high-order finite difference weno
scheme for ideal magnetohydrodynamics on curvilinear meshes. SIAM Journal on Sci-
entific Computing, 40(4):A2631–A2666, 2018.

[26] A. J. Christlieb, X. Feng, D. C. Seal, and Q. Tang. A high-order positivity-preserving
single-stage single-step method for the ideal magnetohydrodynamic equations. Journal
of Computational Physics, 316:218–242, 2016.

[27] A. J. Christlieb, Y. Liu, Q. Tang, and Z. Xu. High order parametrized maximum-
principle-preserving and positivity-preserving weno schemes on unstructured meshes.
Journal of Computational Physics, 281:334–351, 2015.

[28] A. J. Christlieb, Y. Liu, Q. Tang, and Z. Xu. Positivity-preserving finite difference
weighted eno schemes with constrained transport for ideal magnetohydrodynamic equa-
tions. SIAM Journal on Scientific Computing, 37(4):A1825–A1845, 2015.

[29] A. J. Christlieb, J. A. Rossmanith, and Q. Tang. Finite difference weighted essentially
non-oscillatory schemes with constrained transport for ideal magnetohydrodynamics.
Journal of Computational Physics, 268:302–325, 2014.

[30] B. Cockburn and C.-W. Shu. Tvb runge-kutta local projection discontinuous galerkin
finite element method for conservation laws. ii. general framework. Mathematics of
computation, 52(186):411–435, 1989.

[31] W. Dai and P. R. Woodward. A simple finite difference scheme for multidimensional
magnetohydrodynamical equations. Journal of Computational Physics, 142(2):331–369,
1998.

[32] 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.

[33] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz, T. Schnitzer, and M. Wesenberg. Hy-
perbolic divergence cleaning for the mhd equations. Journal of Computational Physics,
175(2):645–673, 2002.

[34] C. R. DeVore. Flux-corrected transport techniques for multidimensional compressible

magnetohydrodynamics. Journal of Computational Physics, 92(1):142–160, 1991.

[35] E. Dorfi. Numerical methods for astrophysical plasmas. Computer Physics Communi-

cations, 43(1):1–15, 1986.

[36] L. Dubuc, F. Cantariti, M. Woodgate, B. Gribben, K. Badcock, and B. Richards. Solu-
tion of the unsteady euler equations using an implicit dual-time method. AIAA journal,
36(8):1417–1424, 1998.

128

[37] C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows-a constrained

transport method. The Astrophysical Journal, 332:659–677, 1988.

[38] 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.

[39] 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.

[40] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time

discretization methods. SIAM review, 43(1):89–112, 2001.

[41] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy. Uniformly high order
accurate essentially non-oscillatory schemes, iii. In Upwind and high-resolution schemes,
pages 218–290. Springer, 1987.

[42] 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.

[43] 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.

[44] M. Hesse, K. Schindler, J. Birn, and M. Kuznetsova. The diffusion region in collisionless

magnetic reconnection. Physics of Plasmas, 6(5):1781–1795, 1999.

[45] G.-S. Jiang. Shu ch.-w.: Efficient implementation of weighted eno schemes. Journal of

Computational Physics, 126(1):202–228, 1996.

[46] G.-S. Jiang and D. Peng. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM

Journal on Scientific computing, 21(6):2126–2143, 2000.

[47] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes. Journal

of computational physics, 126(1):202–228, 1996.

[48] G.-S. Jiang and C.-c. Wu. A high-order weno finite difference scheme for the equations
of ideal magnetohydrodynamics. Journal of Computational Physics, 150(2):561–594,
1999.

129

[49] S. Kawai. Divergence-free-preserving high-order schemes for magnetohydrodynamics:
An artificial magnetic resistivity method. Journal of Computational Physics, 251:292–
318, 2013.

[50] G. Laval, R. Pellat, and M. Vuillemin. Electromagnetic instabilities in a collisionless
plasma;
instabilites electromagnetiques des plasmas sans collisions; ehlektromagnit-
nye neustojchivosti besstolknovitel’nykh plazm; inestabilidad electromagnetica en un
plasma sin colisiones. 1966.

[51] R. J. LeVeque. Wave propagation algorithms for multidimensional hyperbolic systems.

Journal of Computational Physics, 131(2):327–353, 1997.

[52] 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.

[53] X.-D. Liu, S. Osher, and T. Chan. Weighted essentially non-oscillatory schemes. Journal

of computational physics, 115(1):200–212, 1994.

[54] P. Londrillo and L. Del Zanna. High-order upwind schemes for multidimensional mag-

netohydrodynamics. The Astrophysical Journal, 530(1):508, 2000.

[55] 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.

[56] Z. Ma and A. Bhattacharjee. Fast impulsive reconnection and current sheet intensifica-
tion due to electron pressure gradients in semi-collisional plasmas. Geophysical research
letters, 23(13):1673–1676, 1996.

[57] M. Mandt, R. Denton, and J. Drake. Transition to whistler mediated magnetic recon-

nection. Geophysical Research Letters, 21(1):73–76, 1994.

[58] U. M¨uller and L. B¨uhler. Magnetofluiddynamics in channels and containers. Springer

Science & Business Media, 2013.

[59] G. K. Parks. Physics of space plasmas-an introduction. Redwood City, CA, Addison-

Wesley Publishing Co., 1991, 547 p., 1991.

[60] G. K. Parks. Physics of space plasmas: an introduction. CRC Press, 2019.

[61] K. G. Powell. An approximate riemann solver for magnetohydrodynamics (that works
more than one dimension). Technical report, INSTITUTE FOR COMPUTER APPLI-
CATIONS IN SCIENCE AND ENGINEERING HAMPTON VA, 1994.

130

[62] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw. A solution-
adaptive upwind scheme for ideal magnetohydrodynamics. Journal of Computational
Physics, 154(2):284–309, 1999.

[63] E. R. Priest. MHD reconnection. Scholarpedia, 6(2):2371, 2011. revision #136731.

[64] J. A. Rossmanith. An unstaggered, high-resolution constrained transport method for
magnetohydrodynamic flows. SIAM Journal on Scientific Computing, 28(5):1766–1797,
2006.

[65] D. Ryu, F. Miniati, T. Jones, and A. Frank. A divergence-free upwind code for multi-
dimensional magnetohydrodynamic flows. The Astrophysical Journal, 509(1):244, 1998.

[66] D. C. Seal, Q. Tang, Z. Xu, and A. J. Christlieb. An explicit high-order single-stage
single-step positivity-preserving finite difference weno method for the compressible euler
equations. Journal of Scientific Computing, 68(1):171–190, 2016.

[67] M. A. Shay and J. F. Drake. The role of electron dissipation on the rate of collisionless

magnetic reconnection. Geophysical Research Letters, 25(20):3759–3762, 1998.

[68] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock

capturing schemes, 2. 1988.

[69] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-
In Upwind and High-Resolution Schemes, pages 328–374.

ii.

capturing schemes,
Springer, 1989.

[70] M. Torrilhon. Locally divergence-preserving upwind finite volume schemes for magne-
tohydrodynamic equations. SIAM Journal on Scientific Computing, 26(4):1166–1191,
2005.

[71] G. T´oth. The ∇ · B = 0 constraint in shock-capturing magnetohydrodynamics codes.

Journal of Computational Physics, 161(2):605–652, 2000.

[72] B. Udrea. An advanced implicit solver for MHD. University of Washington, 1999.

[73] K. J. Vanden and D. Whitfield. Direct and iterative algorithms for the three-dimensional

euler equations. AIAA journal, 33(5):851–858, 1995.

[74] V. M. Vasyliunas. Theoretical models of magnetic field line merging. Reviews of Geo-

physics, 13(1):303–336, 1975.

[75] 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.

131

[76] T. Xiong, J.-M. Qiu, and Z. Xu. Parametrized positivity preserving flux limiters for the
high order finite difference weno scheme solving compressible euler equations. Journal
of Scientific Computing, 67(3):1066–1088, 2016.

[77] Z. Xu. Parametrized maximum principle preserving flux limiters for high order schemes
solving hyperbolic conservation laws: one-dimensional scalar problem. Mathematics of
Computation, 83(289):2213–2238, 2014.

[78] K. Yee. Numerical solution of initial boundary value problems involving maxwell’s equa-
tions in isotropic media. IEEE Transactions on antennas and propagation, 14(3):302–
307, 1966.

[79] A. L. Zachary, A. Malagoli, and P. Colella. A higher-order godunov method for mul-
tidimensional ideal magnetohydrodynamics. SIAM Journal on Scientific Computing,
15(2):263–284, 1994.

[80] X. Zhang and C.-W. Shu. On positivity-preserving high order discontinuous galerkin
schemes for compressible euler equations on rectangular meshes. Journal of Computa-
tional Physics, 229(23):8918–8934, 2010.

132