GAMERA: A Three-dimensional Finite-volume MHD Solver for Non-orthogonal Curvilinear Geometries

Efficient simulation of plasmas in various contexts often involves the use of meshes that conform to the intrinsic geometry of the system under consideration. We present here a description of a new magnetohydrodynamic code, GAMERA (Grid Agnostic MHD for Extended Research Applications), designed to combine geometric flexibility with high-order spatial reconstruction and constrained transport to maintain the divergence-free magnetic field. GAMERA carries on the legacy of its predecessor, the LFM (Lyon–Fedder–Mobarry), a research code whose use in space physics has spanned three decades. At the time of its initial development, the LFM code had a number of novel features: eighth-order centered spatial differencing, the Partial Donor Cell Method limiter for shock capturing, a non-orthogonal staggered mesh with constrained transport, and conservative averaging-reconstruction for axis singularities. The capability to handle multiple ion species was also added later. GAMERA preserves the core numerical philosophy of LFM while also incorporating numerous algorithmic and computational improvements. The upgrades in the numerical schemes include accurate grid metric calculations using high-order Gaussian quadrature techniques, high-order upwind reconstruction, non-clipping options for interface values, and improved treatment of axis singularities. The improvements in the code implementation include the use of data structures and memory access patterns conducive to aligned vector operations and the implementation of hybrid parallelism, using MPI and OMP. GAMERA is designed to be a portable and easy-to-use code that implements multidimensional MHD simulations in arbitrary non-orthogonal curvilinear geometries on modern supercomputer architectures.


Introduction
Numerical magnetohydrodynamics (MHD), which solves the equations of ideal or non-ideal MHD, is used extensively in research areas such as basic plasma physics, astrophysics, solar physics, geospace environment modeling, planetary sciences and space weather applications.Therefore the development of accurate, multi-dimensional, general-purposed MHD solvers along with comprehensive descriptions of their implementation is important for building powerful computational tools and also for the advancement of scientific understanding, and education of the new generation of scientists.In the past two decades, many general-purposed MHD codes have been developed by research teams from various scientific research communities including but not limited to, the ZEUS code [1], the Athena code [2], the BATS-R-US code [3], the VAC code [4], the Nirvana [5], the RAMSES [6], the PLUTO code [7], the AstroBEAR [8], the FLASH code [9] and the PENCIL code [10].One of the MHD codes that has been widely used within the space physics community is the Lyon-Fedder-Mobarry (LFM) code [11].
The LFM MHD code, initially developed at the Naval Research Laboratory in the early 80s, is one of the pioneers of solving three-dimensional MHD equations in non-orthogonal curvilinear geometry with high-quality advection schemes [11,12,13].The MHD Kernel of the LFM code had a number of novel features: an eighth-order centered spatial reconstruction method to reduce numerical diffusion, a universal-type non-linear limiter to capture shocks [14], a non-orthogonal spherical mesh with high-order constrained transport algorithm to maintain zero magnetic divergence at round-off, a fix for axis singularities and multi-fluid extensions to handle multiple ion species [15].The LFM code has been well known as a global terrestrial magnetosphere model, which has been one of the workhorses for geospace research and space weather applications [16].In the past decade, as the backbone of a whole geospace model, the LFM has been coupled to other first-principles codes, e.g. the magnetosphere-ionosphere coupler/solver code for high-latitude ionospheric electrodynamics [17], the Rice Convection Model of the inner magnetosphere [18], the NCAR thermosphere-ionosphere model for upper atmospheric dynamics [19] and the ionosphere polar wind model [20] for collisionless heavy ion transport.Beyond geospace research, the LFM code has also been adapted to be an inner heliosphere model [21], a global heliosphere model [22], an unmagnetized planetary magnetosphere model for Venus [23], a giant planetary magnetosphere model for Jupiter [24] and basic plasma simulations for magnetic reconnection [25].The code is actively used by research teams in the space physics community, and the list of applications developed based on the LFM code is still growing.
For a code to maintain its place in computational science for almost 40 years is a testament to the quality of its core numerics, but the age of the LFM code makes it difficult to adapt to modern and future architectures.To address these issues, we completed a rewrite of the LFM code from scratch, and re-envisaging it as GAMERA -Grid Agnostic MHD for Extended Research Applications.The goal of the GAMERA project was to rebuild the LFM code with two underlying design principles: retain and improve the high-heritage MHD numerics; and prepare the code for the exascale era.These improvements include 1) flexible grid specifications; 2) high-order grid metric calculations based on Gaussian quadrature; 3) seventh-order upwind spatial reconstruction; 4) non-clipping options in the limiters; 5) higher-order averaging-reconstruction method for axis singularity.While based on the best elements of the LFM legacy, GAMERA is designed to be easily applicable to two-or three-dimensional MHD simulations using general curvilinear computational grids adapted to specific problems.It is highly efficient; for example usable Earth or planetary magnetospheric simulations can be run on a laptop close to real time.
The goal of this paper is to provide a general description of the GAMERA code, including both the numerical schemes and the implementation details.The intent is that the paper will serve as a reference for future users and developers to use, modify and eventually extend the code to their own research applications.The description contains enough details so that the numerical techniques and equations used in the MHD kernel of the GAMERA code can be easily identified, or a functionally equivalent code can be built based on the details described in this paper.An example MATLAB code with the GAMERA algorithms described in this paper is provided in the supplement material for potential users to further understand the numerical algorithms.With a set of standard test problems as references, the quality of the numerical schemes used in GAMERA can be evaluated and improved when necessary.The paper is organized as follows: Section 2 describes the default MHD equations solved in the GAMERA code.Section 3 describes the detailed numerical schemes used in the GAMERA code, including grid discretization, time stepping, spatial reconstruction, flux functions, constrained transport algorithms in non-orthogonal geometry, and details about the code implementation including vectorization and parallelization.Section 4 shows results from five representative test simulation problems to demonstrate the quality of the numeric algorithms, and summary and conclusions are included in Section 5.

The Basic MHD Equations
The default equation set solved in the GAMERA code is the single-fluid, normalized ideal MHD in a semi-conservative form with the plasma energy equation: where ρ is the plasma mass density, u is the plasma bulk velocity, Ī is the unit tensor, P is the plasma thermal pressure, B is the magnetic field, and E = −u × B is the electric field based on the ideal Ohm's law.E P is the plasma energy defined as the sum of kinetic and thermal energy: with γ = 5 3 defined as the ratio of specific heat.The normalization of equations ( 1) -( 4) can be found in Appendix A. The fluid part of the semi-conservative form of the ideal MHD equation set can be written in a compact vector form: where U C is the vector form of the conserved density variables defined as: F U C is the vector form of the fluid part of the ideal MHD equations (ideal hydrodynamics): and M U C is the vector form of the magnetic terms from the Lorentz force defined as: The corresponding state vector of the primitive variables U P is defined as which is mainly used in the reconstruction module.Another important set of variables used in the MHD solver is the volume integrated conserved densities, with a vector form defined as: which are the actually conserved quantities when the MHD equations are evolved.
The use of the plasma energy equation has significant advantages in a number of problems, especially in modeling planetary magnetospheres with strong background magnetic fields and cold ambient plasmas (e.g., β < 10 −6 ).However in using the plasma energy equation instead of the total energy equation, the system of equations is no longer in a fully-conservative form.Non-conservative forms of the energy equation can lead to non-physical results if non-ideal processes occur, whether physical or numerical.However, Lyon et al. [11] showed that the use of plasma energy equation does not change the Rankine-Hugoniot relations to within the numerical truncation error.The result is independent of whether or not the electric field is carried by dissipative processes through the shock.The truncation level of numerical error is not quite as low as the round-off error from the total energy formulation, but it is sufficient to ensure that the behavior of shocks is correct to a good approximation.

Finite-Volume Spatial Discretization
In this section we describe the definitions of the curvilinear computational grids, the plasma and magnetic variables, and basic metric calculations such as face area, face normal unit vector and cell volume used in the non-orthogonal, finite-volume MHD solver.

Grid Definitions
The basic cell structure defined in the finite-volume solver is that of a general hexahedron without the requirement of coordinate orthogonality as shown in Figure 1.Although the physical cells used in the solver can be general hexahedra, the grid structure is logically Cartesian in the computational space.In other words, the grid used in the MHD solver is curvilinear.In the finite volume method, the differential form of the MHD equations ( 1) -( 3) are volume-integrated over individual cells.Thus the quantities then referred to are not the point values, as they would be in a finite difference formulation, but the volume-averaged values within a computational cell.Using Gauss's law, the conservative part of the volume-integrated form of Equation ( 6) becomes a surface integral: where U V is the vector form of the volume-integrated conserved variables as defined in (11).The Lorentz terms M U C in ( 6) are treated in a similar way using an operator splitting technique, which is discussed in Section 3.2.2.
The differential form of the Maxwell's equation ( 4) is integrated over cell interfaces such that the actual magnetic variable evolved in the code is the magnetic flux through cell interfaces, which is discussed in detail in Section 3.5.
While it is possible to solve for the non-Cartesian representations of the velocity and magnetic field using a finite volume technique, it generally leads to geometry-related source terms in the integral form of the MHD equations in which cell integrals must be done explicitly in (12).Therefore in the GAMERA code, all the cell-centered vector components are represented as Cartesian fields in order to avoid such source terms originating from general curvilinear coordinates.In the GAMERA code, the Cartesian components of the cell-centered velocity u (u x , u y , u z ) and magnetic field vectors B (B x , B y , B z ) are defined in the "Base Cartesian Coordinate System".
Figure 1: Left) The grid definition for cell center x C , cell face x F and cell corner x locations in a computational cell.
Right) The staggered grid definition for magnetic fluxes/fields and electric fields.
In the following sections, the Base Cartesian coordinates are denoted as (x, y, z), while the curvilinear coordinates are denoted as (µ, ν, ζ) as shown in Figure 1.In the computational space, the cell-centered index of the curvilinear grid is given by integers i, j and k for the µ-, νand ζ-direction, respectively.Thus the corresponding cell corner and face index is represented using half integers such as i ± 1 2 , j ± 1 2 , k ± 1 2 , respectively.Following these index notations, the spatial locations of the "cell centers" are denoted as x C i,j,k , the "cell corners" are then located at x i± 1 2 ,j± given that the ± 1 2 index indicates displacement halfway between cell centers.The spatial location of the face centers at the µ-, νand ζ-face are denoted as , respectively.Note that in the GAMERA code, the primary grid is the cell corner coordinates x i+ 1 2 ,j+ 1  2 ,k+ 1 2 , while cell centers, face centers and other grid metric calculations are derived from the cell corners.These grid metric calculations include the cell volume, face area and face-normal unit vector.Using the indices defined for the curvilinear grid, the cell volume is denoted as V i,j,k , the face areas at the µ-, νand ζ-interface are denoted as , respectively.Another set of grid metric quantity used in the MHD solver for the calculations of interface fluxes is the unit vectors of the interface normal directions, which is denoted as  grid definition and notations used in the following sections are listed in Table 1.The calculations of cell volume (V ), face-center locations (x F ) and face area (A) use the 12 th -order Gaussian quadrature and are described in Appendix B. The face-normal vector and the associated coordinate transforms are used in computing numerical fluxes through the cell interfaces, which is introduced in Section 3.4.
Based on the finite-volume discretization shown in Figure 1, for a single control volume indexed as (i, j, k), the volume-integrated form of Equation ( 12) becomes: where U V i,j,k is the volume-integrated conservative quantities and F is the numerical flux through the corresponding cell interfaces.

Variable Definitions
The MHD variables used in the GAMERA code are listed in Table 2.The detailed definitions and calculations about these variables are described in the following sections.

Variable Type
Variable Notation Location defined Property Coordinates

Plasma variables
The fluid state vector U V i,j,k of the mass, momentum and energy in the GAMERA code is defined as the volume integral of the conserved densities U C i,j,k over the controlled volume (i, j, k): Here, U C i,j,k are the mean densities (mass density, momentum density, energy density) within the cell (i, j, k).As a result, what is actually conserved in the MHD solver is the volume-integrated quantities U V i,j,k .In the reconstruction algorithm introduced in Section 3.3.1, the high-order reconstruction method for interface value used in the GAMERA code is based on the volume-integrated conserved variables U V i,j,k such that the variation of grid geometry is taken into account in the reconstruction process for non-uniform grids.The limiting step for the left-and right-state at an interface is applied to the primitive state vector U P i,j,k without the impact of variation from grid geometry.The details are described in Section 3.3.1.

Magnetic Fluxes and Fields
The MHD solver adapts the Yee-Grid [26] to non-orthogonal geometries in order to enforce the conservation of the volume-integrated magnetic divergence to round-off error, which is illustrated in the right panel of Figure 1.Therefore the primary electric and magnetic variables are defined on a staggered, non-orthogonal grid that are not co-located with the plasma variables at the cell centers.In the GAMERA code, the primary magnetic field variables evolved by the Faraday's law are the magnetic flux Φ through cell faces.Using the Φ µ component shown in the right panel of Figure 1 as an example, the magnetic flux Φ µ i+ 1 2 ,j,k threading a face in the µ-direction indexed as (i + 1 2 , j, k) is calculated as where B µ i+ 1 2 ,j,k is the mean magnetic field vector at the center of a cell face in the µ-direction.In the GAMERA code, since we track the magnetic flux Φ, the face-centered magnetic field vectors at face centers are nowhere computed directly.Using constrained transport based on the staggered grid, the magnetic fluxes Φ µ,ν,ζ are evolved using the cell-edge centered electric fields (or electric potential) surrounding the face flux through the integral form of the Faraday's law.For example using the Φ ζ i,j,k+ 1 2 component at the (i, j, k + 1 2 ) cell face as an example shown in the right panel of Figure 1, according to the Stokes' theorem, the integral form of the Faraday's law ( 4) is expressed by integrating the electric field along the green contour surrounding where E is the electric field component along the edge and dl is the corresponding length of the cell edge vector.The integral path on the RHS of ( 16) is defined as the closed contour a − b − c − d − a as illustrated in the right panel of Figure 1.In the MHD solver, the contour integral is done explicitly in a piece-wise way, i.e., the quantities E • dl are computed at each cell edges with the electric field E calculated using high-order reconstruction schemes adapted to cell edges.
Using the staggered discretization of the magnetic flux and electric field, it is straightforward to show that the divergence of the magnetic field ∇ • B is conserved in an integral sense [27], as long as ∇ • B = 0 is satisfied initially.In other words the action of an electric field along any cell-edge leaves ∇ • B unchanged during computations.This staggered discretization of the magnetic field allows non-linear limiters which modify the electric field locally to still conserve the solenoidal behavior of the magnetic field, which is essential to the quality of the numerical solutions to MHD equations.
In order to compute the Lorentz force terms −∇ • I B 2 2 − BB on the bulk plasma using magnetic flux functions described in Section 3.4, the Cartesian magnetic field vector denoted as B x,y,z i,j,k (defined in the Base Cartesian coordinates) co-located at cell centers x C i,j,k with plasma variables are needed.To compute cell-centered magnetic field vectors B x,y,z i,j,k using face-centered magnetic fluxes Φ µ,ν,ζ , consider the volume integral of ∇ • (rB) within a controlled volume V given that ∇ • B ≡ 0: where B is the mean magnetic field vector within a finite-volume cell, and V is the corresponding cell volume.Applying Gauss's law to the LHS of equation (17), the volume integral is then converted to the following face integral: Using equations ( 17) and (18) with ∇ • B ≡ 0, the cell centered magnetic fields B is computed as The RHS of equation ( 19) is a face integral using the primary magnetic flux variables Φ µ,ν,ζ .In the GAMERA code, a second order approximation based on the face-centered magnetic flux is used to approximate the face integral on the RHS of equation (19), i.e., the cell-centered magnetic fields B xyz i,j,k using estimations of magnetic flux at cell centers are computed as The above formulation is simple and works for any hexahedral cell defined by its corners.Note that ∇ • B = 0 is required in the above equation calculating the Cartesian components of the magnetic fields located at cell centers.

Time Stepping for Temporal Evolution
This section provides a general description of the numerical algorithms used in the GAMERA code for evolving the MHD equations ( 1)-( 4), including the time-stepping methods, the Courant condition calculations, the operator splitting technique for handling the magnetic terms M(U) in ( 6) and the framework for handling system of equations Figure 2: The main GAMERA MHD algorithm, including the predictor step, the corrector step, the fluid solver and the Maxwell solver.in non-orthogonal, curvilinear geometries.Figure 2 illustrates the main algorithms of the MHD solver including the predictor/corrector steps, the fluid solver for evolving the plasma variables in Equations ( 1)-(3), and the Maxwell solver updating the magnetic flux/fields in Equation (4).
For time-stepping we use the Adams-Bashforth second order scheme as the default algorithm.The predictor step is a linear extrapolation to the half time level (t n+ 1  2 ) based on the state variables of the current (t n ) and previous (t n−1 ) time steps.For the fluid part of the MHD equations, the corrector step is a full flux calculation for the integral form of Equation ( 1)-(3) based on the predicted state using the Partial Interface Method (PIM) developed by Lyon et al. [11], which is denoted as the "Fluid Solver" and introduced in Section 3.2.2.For evolving the integral form of the Faraday's law, the corrector step is the calculation of electric fields along cell edges using a constrained transport method based on the predicted velocity and magnetic flux, which is denoted as the "Maxwell Solver" and described in Section 3.5.Extensions to multi-stage time stepping methods such as Runge-Kutta are possible based on the existing algorithms of the GAMERA code with appropriate modifications.

The Predictor Step
Assuming Q n−1 is a state variable at time step n − 1, Q n is the corresponding state variable at time step n.The predictor step for Q n+ 1 2 at n + 1 2 is computed using the following linear extrapolation in time: where ∆t n−1 and ∆t n are the time steps at t = n − 1 and t = n, respectively.In the GAMERA code, the predictor steps are applied to the primitive state variables U P (for fluid flux/stress calculations), the cell centered Cartesian magnetic field components B x,y,z (for magnetic stress calculations) and the face magnetic fluxes Φ µ,ν,ζ (for electric field calculations), respectively.Once the predictor step is done, the following corrector step is applied to the discrete form of Equation ( 13) for evolving the fluid part of the MHD equations: where F U P,n+ 1 2 is a fluid flux calculation using the primitive variables U P,n+ 1 2 at half time step, is the magnetic stress calculation using both U P,n+ 1 2 and B n+ 1 2 .∆U F is the volume-integrated changes in mass, momentum and energy from the fluid flux/forces, and ∆U B is the corresponding changes originating from the magnetic forces.The detailed algorithms for computing the ∆U F and ∆U B terms are discussed in the Section 3.4.For the integral form of the Faraday's law (16), the corrector step follows where E n+ 1 2 is the edge electric field calculated using the predicted variables Φ n+ 1 2 and U P,n+ 1 2 .The details of the electric field and magnetic flux evolution computations are described in Section 3.5.
This second-order Adams-Bashforth time stepping scheme is simple and robust with the presence of non-linear limiters with the Total variation diminishing (TVD) property [11].The main advantage of using the Adams-Bashforth scheme is the small amount of memory and small number of computations needed for the predictor step.Similarly, a third-order Adams-Bashforth scheme for the predictor is written as which involves two previous time steps of the state vector Since the default Adams-Bashforth scheme is explicit in time, the time step ∆t is determined by the Courant-Friedrichs-Levy (CFL) condition.The explicit time step is calculated as: where N CF L is the Courant Number between 0 and 1, ∆x is an effective minimum "length" of a hexahedron cell as shown in Figure 1.In a computational cell indexed as (i, j, k), the corresponding minimum cell length ∆x is approximated as the cell volume V i,j,k divided by the maximum face area: and v is the maximum wave speed within the computational cell.For the ideal MHD equations, the maximum wave speed is calculated as where |u| is the magnitude of the plasma bulk velocity, V 2 A + C 2 S is the magnetosonic speed calculated using the Alfvén speed V A and the plasma sound speed C S within cell (i, j, k).The Courant number N CF L used in the default GAMERA solver is a function of numerical diffusion introduced in the default non-linear limiter [14], which is described in Appendix C.

The Corrector Step
In the corrector step, we use the "Partial Interface Method" (PIM) described by Lyon et al. [11] in the LFM global magnetospheric model for handling the system of fluid equations in multi-dimensional, non-orthogonal geometries.
The PIM provides a general framework that has the advantage of being readily adapted to arbitrary high-order spatial reconstruction methods as well as underlying numerical flux functions.The default reconstruction algorithms and flux functions used in the GAMERA code are described in the next section.
The PIM algorithm loops over each curvilinear direction and is split into the following steps: Step 1. Start from the predictor U The algorithm of computing ∆U F and ∆U B based on the predictor state is summarized in Figure 3, where steps 2 and 3 are basically in the same vein as a typical one-dimensional TVD scheme for interface states.Note that in the PIM, the reconstruction algorithm is essentially one-dimensional; that is, the provisional values at the cell interfaces are reconstructed along each curvilinear grid coordinate µ, ν, η, and no multi-dimensional corrections are applied to modify the provisional value U interface i + 1  2 are computed in a one-dimensional stencil along the curvilinear direction.In step 4, the default flux functions used in the PIM for both the fluid and magnetic terms are centered flux functions with the following form: developed by Croisille et al. [28].Note that if 6), which is the exact form of the fluid part of the ideal MHD equations.It is only when U L = U R that any numerical dissipation is introduced, although the amount of numerical dissipation introduced by the Gas-Kinetic scheme F GK (•) is implicit.The numerical diffusion of that scheme is not as obvious as a first-order Rusanov scheme and has no simple description as does Rusanov.In the case of PIM, many different numerical flux functions can also be used to replace the default gas-kinetic scheme.For example a global heliosphere implementation of the LFM code used the HLL function function for solving the MHD equations [22].Note that the default GAMERA MHD solver uses a 7 th -order reconstruction method to calculate . With such high-order reconstruction schemes, standard MHD test simulations suggest that the numerical solutions are not sensitive to the choice of the low-order flux function used in the PIM.
In Step 5, the surface integral evaluations for the mass and energy flux are straightforward since the corresponding directions of mass flux and energy flux are normal to cell interfaces.The face integrals for momentum flux calculations in ∆U F require coordinate transform from the local face-normal coordinate system (defined in Section 3.4.3) to the Base Cartesian System in order to evolve the Cartesian components of the plasma momentum vector.This coordinate transformation enables the fluid solver to track the Cartesian components of the plasma momentum without the requirement of grid orthogonality.
When ∆U F and ∆U B are calculated using the PIM algorithm, an operator splitting technique is used to evolve the fluid part of the MHD equations ( 1)-( 3).This operator splitting technique enables solving for the plasma pressure P without implementing the u • ∇ • I B 2 2 − BB term explicitly in the plasma energy equation.The algorithm shown in Figure 2 illustrates the process of updating the plasma state vector using the operator splitting technique described above.The first step of the operator splitting is solving the ideal gas dynamics equations using ∆U F only: and get ρ n+1 and P n+1 with an intermediate velocity update u * , then apply the Lorentz force terms ∆U B to the intermediate velocity u * for the final update of velocity u n+1 : The operator splitting technique is implemented by separating the fluid force terms ∇ • (ρuu + IP ) and the magnetic force terms ∇ • I B 2 2 − BB in the momentum equations, thus the plasma pressure P can be solved directly from the ideal gas dynamics part of the MHD equations ( 1)-( 3) before applying the Lorentz force ∇ • I B 2 2 − BB in the momentum equations, and the u • ∇ • I B 2 2 − BB term is no longer needed in solving the thermal pressure P from the plasma energy equation at each time step.The operator splitting technique simplifies the implementation of the semi-relativistic (Boris) correction [29] and the multi-fluid extension [15] for planetary magnetosphere simulations significantly.

Calculation of Interface States
In this section, we describe the default algorithms for computing the left-state variables (U L ) and right-state variables (U R ) at cell interfaces for flux evaluations in the PIM framework shown in Figure 3.These interface values are computed through two consecutive steps: 1) a high-order reconstruction step and 2) a limiting step.In the reconstruction step, the high-order reconstruction is performed on the volume-integrated state variables U V (mass, momentum, plasma energy) in order to incorporate the variations in grid geometry.After the reconstruction step, the reconstructed U V at cell interfaces are then converted to the conservative densities U C with the geometry variations removed.In the limiting step, the reconstructed U C at cell interface is first converted to the primitive state U P (including density, velocity, pressure and magnetic fields).These interface primitive variables at cell interfaces are then split into left-and right-state using the non-linear limiter algorithm.In the GAMERA code, the reconstruction and the limiting modules are implemented separately in order to make the choice of numerical schemes flexible.In the following sections we introduce the high-order reconstruction schemes and the Partial Donor Cell Method (PDM) limiter implemented as the default in the MHD solver.

High-Order Reconstruction
The interface values used in the GAMERA code for flux calculations are high-order approximations of the primitive variables U P and magnetic fields B xyz at cell interfaces.To incorporate spatial variations originating from curvilinear geometries, the high-order reconstruction scheme operates on the volume-integrated fluid variables U V , V • B xyz defined at cell centers and Φ µνζ at cell faces.In the PIM framework, the reconstruction process for computing interface values is one-dimensional along each curvilinear direction, without multi-dimensional corrections which simplifies the calculation in non-orthogonal geometries.For example, in order to estimate high-order approximations for interface states U V i+ 1 2 along the µ-direction, first define a one-dimensional integral function G(x µ ) along the µ-direction as where U V is the vector form of the volume-integrated conservative variables, x µ is the curvilinear spatial coordinate along the µ-direction in the computational space as shown in Figure 4. Therefore the conserved quantity Figure 4: A one-dimensional, high-order reconstruction stencil for computing The order of the reconstruction method to get interface values depends on how accurately the first derivative term in Equation ( 33) is approximated.An example of an eight-cell reconstruction stencil along the µ-direction is shown in Figure 4.If applying a second order, centered approximation for the numerical derivative ∂ ∂xµ G (x µ ) at the interface i + 1 2 as an example, the RHS of Equation ( 33) is calculated as: Then the conserved state variables U C i+ 1 2 at interface i + 1/2 are calculated as where V i+ 1 2 is an estimation of "cell volume" at the corresponding cell interface i + 1 2 assuming the cell volume V i varies smoothly in the µ-direction.A 2 nd -order centered approximation of interface volume for the 2 nd -order interface reconstruction in (35) is calculated as is then computed for the limiting step described in the next section.The interface magnetic field B i+ 1 2 is calculated using the same algorithm by replacing U V with V i,j,k • B i,j,k .In the original LFM code, the reconstruction module uses the eight-cell stencil shown in Figure 4 with the following 8 th -order centered reconstruction for computing high-order estimations of interface values: where V 8th is an 8 th -order approximation of the interface volume calculated as: This 8 th -order centered reconstruction scheme is chosen as the default method for computing interface state variables in the LFM code based on its high-resolving power of contact discontinuities according to Lyon et al. [11].
In fact, the reconstruction method for estimating a high-order interface value is not necessarily centered.Highorder upwind reconstruction schemes are also implemented in the GAMERA code.Using the second-order centered reconstruction as an example, the interface variable is reconstructed as Combine the above second order approximation for U V i+ 1 2 with a second-order centered approximation for the first derivative of U V evaluated at i + 1 2 , the centered reconstruction becomes the first-order upwind method [30]: Equation ( 40) suggests that upwind reconstruction with order of n − 1 can be derived from a n th -order (n even) centered stencil by canceling the outer most cell in the downwind direction using a numerical first derivative approximated to n th -order accuracy at the cell interface.For example, combining the 8 th order centered interpolation for U V i+ 1 2 with an ∂xµ , the following 7 th -order upwind reconstruction method is obtained: In the above 7 th -order reconstruction, the right-most cell i + 4 shown in Figure 4 is removed from the reconstruction stencil, which makes the reconstruction shifted towards left (the upwind direction).Reconstruction coefficients for left interface states from 1 st -order to 12 th -order are listed in Table 3.The right interface states are calculated using the same reconstruction coefficients by switching the upwind direction of a reconstruction stencil.For the centered reconstruction methods (e.g., 8 th -order), the left-and right-state are computed using the same eight-cell stencil and coefficients due to symmetry; for the upwind reconstruction methods, (e.g., 7 th -order), the left-and right-state are computed using two different seven-cell stencils shifted towards the left and right side of the interface, respectively.In the GAMERA code, the solver uses the 7 th order upwind reconstruction scheme as the default choice, while the original LFM MHD kernel uses the 8 th order centered reconstruction.The 8 th -order reconstruction scheme has a leading truncation error of a dispersive 9 th -order derivative, thus a sudden change in gradient (involving short wavelengths) may excite spurious unphysical oscillations due to the 9 th -derivative dispersion term in the truncation error.However in the 7 th -order reconstruction scheme, the leading truncation-error term is an 8 th -order spatial derivative, which is dissipative instead of dispersive.This is an important improvement in the original LFM algorithm.Numerical experience over the past decade or more has shown that for numerical solutions of convection-dominated fluid flows, the leading truncation  3: Centered and Upwind reconstruction coefficients for U i+ 1  2 , up to 12 th -order.The upwind coefficients are for the left-state interface values with the stencil shifted towards the left.
error in a convection algorithm should be dissipative (an even spatial derivative term) rather than dispersive (an odd spatial derivative term), which should also be of higher order than modeled physical diffusive terms (if any) [33].Thus the odd-order (order ≥ 3) reconstruction schemes are more natural choices satisfying the criteria of 1) dissipative truncation error terms and 2) higher-order than physical diffusion terms.Thus the GAMERA code uses a 7 th -order reconstruction method as the default choice for upwind spatial reconstruction, and the 8 th -order reconstruction method is chosen as the default centered method for reference.The reason for choosing such high-order reconstruction methods as the default in the MHD solver is explained in Appendix D, which is based on a simple measure of optimizing between low numerical diffusion and high computing efficiency as well as the need to resolve both physical and grid structure.

Partial Donor Cell Method
When discontinuities occur within a reconstruction stencil, a problem arises in using the high-order reconstruction schemes since spurious undershoots or overshoots may occur near discontinues.To avoid this problem, nonlinear switchers are used to "correct" the high-order reconstructed values on both sides of the interface when necessary.In the limiting module the Partial Donor Cell Method (PDM) limiter developed by Hain [14] is implemented as the default choice, which does not depend on the numerical order of interface reconstruction.The basic idea of the PDM limiter is that the algorithm monitors sharp discontinuities; if a sharp discontinuity is identified, a "limited" value is used in order to the keep the solution from overshoots/undershoots; otherwise a high-order approximation is always chosen for the interface state to provide a more accurate estimation.The details about the PDM limiter can be found in Hain [14].
Here we provide a simplified description of how the PDM limiter works adapted from Huba [31].
Consider a density structure being advected at a constant velocity V towards the right direction as shown in Figure 5, the PDM limiter determines whether the high-order estimation of the left-state value ρ HO L i+ 1 2 at interface i + 1/2 needs to be "limited" for monotonicity preserving purpose.In Figure 5, ρ i−1 is the density in cell i − 1; ρ i is the density in cell i, ρ HO i+1/2 is the high-order reconstructed density at cell interface i + 1/2, and ρ P DM i+1/2 is the limited value (or the PDM value) of the interface density on the left side which will now be determined.
Assume the density structure is advanced for one time step ∆t, it advects towards the right by a finite distance of V ∆t.After one advection step, the amount of mass entering cell i from cell i − 1 is ρ i−1 V ∆t, which is illustrated by the orange shaded area in Figure 5; the amount of mass leaving cell i and entering cell i + 1 is ρ P DM i+1/2 V ∆t assuming that the left state is the PDM value that guarantees no overshooting/undershooting, which is illustrated by the green shaded area in Figure 5. Therefore the total density change in cell i after one time step is calculated as (ρ i−1 − ρ P DM i+1/2 )V ∆t.Since the maximum density increase allowed in cell i without spurious overshoot is (ρ i−1 − ρ i )∆x, which is denoted by the purple area, the left state of the PDM value at cell interface i + 1/2 is derived by balancing the two quantities: which gives the PDM value: where = V ∆t/∆x which is in the range between 0 and 1.By balancing the density entering and leaving cell i, it is clear that the ρ P DM i+1/2 is the minimum value of the left state density allowed to leave the cell i without spurious overshoot.If the left state density is lower than this value, spurious density overshoot occurs within cell i.By choosing the ρ P DM i+1/2 as the left state, the density on the left-side of the interface is "limited".Note that if = 1, then ρ P DM i+1/2 = ρ i , which is the first-order Donor Cell method with excessive amount of numerical diffusion for non-linear problems.By choosing < 1, the above scheme takes a portion of the Donor Cell solution as the limited value and reduces the amount of numerical diffusion significantly.
To achieve a high-order approximation for the left state, it is not necessary to use the PDM value all the time.As shown in Figure 5, there are three density values to choose as the left state at interface i + 1/2 : ρ i (first-order), ρ HO i+1/2 (7 th -order reconstruction as the default) and ρ P DM i+1/2 (the limited value).The rationale for choosing the final left state value is as follows.In general, one would want to use the high-order reconstructed value because it provides the best accuracy for estimating the left state.Since ρ P DM i+1/2 is the minimum amount of density allowed for a left state leaving cell i, as long as ρ P DM i+1/2 < ρ HO i+1/2 < ρ i , the high-order value ρ HO i+1/2 is the correct choice for the left state.However, when ρ HO i+1/2 < ρ P DM i+1/2 < ρ i , the PDM value ρ P DM i+1/2 must be chosen as the left value in order to avoid overshoots.In other words, the left state is always chosen as the median value of the three interface values: ρ i , ρ HO i+1/2 , and ρ P DM i+1/2 .The PDM limiter is based on balancing the amount of flux moving through a controlled volume rather than on constraining the slope of the underlying reconstruction polynomial for computing interface values, which is similar to the "Universal Limiter" developed by Leonard and Niknafs [33].Thus the interface value ρ HO i+1/2 can be replaced by arbitrary high-order reconstructed values.The right-state density is derived in the same way through the advection of a density structure towards the left.The performance of the PDM limiter combined with different orders of reconstruction methods is demonstrated in Appendix D based on quantitative comparisons of 1-D linear advection test results.
For the linear advection equation, it is straightforward to show that the PDM limiter is TVD [14], thus it tends to "clip" local extrema as typical TVD limiters do.An optional non-clipping algorithm is implemented in the limiting module to preserve local extrema for simulations that are sensitive to the numerical errors introduced due to clipping of local extrema, e.g., the Harris current sheet equilibrium simulation.A simple non-clipping algorithm developed by Leonard and Niknafs [33] is implemented in the limiting process in order to distinguish between artificial numerical peaks and true physical extrema.Figure 6 shows the difference between a physical extremum and a numerical extremum within one reconstruction stencil centered at cell i.If a local extremum is associated with short-wavelength oscillations, the limited value is chosen for the interface state; however, if the curvature of a peak follows the shape of a local extremum, the unlimited high-order value is chosen to avoid clipping.Using the stencil for calculating the left-state U i+ 1  2 at the interface i + 1 2 shown in Figure 6 as an example.The calculation of the local extrema indicator (LEI) is illustrated in Figure 6, by first computing the differences between each pair of consecutive cell center values: When the non-clipping option is switched on, after the limiting step, the final interface state is selected based on the value of the LEI, if LEI = T RU E, the unlimited, high-order approximation (e.g., the ρ HO i+1/2 ) is used for the interface state regardless of the ρ P DM i+1/2 constraint; otherwise if LEI = F ALSE, the limited value is chosen to avoid possible spurious overshooting/undershooting.
Combining the reconstruction and the limiting step, the algorithm to compute the left and right-state variables along one-dimensional stencils (e.g., in the µ-direction) is summarized as follows: Step 1. Compute the volume-integrated conserved quantities in each controlled volume using Step 2. Compute a high-order (default 7 th -order) reconstructed volume-integrated quantity U V i+ 1 2 ,j,k ; Step 3. Estimate a high-order (default 8 th -order) volume V i+ 1 2 ,j,k at the interface location where Step 5. Split the interface primitive variable U P i+ 1 2 ,j,k into the left-and right-state variables using the PDM limiter.
Step 6 (optional).Apply the non-clipping sweep to decide whether the PDM limiter is needed or not.

Gas-Kinetic Flux Functions
The default flux functions for the fluid part of the MHD equations are based on one-dimensional gas-kinetic schemes for the Euler equations in the face-normal coordinate system (detailed in Section 3.4.3).In a local, orthogonal coordinate system (x n , x 1 , x 2 ) with x n the direction normal to a cell interface and x 1 , x 2 the two corresponding orthogonal directions tangential to the cell interface, the one-dimensional mass, momentum and plasma energy equations ( 1)-(3) in this (x n , x 1 , x 2 ) coordinate system are written as: where E P is the plasma energy defined in (5), and B 2 /2 is the magnetic energy calculated as The first part on the RHS of Equation ( 46) corresponds to the ideal gas dynamics with velocity components in the (x n , x 1 , x 2 ) coordinate system, which gives the ∆U F terms in Equation (22).The second part on the RHS of Equation ( 46) is the Lorentz force ∇ • I B 2 2 − BB in the momentum equation ( 2), which gives the ∆U B terms in Equation ( 22).The third part on the RHS of Equation ( 46) is the work done by the Lorentz force u • ∇ • I B 2 2 − BB in the plasma energy equation (3), which is not computed explicitly for solving the plasma pressure P due to the operator splitting technique as discussed in Section 3.2.2.The gas-hydro flux (fluid stress) and the magneto-hydro flux (magnetic stress) for the above one-dimensional MHD equations are calculated separately, which simplifies the implementation of the Semi-relativistic (Boris) correction for magnetospheric simulations.The default solver uses a Maxwellian-based gas-kinetic flux function to compute the fluid fluxes terms ∆U F and uses another Maxwellian-based magneto-gas kinetic flux function to compute the magnetic stresses terms ∆U B .

Fluid Fluxes
The default flux functions used in the solver for the fluid flux calculations are gas-kinetic schemes based on Maxwellian distributions, which considers the fluid on each side of a cell interface as a microscopic distributions of non-interacting particles [28].These particles are allowed to free-stream across the interface for one time step, then the macroscopic conserved variables including mass, momentum, and energy are derived from the zeroth, first and second moment integral of the corresponding distribution functions on both sides of the interface.These updated macroscopic variables are then used to construct new distribution functions for the next time step.
Using the left primitive state variables U P,L = ρ L , u L , P L T calculated in the reconstruction module, the distribution function f L (v) for the plasma population on the left side of the cell interface is a Maxwell-Boltzmann distribution function: Similarly, the corresponding distribution function f R (v) for the plasma population on the right side of the cell interface is where ρ L and ρ R are the left and right state of plasma density, P L and P R are the left and right state of plasma pressure, u L xn and u R xn are the left and right plasma bulk velocity in the x n -direction normal to cell interfaces, respectively.The calculations of the u L xn and u R xn from the interface velocities u L and u R is described in the next section.For one-dimensional MHD flows described by Equation (46) with particles free-streaming in the x n -direction, the moment integrals in this direction determine the form of flux functions.Thus in the local (x n , x 1 , x 2 ) coordinate system the mass flux, the fluid part of the momentum flux and the energy flux across the cell interface i + 1  2 are calculated by integrating the corresponding distribution functions for the left-going particles on the right-side of the interface and the right-going particles on the left-side of the interface: where F n FLUID is the vector form of the fluid fluxes at interface n (n = µ, ν, ζ); F ρ is the mass flux, F ρu is the momentum flux vector and F E P is the plasma energy flux calculated in the (x n , x 1 , x 2 ) coordinates.u L x1,2 and u R x1,2 are the corresponding left-and right-interface tangential velocity components transformed into the local (x n , x 1 , x 2 ) coordinate system, respectively.After the fluid fluxes are calculated, the ∆U F term for the correction step ( 22) is computed as where A s is the face area and n s is the corresponding face-normal unit vector.If the left and right states are identical, i.e., and the RHS of one-dimensional ideal gas dynamics equations are recovered through the flux functions (49).The detailed derivations for evaluating of the moment integrals in Equation ( 49) can be found in [34].Here we only give the mathematical expressions of the Maxwellian distribution based gas-kinetic flux scheme for computing fluid fluxes: The zeroth and the first velocity moments of the left and right distribution functions used in Equations ( 51) are calculated as: , the following relationships hold: The mass flux F ρ and energy flux F E P in equations (51) are in the x n -direction normal to cell interfaces: Thus the mass and energy flux are used to compute the surface integrals in Equation (50) directly: However, the momentum flux vector F ρu calculated using equation ( 51) is defined in the local (x n , x 1 , x 2 ) coordinate system: which cannot be used to compute the surface integrals in Equation (50) directly in order to update the base Cartesian components of the momentum (ρu x , ρu y , ρu z ).Thus after evaluating the fluid fluxes in the face-normal coordinate system (x n , x 1 , x 2 ) using the one-dimensional flux function form (51), the momentum flux vector F ρu at cell interface is then rotated back to the base Cartesian coordinate system for updating the Cartesian component of the momentum (ρu x , ρu y , ρu z ).Then the changes to the plasma momentum ∆U F ρu in the base Cartesian system is computed as: where T −1 is a matrix for transforming the momentum flux vector from the (x n , x 1 , x 2 ) system to the base Cartesian system (x, y, z).The detailed calculation of the transform matrix T −1 is described in Section 3.4.3.This rotation method for evaluating momentum fluxes in the base Cartesian coordinate system (x, y, z) does not require the orthogonality of the physical grid and works with any microscopic distribution functions in the gas-kinetic flux scheme.
The gas-hydro distribution function for the fluid flux calculations does not depend on the local Alfvén speed, therefore the numerical diffusion speed near discontinuities, although non-explicit, is only related to the bulk flow and plasma thermal speed.In very low-beta plasma simulations (e.g., β < 10 −4 ) such that magnetosonic waves dominate the dynamics, the use of thermal speed in the gas-hydro distribution functions above is not adequate for describing the behavior of plasma associated with Alfvén waves.Thus additional numerical diffusion is needed to damp spurious oscillations at the Alfvén speed when the left and right states are not equal.For gas-kinetic fluid flux functions, the numerical diffusion speed for the fluid flux is adjusted based on the average Alfvén speed at the interface, which is simply implemented via incorporating an additional diffusion term that follows a Rusanov type of numerical flux calculation centered at cell interfaces: where V A i+ 1 2 is chosen to be the average Alfvén speed V A at the cell interface i + 1 2 : Note that the change of the momentum flux vector ∆U F ρu in Equation (62) has already been transformed back to the base (x, y, z) system.Thus the diffusion terms for the momentum fluxes are computed using the Cartesian components of the interface bulk velocities (u L and u R ) directly.The additional diffusion terms in Equation (62) only apply to the interfaces where the left-state and right-state variables are not equal.In smooth flow regions, the left-and right-state variables are equal or the differences are negligibly small, no numerical diffusion is introduced to the interface fluxes through the above flux function.

Magnetic Stresses
The calculations of magnetic stress terms go through a similar process as computing the fluid fluxes, using similar Maxwellian-based, magneto-gas kinetic distribution functions for the magnetic fields on the left-and right-side of the interfaces [34]: where /2 and Φ 0s = B 0 • n s are the magnetic energy and the magnetic flux of the background field at interface s, respectively.The background magnetic field components B 0 on cell interfaces used in Equation (71) are computed using a 12 th -order 2-D Gaussian quadrature rather than point evaluations to improve the implementation of force-free background magnetic fields, i.e., ∇ × B 0 = 0.For example, at µ-faces: The MHD solver in the GAMERA code only tracks Cartesian components of the plasma velocity (u x , u y , u z ) regardless of the grid geometry.However, as discussed in the previous section, the calculations of fluid flux are done in the local face-normal coordinate system (x n , x 1 , x 2 ), which is usually not the same as the base Cartesian system.Thus after the interface reconstruction of the Cartesian velocity components (u x , u y , u z ) and magnetic field (B x , B y , B z ) in the µ-, νand ζ-directions, corresponding coordinate transforms at cell interfaces are needed to rotate the left and right velocity and magnetic field from the base Cartesian coordinate system (x, y, z) into the face-normal coordinate system (x n , x 1 , x 2 ), in order to evaluate the numerical fluxes using the Gas-Kinetic flux functions described in Section 3.4.For updating the plasma momentum after computing the momentum flux vector in local face-normal coordinate systems at each interface, inverse transforms are needed to rotate the face-integrated momentum flux from the (x n , x 1 , x 2 ) system to the base Cartesian system (x, y, z).Then the Cartesian components of the plasma momentum (ρu x , ρu y , ρu z ) are updated as described in the previous section.This transform -inverse transform algorithm enables the MHD solver to evolve the plasma velocity components in the base Cartesian coordinate system, which is independent of the curvilinear coordinate system used to define the computational grid.Therefore orthogonality of the curvilinear grid is not required in the MHD solver, and the numerical grid can be adapted to simulate problem specific flow patterns such as planetary magnetospheres.
Using the orthogonal, µ-face normal coordinate system (n µ ,n 1 ,n 2 ) at cell interfaces indexed by (i + 1 2 , j, k) illustrated in Figure 7 as an example, where n µ is the unit vector normal to the µ-face, n 1 and n 2 are the two unit vectors perpendicular to n µ .The n 1 ,n 2 and n µ vectors from an orthogonal system, which are calculated using the four corner grid points forming the µ-face shown in Figure 7: Figure 7: The local µ-face normal coordinate system (n µ ,n 1 ,n 2 ) in a finite volume cell indexed as (i, j, k).
the transform of the interface velocity vectors from the base Cartesian (x, y, z) coordinate system to the face-normal (n 1 , n 2 , n µ ) coordinate system is done as follows: where u n1,n2,nµ = u n1 , u n2 , u nµ is the velocity vector in the face-normal coordinate system (n 1 , n 2 , n µ ) at µ-faces, and u x,y,z = (u x , u y , u z ) is the corresponding velocity vector in the base Cartesian coordinate system (x, y, z).Since the transform matrix T in the above equation is an orthogonal matrix, the inverse transform to rotate the stress F FLUID from the (n 1 , n 2 , n µ ) system to the base (x, y, z) system is simply done using the transpose of T: where F n1,n2,nµ ρu is the momentum flux vector in the (n 1 , n 2 , n µ ) coordinate system described in Equation (60), and F x,y,z ρu is the corresponding vector of momentum flux in the base Cartesian coordinate system (x, y, z).With the coordinate transform, the surface integral of the fluid stress for ∆U F ρu is simply calculated using (61).However the surface integrals of mass flux and energy flux are in the n µ direction which are used to compute the surface integrals directly without coordinate transforms as shown in ( 58) and (59).For the magnetic stresses, since only the zeroth moment of the distribution function is involved, the stress calculations on cell interfaces are done directly using the x-, yand z-components of the interface magnetic fields and the n 1 , n 2 , n µ vectors in the base Cartesian system as shown in Equation (68), without the above transform/reverse transform procedures.

Constrained Transport for Non-Orthogonal Grids
In this section we describe the numerical schemes for calculating the electric fields defined at cell edges for evolving the magnetic flux Φ through Faraday's law.The calculation of the electric field is handled through the implementation of a constrained transport method based on the same high-order reconstruction schemes (described in Section 3.3.1)adapted to non-orthogonal, staggered grids.As a main part of the corrector step illustrated in Figure 2, the predicted velocity u i+ 1 2 and magnetic flux Φ i+ 1 2 are used in the Maxwell solver for computing the electric fields.Figure 8 shows the algorithms of computing the electric field defined at cell edges using the predicted bulk velocity u n+ 1 2 and magnetic flux Φ n+ 1 2 , with details described in the following sections.

Calculation of Electric Fields
Figure 9 shows a schematic of the grid geometry for computing the electric field component E ζ at ζ-edges.The E µ and E ν components along µand ν-edges are calculated using the same algorithm as E ζ , with corresponding rotations of the computational grid.Using the example shown in Figure 9, the electric field E ζ is located at the cell edge indexed as i + 1 2 , j + 1 2 , k , which is surrounded by four neighboring control volumes whose cell centers are indexed as (i, j, k), (i + 1, j, k), (i, j + 1, k) and (i + 1, j + 1, k), respectively.To get a high-order estimation of the electric field component E ζ at the ζ-edge, both face-centered magnetic flux and cell-centered plasma velocity are reconstructed at cell edges using the high-order reconstruction method described in Section 3.3.1.For the velocity u i + 1 2 , j + 1 2 , k at cell edges, this is done by first doing a high order reconstruction to µ-interfaces and then reconstructing these µ-interface values in the ν-direction towards the edge (or corner as it appears in the top view of the 2D slice shown in the right panel of Figure 9): Figure 9: A schematic showing locations of various quantities needed for calculation of the electric field at cell edges.
The default choice for reconstructing edge velocity is an 8 th -order centered scheme described in Equation (37).The estimation of magnetic field vectors at the cell edge is more complicated than computing the velocity vectors, which requires two reconstruction steps for both magnetic flux and face-normal vectors.The first step is to reconstruct the magnetic flux at cell edges.Since magnetic flux Φ µ and Φ ν are already defined on cell interfaces, it only requires one single reconstruction to the corner along the νand µ-direction, respectively.Using the Φ µ component of magnetic flux along the ν-direction as an example shown in Figure 10, after applying the one-dimensional reconstruction algorithm to the magnetic flux Φ µ i+ 1 2 ,j,k in the ν-direction, two "interface" states of magnetic flux, Φ are computed at the cell edge i + 1 2 , j + 1 2 , k .To estimate the magnetic field strength at cell edges, an eighth-order interpolation is applied to the face area A µ i+ 1 2 ,j,k along the ν-direction to find an estimation for the area A µ i+ 1  2 ,j+ 1 2 ,k at the cell edge for the estimated interface fluxes Φ (µ,L) 2 ,k , as shown by the shaded area in Figure 10.Then the magnetic field strength in the µ-direction at the cell edge i + 1 2 , j + 1 2 , k is calculated as the average of the left and right state of the magnetic flux divided by the estimated interface area: Figure 10: The one-dimensional reconstruction stencil for Φ µ and n µ in the ν-direction for estimating B µ avg .
The next step is to estimate the direction of B µ avg at the cell edge i + 1 2 , j + 1 2 , k , which is denoted as μinterp in Figure 10.This is done by interpolating the face-normal unit vectors nµ to cell edge using the same 8 th -order reconstruction algorithm.The expression for the nµ vectors are defined in equation ( 75), and the high-order reconstructed nµ vector at the cell edge is denoted as μinterp 2 ,k .Similarly, the average magnetic field in the ν-direction is calculated by applying the same reconstruction algorithm on the magnetic flux component Φ ν i,j+ 1   2 ,k and the corresponding ν-face area A ν i,j+ 1 2 ,k the µ-direction: The direction of B ν avg is also interpolated from the face-normal nν vectors using the same high-order interpolation schemes.The high-order interpolated nν vector at the cell edge is denoted as νinterp i+ 1  2 ,j+ 1 2 ,k .Using the reconstructed edge velocity and average magnetic field vector, the calculation of the ζ-component of the electric field is basically a Rusanov scheme adapted to a cell corner.With the edge velocity vector u i+ 1 2 ,j+ 1 2 ,k and mean magnetic fields B µ,ν avg , the electric field component E ζ at the ζ-edge is calculated as: where j ζ is the component of ∇ × B along the ζ-edge computed using the left-and right-states of the µand ζ-edge magnetic field as shown in the right panel of Figure 9, and η A is a numerical resistivity dealing with the propagation of Alfvén waves near discontinuities.If computing the edge electric field using Equation (83) without the "resistive" term, simulations tend to develop spurious oscillations in MHD flows where Alfvén waves are dominant, especially when discontinuities occur.We should note that this "resistive" term η A j in Equation ( 83) only operates when the limiter detects a discontinuous situation.When the flow is smooth, B L µ = B R µ and B L ν = B R ν , the j ζ component in the diffusion term is essentially zero.Therefore in smooth regions the electric field E ζ is just a high-order approximation of −u × B computed at cell edges.
In Equation ( 83), the interpolated velocity u i+ 1 2 ,j+ 1 2 ,k and the average magnetic field B avg at the ζ-edge are not in the same coordinate system: the velocity u i+ 1 2 ,j+ 1 2 ,k is in the base Cartesian coordinate system (x, y, z), while the magnetic field B avg is in a local non-orthogonal coordinate system μinterp Then calculate the direction e 1 normal to the ζ vector: where ,k are the edge-centered positions A and B at two neighboring ζ-edges shown in Figure 11: then the third direction e 2 , orthogonal to both e 1 and e ζ , is obtained using the following cross product: (88) Since the reconstructed edge velocity u i+ 1 2 ,j+ 1 2 ,k is in the base Cartesian system, it is straightforward to map the reconstructed edge velocity vector u i+ 1 2 ,j+ 1 2 ,k into the new edge-aligned coordinate system (e 1 , e 2 , e ζ ): Figure 11: The edge-aligned coordinate system (e 1 , e 2 , e ζ ) for computing the electric field.
where u e1 and u e2 are the velocity components in the directions of e 1 and e 2 , respectively.The u e ζ component is along the ζ-edge which does not contribute to the calculation of the E ζ component.The edge magnetic field components B µ avg and B ν avg calculated using Equation ( 81) and ( 82) are in a non-orthogonal local coordinate system μinterp 2 ,k , ζ , thus the coordinate transform is more complicated than the velocity components, which requires solving the following 2 × 2 system: where ξ µ and ξ ν are the projections of the directions of magnetic field B µ,ν avg in the (e 1 , e 2 , e ξ ) coordinate system: After obtaining the average magnetic field components B e1 avg and B e2 avg in the new edge-aligned coordinate system (e 1 , e 2 , e ζ ), the electric field E ζ along the edge e ζ is calculated based on Equation ( 83): (95) Then the electric potential along the ζ-edge is calculated as which is the edge centered electric field E ζ multiplied by the length of the corresponding ζ edge.
2 ,k is the actual component used in the integral form of the Faraday's law to evolve the face-centered magnetic fluxes Φ.The diffusive speed v D in Eq. ( 95) is defined as The magnitude of the diffusive speed is usually set to reflect the magnitude of the local convection speed and Alfvén speed.The µ interp 2 ,k term originates from the ∇× operation in the coordinate transform of the non-orthogonal geometry when the diffusive current j ζ in Equation (83) is computed.

Evolution of Magnetic fluxes
The corresponding ν-edge and µ-edge aligned orthogonal coordinate systems are calculated in a similar way as shown in the ζ-edge.The electric potential E µ and E ν along the µ-edge and ν-edge are also computed in the same way using the edge-aligned coordinate system.These electric potential are calculated using the same subroutine as for E ζ through proper rotation of the computational grid.Once all the electric fields are calculated the magnetic flux threading a face can be updated via Faraday's law: Based on Equation (98), it is straightforward to show that the divergence of the magnetic field defined by the flux follows ensuring that the volume integrated ∇ • B is unmodified during the magnetic field evolution step, regardless of changes in the local electric fields.Therefore, as long as the initial divergence of the magnetic field is zero, the Maxwell solver keeps the ∇ • B term to round-off error automatically.
In this section, we show one set of test simulations of two-dimensional linear advection and four sets of simulation results from two-dimensional standard MHD test problems in both Cartesian and non-Cartesian geometries, including field-loop advection, circularly-polarized non-linear Alfvén waves, the Orszag-Tang vortex and spherical blast waves in strong magnetic fields.The 2-D linear advection tests demonstrate the the choice of high-order reconstruction methods discussed in Section 3.3.1.Both the field-loop advection and non-linear Alfvén wave simulations provide quantitative assessment for the numerical solutions, while the Orszag-Tang and blast wave simulations are less quantitative.The latter two sets of test simulations are used to demonstrate the effectiveness of the numerical schemes on handling highly non-linear MHD flows in non-orthogonal, distorted grid geometries.

Circular Advection
We use the circular advection of a slotted cylinder problem used by Colella et al. [35] as a multi-dimensional linear test to show the effectiveness of the advection scheme, which is originally developed by Zalesak.[36].This test problem is also used to demonstrate the necessity of high-order reconstruction method in the numerical algorithms.In this two-dimensional circular advection test, only the mass continuity equation is solved: where v 0 is a time-stationary circular velocity field defined as v 0 = −2πω θ with ω = 1 the angular velocity and r = (x − 0.5) 2 + (y − 0.5) 2 as shown in Figure 12a.The simulation domain is 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 with periodic boundary conditions in both the xand y-direction.The initial density distribution at t = 0 is defined as a slotted cylinder of radius R 0 = 0.15 centered at (x, y) = (0.5, 0.75), with slot width W = 0.06 and slot height H = 0.1: At t = 1, the center of the slotted cylinder returns to the initial location.The initial solution on a uniform Cartesian mesh with 128 × 128 cells is plotted in Figure 12a.The solution of ρ (at t = 1) using the second-order centered reconstruction method introduced in Section 3.3.1 is shown in Figure 12b, while the solution of ρ at t = 1 using the default seventh-order reconstruction scheme is shown in Figure 12c.At t = 1, although the second-order reconstruction scheme preserves the basic shape of the slotted cylinder, the initially sharp edges of the cylinder is smeared significantly with noticeable fill-in of the slot.When using the seventh-order reconstruction scheme, the density distribution at t = 1 is preserved more accurately compared to the second-order case, with sharper edges of the cylinder and an empty slot.
Figure 12: a) The initial distribution of a slotted cylinder in Cartesian geometry with 128 × 128 cells at t = 0; b) the simulated density distribution at t = 1 using a 2 nd -order reconstruction scheme; and c) the the simulated density distribution at t = 1 using a 7 th -order reconstruction scheme.
Figure 13: The comparison of density profiles at y = 0.75.The blue profile is the reference density distribution at t = 0, the green and the red profiles are the corresponding density distribution at t = 1 using the 2 nd -and 7 th -order schemes, respectively.
Figure 14: The circular advection results using three different curvilinear grids.
Figure 13 shows more quantitative comparisons of the circular advection results using line profiles of simulated ρ at y = 0.75.The blue profile shows the initial density distribution at t = 0 as a reference, while the green and the red profiles show the corresponding density profiles at t = 1 using the 2 nd -and the 7 th -order reconstruction scheme, respectively.The comparison shows that the 2 nd -order reconstruction scheme introduces a significant amount of numerical diffusion.At t = 1, the profile of the slotted density is smeared using the 2 nd -order reconstruction, which has approximately 14 cells resolving the contact discontinuity at the edge of the cylinder (x = 0.5 ± R 0 ) and a filled slot with density ≈1.3 (|x − 0.5| ≤ R 0 ).On the other hand, the 7 th -order reconstruction scheme preserves the initial profile well, with only 4 cells resolving the sharp edge of the cylinder and an empty slot (density ≈1.0).The comparison of the density profiles between the 2 nd -and the 7 th -order reconstruction scheme demonstrates the necessity of using very-high order reconstruction methods in order to reduce the numerical diffusion and resolve sharp contact discontinuity in multi-dimensional flow simulations.
Figure 15: The spatial distribution of magnetic energy at t = 4 from four curvilinear grids.

Non-linear Alfvén Wave
We run the non-linearly polarized circular Alfvén wave simulation from Tóth where B 0 = 1.0 is the strength of the initial magnetic field parallel to the wave vector.
Figure 17a shows the simulated B ⊥ component along the direction of wave propagation x at t = 5.0, using the default scheme with the 7 th -order spatial reconstruction and the PDM limiter.Six grid resolutions with N × N 2 cells are used to test the convergence of the solution at t = 5.0, by increasing the grid parameter N from 8 to 256.When the grid resolution is 64 × 32 and above, the solution of B ⊥ along x starts to converge to the analytical solution B ⊥ = 0.1 sin 2πx .Figure 17b shows the average L1 error norm as a function of grid size N using two different reconstruction methods (2 nd -and 7 th -order reconstruction), We find a 2 nd -order convergence rate based on this smooth but nonlinear Alfvén wave simulation, regardless of the size of the reconstruction stencil, suggesting that the formal order of convergence for the MHD solver is two.This is expected since no high-order quadrature methods are used in the evaluation of face fluxes in the finite-volume solver as discussed in Section 3.4.However, the formal order of convergence (2 nd -order) does not mean that the high-order reconstruction (7 th -order and above) is not necessary in the solver.In fact, it is important to note that when using the default 7 th -order reconstruction scheme, the magnitude of average L1 error is approximately a factor of ten lower than that from a 2 nd -order reconstruction scheme, suggesting the necessity of using high-order reconstructions for achieving highly-accurate solutions with moderate number of computational cells.This reduction in average L1 error is consistent with the comparisons in Figure 14, which shows significant improvement in preserving the shape of initial density structures.Note that when the non-clipping option is switched on, the average L1 error is further reduced by approximate a factor of 3 due to the preserving of local extrema.
Figure 18 compares the spatial distributions of B ⊥ at t = 5.0 from two sets of non-linear Alfvén wave simulations using different computational grids.Figure 18a is from a standard, orthogonal Cartesian grid, and Figure 18b is from a distorted, non-orthogonal Cartesian grid defined in Equation ( 102) -( 103) with w 0 = 0.1.Both simulations use 256 × 128 computational cells.Compared to the simulation using a Cartesian grid, the simulated wavefronts in the distorted, non-orthogonal grid are perpendicular to the direction of propagation, without noticeable distortions.This comparison shown in Figure 18 also illustrates the capability of the numerical schemes handling smooth MHD solutions in the non-linear regime when using non-orthogonal grids.

Orszag-Tang Vortex
We use the Orszag-Tang vortex simulation to test how robust the code is for handling the formation and interaction of non-linear MHD shocks in non-orthogonal, curvilinear geometry.The Orszag-Tang vortex is a well-known twodimensional problem for testing the transition to supersonic MHD turbulence, which is a very common test of numerical    MHD schemes used in many previous studies as a basis for consistent comparison of codes.In this test, we focus on the effectiveness of the numerical schemes in handling non-orthogonal grids and compare the results in non-orthogonal geometries with that from a standard uniform Cartesian grid.For simulating the Orszag-Tang vortex, we use a 2-D square simulation domain with 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1.The initial density ρ = 25  36 π and the initial pressure is P = 5  12 π distributed uniformly in the simulation domain, with γ = 5 3 .The initial velocities are set as v x = − sin(2πy) and v y = sin(2πx).The magnetic fluxes Φ are initialized using the following vector potential function A: where B 0 = 1.0.The above vector potential function gives B x = −B 0 sin(2πy) and B y = B 0 sin(4πx) initially.The boundary conditions are periodic everywhere.All the test simulations shown in Figure 19 have a computational grid with 512 × 512 cells.The reconstruction method is the default 7 th -order scheme without the non-clipping option switched on.
Figure 19a shows the simulated spatial distributions of plasma pressure at simulation time t = 0.48 in a uniform Cartesian grid, when the MHD shocks start to interact near the center of the simulation domain.Figure 19b-d shows the corresponding spatial distributions of plasma pressure at the same simulation time using three sets of non-orthogonal, distorted Cartesian grids as defined in Equation ( 102)-( 103) with increasing w 0 .Although these comparisons are less quantitative, they do suggest that the simulated distributions of plasma pressure are not significantly affected by the choice of the computational grid.Both smooth structures and shocks at t = 0.48 using non-uniform, non-orthogonal grids are remarkably similar to those in Figure 19a using a uniform Cartesian grid.The effectiveness of the MHD solver is also illustrated in more quantitative comparisons of using line profiles.Figure 20 shows the comparisons of the simulated plasma pressure profiles (at t = 0.48) along the y-direction with x = 0.5.Although the simulated pressure exhibits some minor differences at the shock transition regions, possibly due to the inhomogeneity of the grid, very little difference occurs in the smooth regions.The Orszag-Tang vortex simulations demonstrate the capability of the GAMERA code handling MHD shock propagation and interaction, without the requirement of the orthogonality of the computation geometry.An animation of the Orszag-Tang vortex simulation using the four different grid geometries is included in the supplemental material.

Spherical Blast Wave
We use MHD blast wave simulations to test the robustness the numeric schemes in handling multi-dimensional, strong MHD shocks and rarefactions.In a two-dimensional simulation domain in the x-y plane, The initial plasma has uniform density ρ = 1, thermal pressure P = 0.1 and zero velocity u = 0, with γ = 5 3 .Within the spatial region r = x 2 + y 2 < 0.1, the thermal pressure P is set to 10.0 (that is, 100 times greater than the ambient plasma pressure).The initial magnetic field is along the diagonal in the x-y plane with B x = 1 √ 2 , B y = 1 √ 2 and B z = 0. We use four different computational grids for the MHD blast wave simulation with the same initial conditions.These grids include a standard Cartesian grid (Figure 21a), a cylindrical polar grid (Figure 21b), a quarter of a ring (Figure 21c) and a distorted Cartesian grid (Figure 21d) as used in the Orszag-Tang simulations.Each grid has a resolution of 512 × 512 cells with the center of the blast wave located at the origin.The Cartesian grid is used as a reference, while the non-Cartesian grids are used to test the propagations of strong MHD shocks and rarefaction in non-orthogonal, curvilinear geometries.The cylindrical polar grid is singular at the pole, which is treated using a conservative averagingreconstruction technique developed by Zhang et al. [37] without changing the numerical schemes of the MHD solver.An animation of the MHD blast wave simulation in the four different grids is included in the supplemental material.As shown in the animation, at early simulation times (e.g., t < 0.05), the out-going blast waves are approximately spherical and show no grid alignment effects due to non-orthogonal curvilinear geometries.The results shown in Figure 21 are the spatial distributions of plasma density at t = 0.18.Although the spherical blast wave simulations are not very quantitative, it shows that the spatial distributions of plasma density in the blast wave simulation at a later time are not sensitive to the choice of the curvilinear grid.In the simulation results shown in Figure 21, both smooth structures and shocks in the non-Cartesian geometries (Figure 21b-d) are remarkably similar to that in the reference Cartesian grid (Figure 21a) without noticeable differences.Note that the Richtmyer-Meshkov instability is suppressed by the strong magnetic field, and no fingers are evident in the interior of the field-aligned bubble as those in a hydrodynamic blast wave simulation.

Summary
The GAMERA code is a re-invention of the LFM MHD kernel with significant upgrades in both numerical schemes and software implementation.The improvements in numerical schemes include 1) 12 th -order grid metric calculations using Gaussian quadrature; 2) high-order upwind reconstruction method; 3) the PDM limiter with an extremum-preserving option; 4) background magnetic field splitting for very low-β conditions.The Ring-Average technique for spherical axis singularity described by Zhang et al. [37] is also implemented in the GAMERA code for high-resolution simulations in general geometry with axis singularities.The improvements in software implementation include 1) modern Fortran implementation with minimum external library dependence; 2) core computational routines that act on vector-length aligned blocks of memory; 3) extensive use of OpenMP threading within a socket allowing MPI ranks distributed across sockets; 4) High performance Computing (HPC) -friendly data structures, namely the use of large, contiguous arrays and contiguous memory access patterns.The GAMERA code is designed to be easily applicable to multi-dimensional MHD flow simulations in non-orthogonal curvilinear grids adapted to specific conditions.Current applications of the GAMERA code includes global simulations of planetary magnetospheres, the inner heliosphere and solar wind, and local simulations of basic plasma physics applications, e.g., current sheets.
This paper provides a comprehensive description of the numerical recipes used in the GAMERA code.Our hope in the development of this paper has been to give some context for the key questions and solutions that need to be addressed in order to design a large scale simulation code for MHD problems.The numerical recipes described in this paper are useful as a reference, but it also shows that, when compared with other codes, these recipes are not the same.The nuances of disparate application areas will lead to different codes that excel in their tailored regimes.Although the numerical considerations of the LFM, inherited by GAMERA, are tailored to heliospheric environments we have designed GAMERA to be more general.We expect GAMERA to be a useful tool in any application where geometric flexibility is important.

Acknowledgement
The research was supported by the National Center for Atmospheric Research (NCAR) and the Johns Hopkins University Applied Physics Laboratory Independent Research & Development funds.We would like to acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX)provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation (NSF).Michael Wiltberger was serving at the NSF during the production of the this paper.Any opinion, findings, or conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.

A. Normalization of the ideal MHD equations
The mass and momentum equations in SI units follows where µ 0 = 4π × 10 −7 T • m/A.To normalize the MHD variables in the above equations, define where ρ 0 , u 0 , P 0 and B 0 are the normalization constants for mass, velocity, pressure and magnetic field, respectively.ρ, ũ, P , B are the corresponding dimensionless variables evolved by the MHD solver.For time and spatial coordinates, let t = t 0 t and x = x 0 x with t and ∇ the dimensionless temporal and spatial differential operators: First substitute (111) and (112) into the mass equation ( 109): If we choose the following normalization between time and velocity: the dimensionless mass equation is written as: Next substitute (111) -( 114) into the momentum equation ( 110) and given that x 0 = u 0 t 0 from the normalized mass equation: The above equation gives the following normalization for the plasma pressure P 0 and the magnetic field B 0 : Thus the dimensionless momentum equation is written as: It is straightforward to show that with the normalization relations (118), ( 122) and (123), the dimensionless plasma energy equation and Faraday's law are written as: where ẼP = 1 2 ρũ 2 + P γ−1 is the normalized plasma energy.To summarize, Equations ( 119), ( 124), ( 125) and (126) forms a set of dimensionless MHD equations with the following normalizations:

B. The Grid Metric Calculations
As described in Section 3.1, the primary grid is the cell corner coordinates x i± 1 2 ,j± 1 2 ,k± 1 2 .Thus the corresponding cell center positions x C i,j,k for cell (i, j, k) as shown in Figure 1 are computed using the eight cell corners forming a hexahedron as: In the GAMERA code, the face center locations x µ , x ν and x ζ are used directly in the calculation of magnetic field vectors.Thus they are calculated using high-order Gaussian quadrature, which gives better definitions of the "barycenter" of a cell interface especially when the cells are either very distorted or degenerated to tetrahedrons, e.g., spherical grid cells near the axis.For example, when using two-dimensional Gaussian quadratures, the µ-face center locations x µ i+ 1 2 ,j,k are calculated as:

C. The General Form of the PDM Operator
For non-linear problems, the form of the PDM operator M + using a 2 nd -order centered reconstruction is written as [14] M where and A is the parameter determines the amount of numerical diffusion.The direction of sweeping is sign(v), thus Equation (138) gives both the left-and right-state at interface i + 1 2 .If A = 0, the above operator M + becomes the donor cell method which is one of the most diffusive first-order method.For linear problems, Hain [14] have shown that the monotonicity condition for the PDM parameter A is related to the Courant number N CF L : In principle any value of A > 0 could be used in the PDM limiter.Larger values of A lead to more aggressive limiting as indicated by (138).However, the allowable Courant number decreases with increasing A as shown in equation ( 139).
In practice, we find that there is little improvement in resolving square wave profiles if A > 4. Thus in the GAMERA solver, the default value for the PDM value A is chosen to be 4.0, together with a fixed Courant number N CF L = 0.3 which satisfies (139).The PDM limiter can be extended to use arbitrary high-order spatial reconstruction instead of using m + in Equation (138) as shown in Figure 22.For a high-order interface value f HO reconstructed at i + 1 2 , the left-state (v > 0) at interface i + 1 2 is calculated using the following equation: Figure 22: The calculation of left-state at interface i + 1 2 using the PDM limiter.
Here the f HO is a "clipped" version of the high-order reconstruction value calculated as f HO = median(f i , f HO , f i+1 ).
It is very straightforward to show that if A = 0, the above equation becomes which is the first-order Donor cell method.The numerical diffusion is reduced significantly as A is increased from 0 to 4.0.

D. The Choice of Reconstruction Order
The PDM limiter can be combined with an arbitrarily high order reconstruction scheme as discussed in Section 3.3.1.
In the GAMERA code, the default choices of the reconstruction method is the 7 th -order upwind reconstruction (the 8 th -order centered reconstruction is the one used in the original LFM, which is also implemented in the GAMERA code for reference).The reason for choosing such high-order reconstruction schemes as the default is based on achieving low numerical diffusion with reasonable amount of computing resource, which is important for 3-D global-scale simulations of large space plasma systems such as planetary magnetospheres and the heliosphere.Here we show two sets of numerical solutions of the linear advection equation to justify the necessity of the high-order reconstruction implemented in the GAMERA code.

1-D Linear Advection of Four Shapes in Non-uniform Grid
The 1-D linear advection equation solved in this section follows where u = 1.0 and x is defined on [−1, 1].The computational domain is non-uniform with where i = 0, 1, 2, ...N , and N is the total number of cell centers in the x-direction.Periodic boundary conditions are imposed at x = 1 and x = −1, thus the initial profile of f goes back to exactly the same location at t = 2.The red profile in the top panel of Figure 23 shows the initial profile of f at t = 0, with the gray vertical lines showing variation of the non-uniform grid geometry.Four distinctive shapes are used in the initial profile: a Gaussian peak, a square wave, a triangle and a half-circle.This linear advection test simulation has been used extensively in previous studies to test the quality of advection schemes in fluid simulations.
As shown in the initial profile, smooth resolution is tested by a narrow Gaussian wave with a width of 2.4∆x centered at x = −0.75, while discontinuity resolution is tested by a square wave centered at x = −0.25.Discontinuity in the first derivative is tested using a triangular profile centered at x = 0.25.Although there is no discontinuity in value, the half-circle profile centered at x = 0.75 is very challenging since it combines sudden and gradual changes in gradient with limited number of cells.In addition, there is a discontinuity in curvature at each side of the half-circle profile, which makes it more challenging to resolve the profile especially in non-uniform geometry.In the following simulations, the PDM parameter A described in Appendix C is chosen to be 4.0 and the Courant number is 0.3, which are the same as the default values used in the MHD solver.The black profiles in Figure 23 are the simulated distributions of f at t = 2.0 using six different orders of upwind reconstruction (e.g., 1 st , 3 rd , ..., 11 th -order).The 1 st -order scheme, which is basically the Donor cell method, is over-diffusive -all the four profiles are smeared significantly as Error functions at t = 2.0.As the order of reconstruction increases to 7 th -order, the algorithm starts to resolve the full narrow Gaussian peak within 9 cells.Above 7th-order, the improvement on the Gaussian peak is small.The resolution of the square wave is also improved as the order of reconstruction increases.Above 5 th -order, the resolution of the sharp transition of the right-side of the square wave remains at 4 cells.In other words, the improvement on resolving the square wave is small.Given that the total computing time (including both reconstruction and limiting) of the 7 th -order scheme is only about 10% more than the 5 th -order scheme, we choose the 7 th -order scheme as the default choice for spatial reconstruction in the GAMERA code.A more quantitative analysis on the effective numerical diffusion coefficient based on the advection of a square wave is in the next section.

The Effective Diffusion Coefficient
When the one-dimensional linear advection equation ( 142) is solved numerically, the effective equation solved by numerical schemes is written as where u is the advection speed and D is an effective diffusion coefficient, which depends on both the numerical scheme and the grid resolution.In general, the above equation is a mixed advection-diffusion equation with an analytical solution written as When the diffusion coefficient D ≡ 0, the solution becomes We use the numerical solution to a special case of the above linear advection-diffusion equation to approximately quantify the amount of effective numerical diffusion introduced by different reconstruction methods.If the linear advection-diffusion equation ( 144) is solved numerically in a semi-infinite domain x ∈ [0, +∞) with an initial condition f (x, t = 0) = 0 and a boundary condition f (x = 0, t) = 1, the analytical solution (145) of the advection-diffusion equation with constant u and D is written as: where erfc(•) is the complimentary error function.After solving the linear advection equation (142) numerically with f (x, t = 0) = 0 and f (x = 0, t) = 1, the numerical diffusion coefficient D is approximated by fitting the numerical solution f i to the analytical solution (147).To simplify the calculations, we solve the linear advection equation (142) in a uniform computational domain defined at [0, 4] with 256 cells.At t = 2, the front of the step function is located at x = 2, and a least square fit process is applied to the numerical solution f i to estimate the effective diffusion coefficient D. Figure 24 shows the numerical solutions of the step function using reconstruction schemes from 1 st -order to 12 th -order, together with the effective diffusion coefficient √ D calculated for each profile.The numerical diffusion width W listed in each panel is defined as the number of computational cells in the transition region of the step function (0.01 < f i < 0.99).If the effective diffusion coefficient D is exactly zero, the numerical diffusion width W is expected to be 0. The summary of the numerical solution to the advection equation ( 142) is listed in Table 6.Reconstruction schemes from 1 st -order to 12 th -order are tested with the analytical solution to the advection-diffusion equation (144).The effective diffusion coefficient ( √ D) is listed in the second column, which decreases from 0.0604 to 0.0051 as the order of reconstruction increases from 1 to 12.The numerical diffusion width (W ) is listed in the third column of Table 6.When the order of reconstruction is higher than 7 th -order, the numerical diffusion width remains at 4, suggesting that the improvement on the slope of the square wave is small above 7 th -order.This quantitative comparison is consistent with the linear advection test simulations shown in Figure 23.The third column shows the relative grid Reynolds number in the simulation normalized by the effective grid Reynolds number (uL/D) in the 1 st -order solution.When using a 7 th -order reconstruction scheme, the effective grid Reynolds number is enhanced by approximately two orders of magnitude compared to the 1 st -order solution.Above 7 th -order, there is still improvement to the relative grid Reynolds number but the improvement is small.If we use a simple measure (R/Order) to quantify the efficient of the reconstruction methods, it peaks at 7 th -order as shown in the last column of Table 4.

Step 3 .Step 4 .Step 5 .
Split the provisional values into left-and right-interface states (U Evaluate numerical flux F(•) and M(•) at cell interfaces using left-and right-states; Integrate numerical flux F and M through µ-interfaces and add to ∆U F and ∆U B ; Step 6. Repeat Step 2-5 in the νand ζ-direction.

Figure 5 :
Figure 5: The calculation of the PDM value on the left-side of the interface i + 1 2 , using the advection of density in a one-dimensional stencil towards the right.The figure is adapted from Huba [2003].

Figure 6 :
Figure 6: A local maximum versus numerical peaks within a stencil of five computational cells.

Figure 8 :
Figure 8: The algorithm for calculating edge-aligned electric field using the predictors.

i+ 1 2 ,j+ 1 2 1 2 ,j+ 1 2
,k , ζ as illustrated by the green arrows in Figure 11.Therefore to calculate the cross product between the velocity and magnetic field at cell edges, we define a new orthogonal coordinate system (e 1 , e 2 , e ζ ) with one axis e ζ aligned with the direction of the cell edge and transform u i+ ,k and B µ,ν avg into this new coordinate system.Using the ζ-edge aligned orthogonal system shown in Figure 11 as an example, the local coordinate system (e 1 , e 2 , e ζ ) for computing the ζ-component of the electric field using Equation (83) is shown in red, where E ζ is along the direction of e ζ .For the ζ-edge indexed by i + 1 2 , j + 1 2 , k the (e 1 , e 2 , e ζ ) coordinate system is computed as follows.First compute the ζ-edge unit vector e ζ :

[ 4 ]
to test the performance of the numerical MHD schemes in nonlinear regime.Since the Alfvén wave simulation is smooth, it is also a reasonable test for convergence of the MHD solver.The simulation domain is chosen to be a 2-D Cartesian box with 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1 cos α , where α = π 3 is the angle at which the Alfvén wave propagates with respect to the x-axis.Thus the waves do not propagate along the diagonal of the simulation domain, which guarantees the multi-dimensional nature of the test simulation by making the numerical flux different in the xand y-direction.The initial conditions are ρ = 1, P = 0.1, u ⊥ = 0.1 sin 2πx , B ⊥ = 0.1 sin 2πx , and u z = B z = 0.1 cos 2πx with γ = 5 3 and x = (x cos α + y sin α), where u ⊥ and B ⊥ are the components of velocity and magnetic field perpendicular to the wavevector.The B and B ⊥ components are calculated using cell centered B x and B y via B ⊥ = B y cos α − B x sin α, and B = B x cos α + B y sin α.The initial magnetic fluxes Φ are calculated using the following vector potential A:

Figure 16 :
Figure 16: The decay of magnetic energy as a function of time derived from the four field-loop advection simulations using different grids.

Figure 17 :
Figure 17: a) The distribution of B ⊥ versus x at t=5.0, computed from five simulations with increasing resolution.b) the average L1 error as a function of grid resolution N .

Figure 18 :
Figure 18: The spatial distribution of B ⊥ at t = 5.0 derived from two simulations in a) Cartesian geometry and b) distorted Cartesian geometry with 256 × 128 cells.

Figure 19 :
Figure 19: The spatial distribution of plasma pressure P at t = 0.48 in four Orszag-Tang simulations using different grids.Each simulation has 512 × 512 cells.

Figure 20 :
Figure 20: The line profiles of plasma pressure at x = 0.5 from four Orszag-Tang simulations using different grids, taken from same simulation snapshots with t = 0.48.

Figure 21 :
Figure 21: The spatial distribution of plasma pressure P at t = 0.48 in four MHD spherical blast wave simulations using different grids.Each simulation has 512 × 512 cells.

Figure 23 :
Figure 23: 1-D linear advection of four shapes in a non-uniform grid using different orders of reconstruction.The vertical gray line shows the spatial distribution of the computational grid.The red profile is the initial condition at t = 0.The black profiles are simulation results at t = 2.0 with different orders of spatial reconstruction.

1 2 Table 1 :
The grid definitions and index notations used in the MHD solver.

Table 2 :
The variable definitions in the MHD solver.