Divergence-Free Magnetohydrodynamics on Conformally Moving, Adaptive Meshes Using a Vector Potential Method

TL;DR

This paper introduces a novel divergence-free magnetohydrodynamics method that employs a vector potential to evolve magnetic fields on adaptive, conformally moving meshes, ensuring high accuracy and maintaining divergence-free conditions.

Contribution

The paper presents a new vector potential-based approach for MHD simulations on adaptive, moving meshes, improving divergence control and computational flexibility.

Findings

Maintains divergence-free magnetic fields to machine precision.

Successfully tested on a variety of MHD problems, including relativistic cases.

Produces accurate results across different grid configurations.

Abstract

We present a new method for evolving the equations of magnetohydrodynamics (both Newtonian and relativistic) that is capable of maintaining a divergence-free magnetic field ($\nabla \cdot \mathbf{B} = 0$) on adaptively refined, conformally moving meshes. The method relies on evolving the magnetic vector potential and then using it to reconstruct the magnetic fields. The advantage of this approach is that the vector potential is not subject to a constraint equation in the same way the magnetic field is, and so can be refined and moved in a straightforward way. We test this new method against a wide array of problems from simple Alfven waves on a uniform grid to general relativistic MHD simulations of black hole accretion on a nested, spherical-polar grid. We find that the code produces accurate results and in all cases maintains a divergence-free magnetic field to machine precision.

Figures

Tables

Table 1. Table 1 : Steps taken during a single evolution step of Cosmos++ .

Evolution algorithm

1. calculate timestep (a) Each physics package is polled for its smallest timestep (b) All packages, in all zones, on all mesh levels, are advanced using the smallest overall timestep 2. move mesh 3. refine mesh 4. second-order time-stepping loop (a) invert conserved fields to primitives i. Use current $𝐀$ to recover $𝐁$ (b) HRSC update i. construct face-centered fluxes ${𝐅}_R$ and ${𝐅}_L$ A. Use magnetic fluxes to construct emfs used to update $𝐀$ ii. update conserved fields from fluxes iii. update conserved fields from source terms 5. final reconstruction of updated conserved fields

Table 2. Table 2 : Refinement algorithms in Cosmos++ .

Strict refinement algorithm	Fast refinement algorithm
tag zones based on refinement criteria	tag zones based on refinement criteria
for each level $i$ of refinement, from highest to lowest	for each zone $j$ tagged for refinement
for each zone $j$ tagged for refinement	doRefine = true
for each neighbor zone $k$	for each neighbor zone $k$
if $j$ is already more refined than $k$	if $j$ is already more refined than $k$
refine($k$)	doRefine = false
refine($j$)	break
	if doRefine == true
	refine($j$)

Table 3. Table 3 : L 𝐿 L -1 Norm Errors for Circularly Polarized Alfvén Wave

$N$	${\Vert E\left(D\right)\Vert }_1$	${\Vert E\left(ℰ\right)\Vert }_1$	${\Vert E\left({ℬ}^x\right)\Vert }_1$	${\Vert E\left({ℬ}^y\right)\Vert }_1$	${\Vert E\left({ℬ}^z\right)\Vert }_1$
Traveling
8	$9.07\times {10}^{-3}$	$2.57\times {10}^{-3}$	$1.77\times {10}^{-2}$	$2.33\times {10}^{-2}$	$2.41\times {10}^{-2}$
16	$1.84\times {10}^{-3}$	$5.30\times {10}^{-4}$	$3.98\times {10}^{-3}$	$5.57\times {10}^{-3}$	$5.58\times {10}^{-3}$
32	$4.31\times {10}^{-4}$	$1.29\times {10}^{-4}$	$9.64\times {10}^{-4}$	$1.36\times {10}^{-3}$	$1.36\times {10}^{-3}$
64	$1.06\times {10}^{-4}$	$3.21\times {10}^{-5}$	$2.39\times {10}^{-4}$	$3.39\times {10}^{-4}$	$3.38\times {10}^{-4}$
128	$2.64\times {10}^{-5}$	$8.00\times {10}^{-6}$	$5.97\times {10}^{-5}$	$8.46\times {10}^{-5}$	$8.43\times {10}^{-5}$
Standing
8	$4.40\times {10}^{-3}$	$2.37\times {10}^{-3}$	$2.28\times {10}^{-3}$	$5.92\times {10}^{-3}$	$5.96\times {10}^{-3}$
16	$8.29\times {10}^{-4}$	$4.74\times {10}^{-4}$	$4.42\times {10}^{-4}$	$1.09\times {10}^{-3}$	$1.09\times {10}^{-3}$
32	$1.92\times {10}^{-4}$	$1.11\times {10}^{-4}$	$1.01\times {10}^{-4}$	$2.41\times {10}^{-4}$	$2.41\times {10}^{-4}$
64	$4.69\times {10}^{-5}$	$2.72\times {10}^{-5}$	$2.46\times {10}^{-5}$	$5.85\times {10}^{-5}$	$5.83\times {10}^{-5}$
128	$1.17\times {10}^{-5}$	$6.76\times {10}^{-6}$	$6.12\times {10}^{-6}$	$1.45\times {10}^{-5}$	$1.45\times {10}^{-5}$

Table 4. Table 4 : L 𝐿 L -1 Norm Errors for MHD Vortex

$N$	${\Vert E\left(D\right)\Vert }_1$	${\Vert E\left(E\right)\Vert }_1$	${\Vert E\left(s_x\right)\Vert }_1$	${\Vert E\left(s_y\right)\Vert }_1$	${\Vert E\left(B^x\right)\Vert }_1$	${\Vert E\left(B^y\right)\Vert }_1$
50	$1.54\times {10}^{-4}$	$7.43\times {10}^{-4}$	$5.39\times {10}^{-4}$	$5.55\times {10}^{-4}$	$2.16\times {10}^{-3}$	$1.58\times {10}^{-3}$
100	$4.01\times {10}^{-5}$	$1.79\times {10}^{-4}$	$1.34\times {10}^{-4}$	$1.32\times {10}^{-4}$	$5.37\times {10}^{-4}$	$4.00\times {10}^{-4}$
200	$1.01\times {10}^{-5}$	$4.47\times {10}^{-5}$	$3.36\times {10}^{-5}$	$3.24\times {10}^{-5}$	$1.34\times {10}^{-4}$	$1.00\times {10}^{-4}$
400	$2.54\times {10}^{-6}$	$1.12\times {10}^{-5}$	$8.41\times {10}^{-6}$	$8.08\times {10}^{-6}$	$3.35\times {10}^{-5}$	$2.52\times {10}^{-5}$

Table 5. Table 5 : L 𝐿 L -1 Norm Errors for MHD Vortex with Mesh Motion

$N$	${\Vert E\left(D\right)\Vert }_1$	${\Vert E\left(E\right)\Vert }_1$	${\Vert E\left(s_x\right)\Vert }_1$	${\Vert E\left(s_y\right)\Vert }_1$	${\Vert E\left(B^x\right)\Vert }_1$	${\Vert E\left(B^y\right)\Vert }_1$
50	$1.02\times {10}^{-4}$	$1.26\times {10}^{-3}$	$8.00\times {10}^{-4}$	$7.96\times {10}^{-4}$	$1.57\times {10}^{-3}$	$1.35\times {10}^{-3}$
100	$2.67\times {10}^{-5}$	$3.27\times {10}^{-4}$	$2.07\times {10}^{-4}$	$2.08\times {10}^{-4}$	$4.00\times {10}^{-4}$	$3.43\times {10}^{-4}$
200	$6.74\times {10}^{-6}$	$8.23\times {10}^{-5}$	$5.23\times {10}^{-5}$	$5.25\times {10}^{-5}$	$1.01\times {10}^{-4}$	$8.64\times {10}^{-5}$

Equations

$${\mathcal{T}}^{\alpha \beta }=\rho v^{\alpha }{v}^{\beta }+(P+P_B)g^{\alpha \beta }-b^{\alpha }{b}^{\beta }\text{leavevmode},$$

$${\partial }_i{v}^i=\frac{\partial v^1}{\partial x^1}+\frac{\partial v^2}{\partial x^2}+\frac{\partial v^3}{\partial x^3}\text{leavevmode}.$$

$${\partial }_t(\sqrt{g}\rho )+{\partial }_i\left[\sqrt{g}\rho (v^i-V_g^i)\right]+\sqrt{g}\rho {\partial }_i{V}_g^i=0\text{leavevmode},$$

$${\partial }_t(\sqrt{g}E)+{\partial }_i\left\{\sqrt{g}E(v^i-V_g^i)+\sqrt{g}\left[(P+P_B)v^i-b^i{b}_j{v}^j\right]\right\}+\sqrt{g}E{\partial }_i{V}_g^i=-\sqrt{g}\rho v^i{\partial }_i\varphi \text{leavevmode},$$

$${\partial }_t(\sqrt{g}{s}_j)+{\partial }_i\left\{\sqrt{g}{s}_j(v^i-V_g^i)+\sqrt{g}\left[(P+P_B)g^{ij}-b^i{b}^j\right]\right\}+\sqrt{g}{s}_j{\partial }_i{V}_g^i=\sqrt{g}{\mathcal{T}}^{ik}{\Gamma }_{ikj}-\sqrt{g}\rho {\partial }_j\varphi \text{leavevmode},$$

$${\partial }_t{B}^j+{\partial }_i\left[B^j(v^i-V_g^i)-B^i{v}^j\right]+B^j{\partial }_i{V}_g^i=0\text{leavevmode},$$

$$T^{\alpha \beta }=(\rho h+2P_B/c^2)u^{\alpha }{u}^{\beta }+(P+P_B)g^{\alpha \beta }-b^{\alpha }{b}^{\beta }\text{leavevmode}.$$

$${\partial }_{t}D+{\partial }_i\left[D(V^i-V_g^i)\right]+D{\partial }_i{V}_g^i=0\text{leavevmode},$$

$${\partial }_t\mathcal{E}+{\partial }_i\left[\mathcal{E}(V^i-V_g^i)\right]+{\partial }_i\left[\sqrt{-g}(P+P_B)V^i\right]+\mathcal{E}{\partial }_i{V}_g^i=-\sqrt{-g}{T}_{\sigma }^{\mu }{\Gamma }_{\mu 0}^{\sigma }\text{leavevmode},$$

$${\partial }_t{\mathcal{S}}_j+{\partial }_i\left[{\mathcal{S}}_j(V^i-V_g^i)\right]+{\partial }_i\left[\sqrt{-g}(P+P_B)\text{leavevmode}{g}_j^0\text{leavevmode}{V}^i\right]+{\mathcal{S}}_j{\partial }_i{V}_g^i=\sqrt{-g}{T}_{\sigma }^{\mu }{\Gamma }_{\mu j}^{\sigma }\text{leavevmode}.$$

$${\partial }_t{\mathcal{B}}^j+{\partial }_i\left[{\mathcal{B}}^j(V^i-V_g^i)-{\mathcal{B}}^i{V}^j\right]+{\mathcal{B}}^j{\partial }_i{V}_g^i=0\text{leavevmode},$$

$${\partial }_t\mathbf{U}(\mathbf{P})+{\partial }_i{\mathbf{F}}^i(\mathbf{P})=\mathbf{S}(\mathbf{P})$$

$${\int }_V{\partial }_t\mathbf{U}dV=-{\oint }_S{\mathbf{F}}^{i}d{A}_i+{\int }_V\mathbf{S}dV\text{leavevmode}.$$

$${\mathbf{U}}^{n+1}={\mathbf{U}}^n-\frac{\Delta t}{V}\sum _{faces}\left({\mathbf{F}}^i{A}_i\right)+\Delta t\mathbf{S}\text{leavevmode},$$

$$\mathbf{F}=\frac{c_{min}{\mathbf{F}}_R+c_{max}{\mathbf{F}}_L-c_{max}{c}_{min}({\mathbf{U}}_R-{\mathbf{U}}_L)}{c_{max}+c_{min}}\text{leavevmode},$$

$${\partial }_t(\sqrt{-g}{T}_{\alpha }^0)+{\partial }_i(\sqrt{-g}{T}_{\alpha }^0{V}^i)$$

$${\partial }_t(\sqrt{-g}{T}_{\alpha }^0)+{\partial }_i(\sqrt{-g}{T}_{\alpha }^0(V^i-V_g^i))+\sqrt{-g}{T}_{\alpha }^0{\partial }_i{V}_g^i\text{leavevmode}.$$

$$\frac{\partial \vec{\mathcal{A}}}{\partial t}-({\mathbf{V}}_g\cdot \nabla )\vec{\mathcal{A}}$$

$${\partial }_t{\mathcal{A}}_i-V_g^j{\partial }_j{\mathcal{A}}_i$$

$$=-{\mathcal{E}}_i+{\partial }_i\left[\frac{1}{g^{00}}(\Phi -g^{0j}{\mathcal{A}}_j)\right]\text{leavevmode},$$

$${\left({\mathcal{E}}_3\right)}_{i-1/2,j-1/2,k}=\frac{1}{4}\left[{\left(F_{B^1}^2\right)}_{i,j-1/2,k}-{\left(F_{B^2}^1\right)}_{i-1/2,j,k}+{\left(F_{B^1}^2\right)}_{i-1,j-1/2,k}-{\left(F_{B^2}^1\right)}_{i-1/2,j-1,k}\right]\text{leavevmode},$$

$$\Phi =g^{0j}{\mathcal{A}}_j$$

$${\partial }_t{\mathcal{A}}_i-V_g^j{\partial }_j{\mathcal{A}}_i=-{\mathcal{E}}_i\text{leavevmode}.$$

$${\partial }_t(\sqrt{-g}\Phi )-V_g^i{\partial }_i(\sqrt{-g}\Phi )+{\partial }_i\left(\sqrt{-g}{\mathcal{A}}^i+\frac{g^{0i}}{g^{00}}\sqrt{-g}\Phi \right)=0\text{leavevmode}.$$

$$\mathbf{B}$$

$${\mathcal{B}}^i$$

$$(\mathrm{Newtonian})$$

$$B^1$$

$$B^2$$

$$B^3$$

$$\nabla \cdot \mathbf{B}=\nabla \cdot (\nabla \times \mathbf{A})=0\text{leavevmode},$$

$${\partial }_i{\mathcal{B}}^i$$

$$B_{i-\frac{1}{2},j,k}^1$$

Full text

††journal: Journal of Computational Physics

1 Introduction

In magnetohydrodynamics (MHD), the time component of Maxwell’s equations, known independently as Gauss’ law of magnetism, provides a constraint condition. Mathematically, the constraint is $\nabla\cdotp\mathbf{B}=0$, where $\mathbf{B}$ is the magnetic field vector. Simply put, this says that a naturally occurring magnetic field should be divergence free. Physically, this is equivalent to saying that there are no magnetic monopoles.

However, most numerical implementations of MHD are not guaranteed to abide by the $\nabla\cdotp\mathbf{B}=0$ constraint, even if the initial field is divergence free. Thus, the magnetic field is allowed to build up unphysical divergence (monopoles) during the evolution. This can potentially lead to multiple numerical problems, including nonlinear instabilities, potentially driving pressures and densities to unphysical values [1]. More critically, for conservative schemes, the build up of $\nabla\cdot\mathbf{B}$ can generate anomalous forces parallel to the field [2] or cause an extra compressive component to arise within the magnetic field itself [3].

Several approaches for controlling the build up of $\nabla\cdotp\mathbf{B}$ have been proposed [e.g. 2, 4, 5, 6, 3, 7]. Some aim only to minimize it by adding supplemental fields that dampen or otherwise decay the accumulation of monopoles [8, 9]. Although straightforward to implement in some numerical schemes, these approaches formally break the conservative nature of the MHD equations [10], which may be undesirable for certain schemes.

Other approaches [4, 5, 6, 3] aim to carefully evolve the magnetic field in such a way that after each update cycle, it is again divergence free, at least for one chosen numerical representation of the divergence [see 7, for an exploration of some of the different implementations for Godunov codes]. One such approach is known as Constrained Transport (CT), so named because of its aim to constrain the evolution of the magnetic fields. While many of these approaches work well for uniform, static meshes, they can be challenging to extend to unstructured or refined meshes [10]. However, Balsara & Dumser [11] have recently discovered a projection operator that can be coupled with a cell-centered WENO algorithm to preserve a divergence-free magnetic field on unstructured meshes. Additional work [12, 13] has focused on finding suitable multi-dimensional Riemann solvers to obtain edge-centered electric fields, which are then used to update face-centered representations of the magnetic field.

In this paper we take a different approach, explored in some previous codes [e.g. 14, 15], to evolve the magnetic vector potential rather than the magnetic field itself. Since the magnetic field can then be reconstructed as needed by taking the curl of the vector potential, and since the divergence of a curl is always zero, this procedure is guaranteed to maintain a divergence-free field. The benefit of this approach for our purposes is that, since there are no constraints on the vector potential, it can be refined, de-refined, and relocated without introducing divergence errors.

There has been some concern, though, expressed in the literature about the use of the vector potential (VP) approach, mainly that it may produce spurious forces or current reversals near discontinuities [4]. The trouble comes from the calculation of the Lorentz force, $\mathbf{J}\times\mathbf{B}=(\nabla\times\mathbf{B})\times\mathbf{B}/4\pi$. The thinking is that, because the Lorentz force requires taking a derivative of the magnetic field (to get the electric current), that equates to a second derivative of the vector potential. While Choptuik [16] showed that the order of accuracy of derivatives is the same as the order of the underlying reconstruction of the field for smooth functions, so that the order of accuracy of the method can be preserved, the concern is in how the VP method behaves when the field is not smooth, i.e. near discontinuities. The expectation of Evans & Hawley [4] was that these methods would produce spurious current reversals in such locations. However, as we show, we do not see any evidence of spurious currents nor forces in our tests.

This paper will focus on the development of our magnetic VP method, specifically as implemented in the Cosmos++ computational astrophysics code [17]. The paper is structured as follows: Section 2 provides an overview of the capabilities and general structure of Cosmos++; Section 3 details the implementation of the new VP method within Cosmos++; Section 4 presents the results of well known MHD and astrophysical test problems generated using the updated code; and conclusions are drawn in Section 5.

2 Cosmos++

Cosmos++ is a parallel, multi-dimensional, fully covariant, modern, object-oriented (C++), radiation magnetohydrodynamics (RMHD) code for both Newtonian and general relativistic astrophysical and cosmological applications. Cosmos++ utilizes unstructured meshes with adaptive ($h$-) refinement [17], moving-mesh ($r$-refinement) [18], and adaptive order ($p$-refinement) [19] capabilities, enabling it to evolve fluid systems over a wide range of spatial scales with targeted precision. Cosmos++ can, in principle, support triangular/tetrahedral meshes, though it currently only offers quadrilateral/hexahedral options. Nevertheless, they are unstructured in the sense that Cosmos++ supports arbitrarily reduced and enhanced nodal connectivity. The mesh refinement is done on a local, cell-by-cell basis as illustrated in Fig. 2.1 for a simple example of a 3-level mesh tracking the dense gas front created by a cylindrical blast wave plowing into an ambient medium. A typical block-refinement algorithm would be much less efficient on problems such as this. Cosmos++ can also use adaptive and moving mesh refinement [20] or adaptive mesh and adaptive order [19] in conjunction for even greater precision.

2.1 Equations of Magnetohydrodynamics

Since Cosmos++ is designed to handle problems in both Newtonian and general relativistic gravity, we present the MHD equations for both, though they are intentionally written to look quite similar.

2.1.1 Newtonian

Introducing an effective Newtonian MHD stress energy tensor

[TABLE]

where Greek indices run over the four spacetime coordinates $\{0,1,2,3\}$, we can derive the four fluid equations (energy and the three components of momentum) from the conservation of stress energy: $\nabla_{\mu}{\cal T}^{\mu}_{\ \nu}=\partial_{\mu}{\cal T}^{\mu}_{\ \nu}+\Gamma^{\mu}_{\alpha\mu}{\cal T}^{\alpha}_{\ \nu}-\Gamma^{\alpha}_{\mu\nu}{\cal T}^{\mu}_{\ \alpha}=S_{\nu}$, where $\rho$ is the fluid mass density, $v^{i}$ is the fluid velocity ($v^{0}=1$), $P$ is the fluid pressure (for an ideal gas $P=(\Gamma-1)e$, $e=\rho\epsilon$ is the fluid internal energy density, $\Gamma$ is the adiabatic index), $P_{B}=b^{i}b_{i}/2$ is the magnetic pressure, $b^{i}$ is the primitive Newtonian field ($b^{0}=0$), $g_{\alpha\beta}$ is the curvature metric used here to handle curvature associated with non-Cartesian meshes, $\Gamma^{\alpha}_{\mu\nu}$ are the Christoffel symbols representing the metric derivatives, and $S_{\nu}$ includes any source terms. Note that we define our magnetic field such that the factor of $\sqrt{4\pi}$ is absorbed. In addition to energy and momentum, we also require equations for the conservation of mass $\nabla_{\mu}(\rho v^{\mu})=0$ and magnetic induction $\nabla_{\mu}(v^{\mu}b^{\nu}-b^{\mu}v^{\nu})=0$. Note that here and throughout, we use the Einstein summation convention, where repeated indices imply summation, e.g.

[TABLE]

Expanding out the space and time coordinates and including a grid velocity $V_{g}^{i}$ for when moving mesh is active, the Newtonian conservation equations (of mass, energy, and momentum) and the magnetic induction equation can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $s_{i}$ are the covariant momentum components, $B^{i}=\sqrt{g}b^{i}$ is the conserved Newtonian magnetic field, $\sqrt{g}$ is the determinant of the curvature metric, $\phi$ is the gravitational potential, and $E$ is the total energy density, which includes internal, magnetic, and kinetic energy contributions: $E=e+b^{i}b_{i}/2+\rho v^{i}v_{i}/2$. By treating the gravitational energy as a separate source term, the Newtonian and relativistic flux and source terms are easily interchangeable within the code framework. See [21] for a more complete description of the Newtonian capabilities of Cosmos++.

2.1.2 General Relativity

For general relativistic MHD, we follow a very similar procedure to arrive at the appropriate conservation equations. We again start from the contravariant stress energy tensor for an ideal MHD fluid:

[TABLE]

Here $h=1+\epsilon/c^{2}+P/(\rho c^{2})$ is the specific enthalpy, $c$ is the speed of light, $u^{\alpha}=u^{0}V^{\alpha}$ is the contravariant four-velocity, $V^{\alpha}$ is the transport velocity, $b^{\alpha}$ is the magnetic field measured by an observer co-moving with the fluid, and $P_{B}=b^{\alpha}b_{\alpha}/2$.

By defining the boosted mass density, total energy density, momentum density, and magnetic fields as $D=\sqrt{-g}u^{0}\rho=W\rho$, ${\cal E}=-\sqrt{-g}{T}^{0}_{0}$, ${\cal S}_{j}=\sqrt{-g}{T}^{0}_{j}$, and ${\cal B}^{j}=\sqrt{-g}(u^{0}b^{j}-u^{j}b^{0})$, respectively, the equations take on a traditional conservation law formulation

[TABLE]

for mass density,

[TABLE]

for energy,

[TABLE]

for momentum, and

[TABLE]

for magnetic induction. These fields $\{D,{\cal E},{\cal S}_{j},{\cal B}^{j}\}$ serve as the conserved fields for our relativistic MHD scheme, with $\{\rho,\rho\epsilon,V^{i},{\cal B}^{j}/\sqrt{-g}\}$ being the corresponding primitives.

2.2 Evolution

For this work we utilize the fully conservative High Resolution Shock Capturing (HRSC) Godunov scheme [22] to solve the Newtonian or relativistic MHD equations. To do so, we note that all of the equations have already been cast in conservation form

[TABLE]

where $\mathbf{U}(\mathbf{P})$, $\mathbf{F}^{i}(\mathbf{P})$, and $\mathbf{S}(\mathbf{P})$ are the arrays representing the conserved quantities, fluxes, and source terms, respectively. By integrating both sides with respect to volume and applying Gauss’ Theorem, we can rewrite the conservation equation in the form

[TABLE]

We can then discretize the equations using a finite volume representation as

[TABLE]

where $V$ is the cell volume and $A_{i}$ is the area of face $i$. The $\mathbf{F}^{i}$ are determined using the two-speed HLL approximate Riemann solver

[TABLE]

where $\mathbf{F}_{R}=\mathbf{F}(\mathbf{P}_{R})$ and $\mathbf{F}_{L}=\mathbf{F}(\mathbf{P}_{L})$ are the fluxes at the right- and left-hand side of each zone interface and $\mathbf{U}_{R}=\mathbf{U}_{R}(\mathbf{P}_{R})$ and $\mathbf{U}_{L}=\mathbf{U}_{L}(\mathbf{P}_{L})$ are the corresponding conserved quantities, each calculated from the extrapolations of the primitive variables to the zone interface, $\mathbf{P}_{R}$ and $\mathbf{P}_{L}$. Cosmos++ has options for piecewise-linear and piecewise-parabolic extrapolations, both limited by a one-dimensional MUSCL-type limiter with $1\leq\theta\leq 2$, and a multi-dimensional limiter designed to work on unstructured meshes by constructing upwind and downwind tetrahedral elements [17]. In this paper, all tests use the piecewise-parabolic option with the 1D MUSCL-type limiter and $\theta=1.5$. The bounding wave speeds are $c_{\mathrm{max}}\equiv\mathrm{max}(0,\leavevmode\nobreak\ c_{+,R},\leavevmode\nobreak\ c_{+,L})$ and $c_{\mathrm{min}}\equiv-\mathrm{min}(0,\leavevmode\nobreak\ c_{-,R},\leavevmode\nobreak\ c_{-,L})$, where $c_{\pm,R}$ and $c_{\pm,L}$ are the maximum right- and left-going waves speeds. Time integration is handled by a five-stage, strong-stability-preserving Runge-Kutta (SSPRK) method [23, 19].

Table 1 outlines the steps taken during a typical compute cycle. Usually the first physics package that is called is the primitive solver. In this step, a series of coupled equations are solved to extract primitive fields (mass density, internal energy, velocity, and primitive magnetic field) from evolved conserved fields (boost density, total energy, momentum, and conserved magnetic field or vector potential), after which the equation of state is applied to compute thermodynamic quantities, pressure, sound speed, and temperature. For Newtonian systems this procedure is straightforward, but relativity introduces a nonlinear inter-dependency of primitives so their extraction from conserved fields requires special numerical treatments. We have implemented several procedures for doing this, solving one, two, or five dimensional inversion schemes for MHD [24].

2.3 Refinement

2.3.1 Adaptive Mesh Refinement

Adaptive mesh refinement (AMR) in Cosmos++ uses a set of refinement criteria based on the local field, slope, or curvature of any user-specified code variable to actively adapt the mesh in each zone to best match the problem. A simple example of a field refinement criteria might be to refine a zone whenever the fluid density exceeds a given value and to de-refine if the density drops below a different value. A firm requirement in Cosmos++ is that neighboring zones never differ by more than one level of refinement. We have two options in the code for how to enforce this requirement. The first, referred to as “strict,” guarantees that any zone that meets a refinement criteria will indeed refine. In order to ensure that it will not then differ from any neighbor by more than one level of refinement, the refining zone queries the refinement level of all its neighbors and then forces them to refine as necessary to preserve the one level separation. The “fast” option takes the opposite approach. Whenever a zone is tagged for refinement, similar to “strict,” it polls the refinement level of its neighbors. If the refinement of the tagged zone will cause it to differ from any of its neighbors by more than one level of refinement (i.e. if any neighbors are currently at a lower level of refinement), then that zone is no longer tagged for refinement. Fig. 2.2 illustrates the practical difference. In this example, we start with a coarse $8\times 8$ base mesh and define a “circular” region of interest (red zones). These zones are then marked for refinement three separate times, for a total of four mesh levels. Note how the fast method refines “inward” from the edges of the targeted refinement region, whereas the strict method refines “outward.” Incidentally, this example also demonstrates why it is often useful to start a simulation fully refined and then allow it to de-refine as it can. Otherwise, the coarse initial resolution will imprint itself upon the problem as it has here.

De-refinement, on the other hand, solely follows the “fast” approach, i.e. zones only de-refine if after doing so they will not differ by more than one level from any direct neighbor. In this way, we err on the side of keeping greater resolution than we may need. The basic algorithmic steps for each option are presented in Table 2.

Whenever a zone is refined, the conserved fields are interpolated to the child zones in such a way as to maintain the same volume-integrated value as the parent. For de-refinement, the conserved fields in parent zones are set from the volume integral of all child zones. Thus, the conservative nature of our scheme is preserved for both refinement and de-refinement. New primitive fields are assigned to each zone after refinement/de-refinement is complete.

2.3.2 Moving Mesh Refinement

Mesh motion is implemented by a straight-forward replacement of the generic advective terms

[TABLE]

with

[TABLE]

Here $T^{0}_{\alpha}$ is used here to represent any evolved field.

The specification of the grid velocity $\mathbf{V}_{g}$ is a highly problem-dependent choice, and we do not in this paper advocate one method over another. Rather than developing and locking into specific algorithms, we have instead taken the tact of designing the code in such a way that users can easily implement their favorite options. However, whatever choices are made for mesh motion, they will almost certainly require rezoning strategies for mesh optimization in order to prevent tangling and highly distorted cells from forming and degrading the numerical solution. This is, in fact, the more crucial element needed for robust mesh motion. To this end, Cosmos++ has numerous options for mesh relaxation and optimization, including softening procedures generally adapted to Lagrangian flows, potential-based optimal equidistribution methods [25], velocity-based geometric conservation law methods [26], parabolic gradient flow relaxation methods [27], and variational equidistribution and mass optimization methods [28, 29].

Cosmos++ supports grid velocity specification at either the cell centers or nodes, although when specified at the centers, it must be projected onto the nodes in order to guarantee cell separation does not occur. This is done with simple mass weighting. After the node positions are physically relocated according to the nodal velocity and timestep, zone attributes (including face area vectors, cell volumes, and the spacetime metric) are updated for each cell [19, 20].

3 Divergence-Free Evolution of Magnetic Fields

Along with the evolution equations, (3)-(6) for Newtonian or (8)-(11) for general relativity, the time component of Maxwell’s equations provides a constraint condition: $\nabla\cdotp\mathbf{B}=0$ in Newtonian MHD or $\partial_{j}\mathcal{B}^{j}=0$ in general relativity. In both cases, the meaning is the same: there can never be more magnetic field entering a region of space than leaving it, or that the divergence of the magnetic field through a closed surface must always be zero. However, most numerical implementations of the induction equation [(6) or (11)] will not be guaranteed to preserve the $\nabla\cdotp\mathbf{B}=0$ constraint, even if the initial field is divergence free.

The original version of Cosmos++ [17] included a simple parabolic divergence cleanser option [9], intended to dampen any accumulated divergence error. This option was extended in later work [e.g. 30, 31] to include elliptic and hyperbolic elements to filter and propagate away the divergence error. Such methods are generally successful at maintaining $\nabla\cdot\mathbf{B}=0$ to about the level of truncation errors. Once the fully conservative HRSC MHD method was implemented in Cosmos++ [32], we introduced a staggered representation of the magnetic field, with the primary (conserved) fields ($B^{i}$ for Newtonian or ${\cal B}^{i}$ for GR) being treated as area-weighted variables located at face centers (see Fig. 3.1) and implemented a form of CT. The new VP method shares many elements with our staggered CT scheme.

3.1 Vector Potential Evolution

Instead of evolving the magnetic field directly using the induction equation [(6) or (11)], another option is to evolve the magnetic vector potential [33, 15]. In doing so, we replace the evolution equations (6) and (11) with the following

[TABLE]

where $\mathscr{E}_{i}=-\epsilon_{ijk}V^{j}\mathcal{B}^{k}$, is the $i^{\mathrm{th}}$-component of the electric field (or emf), centered on the appropriate cell edge (Fig. 3.1) and $\Phi$ is the electrostatic potential.

The most robust method we have found so far for constructing the edge-centered electric fields is to simply average the surrounding face-centered fluxes recovered from the Riemann solver (being careful of signs), as described in [3]. For example,

[TABLE]

where $F_{B^{1}}^{2}=B^{1}V^{2}-B^{2}V^{1}$ and $F_{B^{2}}^{1}=B^{2}V^{1}-B^{1}V^{2}$ are the properly upwinded Godunov fluxes resulting from equation (15). Similar procedures are used for the other components of $\mathscr{E}_{i}$. Note here that we always use the fluid velocity, $V^{j}$, in defining the electric field and not the relative velocity with the grid motion subtracted out, $(V^{j}-V_{g}^{j})$. This reconstruction is second-order accurate in space and retains the temporal accuracy of whatever time integration is used [3]. This is consistent with the order of all source terms in Eqs. (4), (5), (9), and (10). We will explore higher-order interpolations once we port this method over to our discontinuous Galerkin (DG) version of Cosmos [19]. Another option would be to implement a multi-step or multi-dimensional Riemann solver [33, 12, 13].

The electrostatic potential in Eq. (18)-(19) is dependent on a choice of gauge. For our work, we have adopted the simplest, algebraic EM gauge condition, in which

[TABLE]

and the evolution simplifies to

[TABLE]

Another choice would be the Lorentz gauge [34], which requires an additional evolution equation for $\Phi$:

[TABLE]

This could easily be added to Cosmos++, but has, so far, proven unnecessary.

In Cosmos++, we choose to have each zone store its own bottom-left edge-centered vector potential components in two dimensions or its back-bottom-left edge-centered vector potential components in three dimensions, as illustrated in Fig. 3.2. This allows each zone to only store a single set of edge-centered vector potentials, while still facilitating easy reconstruction of the magnetic field. This is done during the primitive inversion step by taking the curl of the vector potential:

[TABLE]

or in component form, explicitly, as

[TABLE]

Since the divergence of a curl is always zero,

[TABLE]

this procedure is mathematically guaranteed to maintain a divergence-free field. In practice, it preserves the following finite volume representation of divergence:

[TABLE]

where again $V$ is the cell volume and $A_{i}$ is the area of face $i$. When normalized by $\sqrt{P_{B}}/\Delta x$, where $\Delta x$ is a characteristic zone length, this is preserved to round-off error in each zone. the same representation of $\nabla\cdotp\mathbf{B}=0$ as the CT scheme (30). To see this explicitly, we can expand Eq. (26)-(29). First, the magnetic field components at the faces are

[TABLE]

The magnetic field divergence through the surface of the cell at position ($i$, $j$, $k$) is then

[TABLE]

It is easy to see that all of the terms in Eq. (34) cancel out pairwise.

The staggered magnetic field recovered in this way is used for two purposes. The first is to overwrite the appropriate component of the extrapolated primitive field used in the flux reconstruction [Eq. (15)] at each face. The other is to calculate cell-centered values of the primitive magnetic fields as required, for instance, for calculating the magnetic pressure and recovering the internal energy of the gas. In this case we use the volume-averaged fields

[TABLE]

We store both the zone-centered and staggered magnetic fields along with the edge-centered vector potential components.

Occasionally, we may want to solve problems containing symmetries that may allow us to reduce the dimensionality. In two dimensions, the $x_{3}$-component of the magnetic field is divergence free by definition, as all gradients in that direction are assumed to be zero. Looking at the expressions for $B^{1}$ and $B^{2}$ in Eq. (31) and (32), it is clear that there is no reason to evolve $\mathcal{A}^{1}$ nor $\mathcal{A}^{2}$ in two dimensions, as they are not needed. Thus, when considering two-dimensional problems, we only evolve $\mathcal{A}^{3}$ and use it to recover $B^{1}$ and $B^{2}$. $B^{3}$ is updated through the standard induction equation (6). In one dimension, there is no need to evolve any of the vector potential components, as there is no chance of building up divergence in this case. Here, the induction equation (6) is entirely sufficient.

Before moving on, we mention one slight drawback of the VP method: it can make the initialization of problems slightly more difficult. This is because, in cases where one has a starting magnetic field configuration in mind, it may not be an easy matter to find the corresponding vector potential, as there is no straightforward way to invert Eq. (24). For simple magnetic field configurations, as in all the test problems in this paper, this is not a major hurdle, but for complex configurations, this could be challenging.

3.1.1 Mesh Refinement Considerations

For a uniform mesh, the VP method is functionally equivalent to a conventional staggered CT scheme. Effectively, it just changes the order of the spatial and temporal operations in updating the magnetic field. However, because the vector potential itself is not subject to any constraint, we can refine and de-refine the mesh and evolved fields without introducing magnetic monopoles. This is the principle advantage of the VP method. Still, care must be taken in updating the vector potential and magnetic field components near refinement boundaries, where neighbors may be at different refinement levels. We discuss these considerations in this section.

Let us look first at the VP update. The source term for the VP update equation [(18) or (19)] is just the edge-centered electric field. We see from Eq. (20) that this field comes from an average of the face-centered magnetic field fluxes on the four surrounding faces. However, near refinement boundaries, a given edge could have as many as eight (four) surrounding faces in three (two) dimensions or as few as three, with the faces lying at potentially different distances from the edge center. For this reason, we must replace Eq. (20) with one that weights the fluxes appropriately. As a basic rule, we always use the closest set of fluxes in each direction. This can sometimes mean that a given zone does not even use its own fluxes when calculating one (or more) of its edge-centered vector potentials. Fig. 3.3 shows a two-dimensional example where the zone in question only uses fluxes from its more refined neighbors. The corresponding expression for the edge-centered electric field is

[TABLE]

where $w_{T}$ is the sum of all the individual weights

[TABLE]

and $\mathbf{x}$ is the appropriate position vector. For a uniform mesh, the result is identical to Eq. (20).

As mentioned, there is a possibility for “hanging” edges to have fewer than four surrounding faces. Fig. 3.4 illustrates such a case. Here, we simply replace the edge-vector potential in question with the average of the already calculated edge-vector potentials along the neighboring edges

[TABLE]

The physical motivation for doing so is that this procedure guarantees that the sum of the magnetic fluxes, $\Phi_{B}$, through the refined faces equals the magnetic flux through the neighbor’s unrefined face as it should, i.e. $(\Phi_{B})_{0}=(\Phi_{B})_{1}+(\Phi_{B})_{2}$, where

[TABLE]

Following the same principles, we can extend these procedures for updating the vector potential components at refinement boundaries to three dimensions.

The next step is to obtain the face-centered magnetic fields from the edge vector potentials by computing their curl. This is accomplished computationally by taking gradients of the edge vector potentials surrounding each face as in Eq. (31)-(33). Again, though, there are special considerations at refinement boundaries. For example, in the case of hanging edges, there may be no corresponding vector potential at that location. For these cases, we follow the same procedure as for calculating vector potentials for edges that have fewer than four surrounding faces, namely we average the already calculated edge-vector potentials from the neighboring edges as in Eq. (40) before taking any gradients.

3.1.2 Mesh Motion Considerations

As currently implemented, the new VP method requires that magnetic field components lie normal to their respective faces as in Fig. 3.1. For this reason, we limit our considerations in this paper to conformal motion, such as uniform translation of the grid, which preserves the orthogonality, though conformal expansion and contraction could be treated as well. However, since Cosmos++ supports full mesh motion, including mesh distortion, as described in Sec. 2.3.2, we plan to expand the mesh motion capabilities of the VP method in the future. One way to go about that would be to perform the magnetic field updates in an orthogonal reference zone. To construct the reference zone, we could use the properties of the distorted zone (edge lengths and inner products of face areas) to construct a local metric that encapsulates the distortion. We would then use that metric to transform the vector potential and magnetic field components into the orthogonal reference frame, where the magnetic fields would be normal to the zone faces, allowing the VP method to be applied. This would result in a magnetic field that is divergence free in the reference frame. We already follow a similar procedure for solving the hydrodynamic equations in CosmosDG [19], so this would be a natural thing to do when we port the VP method over to that version of the code.

4 Code Tests

In this section, we verify that our new MHD method is able to handle a large suite of standard code tests. We only consider multi-dimensional tests, as the VP method plays no role in one dimension, and one dimensional MHD tests using Cosmos++ were already presented in [17]. We pay particular attention to any accumulation of $\nabla\cdotp\mathbf{B}$ and report the $L$-2 norm error of this quantity

[TABLE]

for each test. We also compare to known analytic and published solutions whenever possible to assess how our method compares and calculate its convergence order. All tests are done with the time integration set to third order. This ensures that our errors will reveal the true spatial accuracy of the method.

4.1 Code Tests without Mesh Refinement

Before considering test problems that involve refinement (adaptive and moving mesh), we consider a few test problems on uniform, static meshes, just to confirm things are working correctly.

4.1.1 Circularly Polarized Alfvén Wave

Since a circularly polarized Alfvén wave in a periodic box has a known analytic solution [7], it provides an excellent starting test and way to confirm the accuracy and convergence order of our code. To engage all components of the VP method, we make this a three-dimensional test following the approach of [35] by initializing the problem in the rotated coordinate system, $(x_{1},x_{2},x_{3})$, where the orientation with respect to the standard $(x,y,z)$ coordinates is given by

[TABLE]

where $\sin\alpha=2/3$ and $\sin\beta=2/\sqrt{5}$. The mesh extends over the range $0\leq x\leq 3$, $0\leq y\leq 1.5$, and $0\leq z\leq 1.5$, resolved by a $2N\times N\times N$ grid, and is periodic in all directions (except the vector potential $\vec{\mathcal{A}}$ is not periodic). The initial solution has mass density $\rho=1$, gas pressure $P=0.1$, adiabatic index $\Gamma=5/3$, and magnetic field $\mathbf{B}=\{1,0.1\sin(kx_{1}),0.1\cos(kx_{1})\}$ in $(x_{1},x_{2},x_{3})$ coordinates, with wave vector $k=2\pi/\lambda$ and wavelength $\lambda=1$. Consistent with our new approach, we actually initialize the magnetic vector potential $\vec{\mathcal{A}}=\{0,[0.1\sin(kx_{1})]/k-0.5x_{3},[0.1\cos(kx_{1})]/k+0.5x_{2}\}$, which ensures $\nabla\cdot\mathbf{B}=0$ initially. However, because of truncation errors associated with recovering the magnetic field from the vector potential, the initial conditions are subject to grid-level magnetic pressure perturbations, which drive secondary compressive waves.

The wave can either be set up to propagate parallel to the $x_{1}$-axis, in which case the corresponding fluid velocity component is $v_{1}=0$ (traveling case), or remain stationary, in which case the fluid velocity component is $v_{1}=1$ (standing case). The other velocity components are $v_{2}=0.1\sin(2\pi x_{1}/\lambda)$ and $v_{3}=0.1\cos(2\pi x_{1}/\lambda)$. Since $v_{2}$, $v_{3}$, and the boundaries are periodic, then for these parameter choices, the traveling wave should return to its initial state after a period of $t=1$. This makes it easy to test the convergence of our method by comparing the final state against the initial one. We first calculate the $L$-1 norm errors of each of our conserved fields, $U_{s}$, as

[TABLE]

where $(U_{s})^{n}_{i,j,k}$ is the final numerical solution and $(U_{s})^{0}_{i,j,k}$ is the initial (analytical) solution for field $s$. Most of these errors are reported in Table 3. In Fig. 4.1, we present a plot of the root mean square of the $L$-1 norm errors

[TABLE]

as a function of the grid resolution, $N$. Our errors compare favorably to those reported [35] and show almost exactly second order convergence, as expected. Note, though, that for the standing wave case, we could only run the test to a time of $t=0.25$ (a quarter of a traveling wave period). For longer run times, the solution became unstable, possibly due to a parametric instability [36, 37], especially at higher resolutions. The $L$-2 norm error, $||\nabla\cdotp\mathbf{B}||_{2}$, never exceeds $1.05\times 10^{-17}$ for the $N=8$ traveling or standing test.

4.1.2 MHD Vortex

The MHD vortex problem [38] is another test that can be used for convergence studies. In fact, as it is an exact solution of the MHD equations in two dimensions, it perhaps gives a better measure of the true error than the circularly polarized Alfvén wave, which is really a one-dimensional problem. The vortex is set up in an $N\times N$ computational domain covering $-5\leq x\leq 5$, $-5\leq y\leq 5$, with periodic boundaries in both directions. The unperturbed background has properties $\{\rho,P,v^{x},v^{y},B^{x},B^{y}\}=\{1,1,1,1,0,0\}$, with $\Gamma=5/3$. The vortex is initialized via the following perturbations

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$. The perturbations for $\{\delta B^{x},\delta B^{y}\}$ have the exact same form, though for the VP method, we initialize the perturbation

[TABLE]

The centrifugal force of the vortex motion, the centripetal force of the magnetic tension, and the magnetic pressure all contribute to the dynamical force balance, along with gas pressure. To maintain the balance, the pressure must have a perturbation of the form

[TABLE]

The problem is run until $t=10$, by which time the vortex returns to its initial location.

We again test convergence by comparing the final state against the initial one, with $L$-1 errors reported in Table 4. Our errors are about an order of magnitude lower than those reported [38] and show almost exactly second-order convergence. Furthermore, the $L$-2 norm error, $||\nabla\cdotp\mathbf{B}||_{2}$, never exceeds $6.98\times 10^{-18}$ for the $N=50$ test.

4.1.3 Current Sheet

The current sheet test [39] is another useful multi-dimensional MHD test. It particularly challenges the robustness of codes in highly magnetized regions. In it, oppositely directed magnetic fields abut one another at multiple current sheets. An oscillating sheer flow is then set up, which forces some level of reconnection within the current sheets. Some MHD schemes go unstable and ultimately fail this test problem.

This relativistic test is set up on a $200\times 200$ grid with $0\leq{x}\leq 2$, $0\leq{y}\leq 2$, uniform mass density $\rho=1$ and gas pressure $P=\beta/2$. We test values in the range $10^{-3}\leq\beta\leq 10^{-1}$. The velocity is initialized as $V^{x}=V_{0}\cos(\pi y)$ and $V^{y}=V^{z}=0$ with amplitude $V_{0}=0.2$. The magnetic field is initialized as $\mathcal{B}^{x}=\mathcal{B}^{z}=0$, $\mathcal{A}_{x}=\mathcal{A}_{y}=0$, $\mathcal{B}^{y}=-1$ and $\mathcal{A}_{z}=\mathcal{B}^{y}x$ for $0.5\leq{x}\leq 1.5$, $\mathcal{B}^{y}=1$ and $\mathcal{A}_{z}=\mathcal{B}^{y}(x-1)$ for $x<0.5$, and $\mathcal{B}^{y}=1$ and $\mathcal{A}_{z}=\mathcal{B}^{y}(x-3)$ for $x>1.5$.

The evolution of this test for $\beta=0.01$ is shown in Fig. 4.2. Along the current sheets, magnetic islands begin to form as grid-scale reconnection occurs. As the simulation proceeds, smaller magnetic islands merge into larger structures until a stable configuration of four islands is reached. At no point does the $L$-2 norm error, $||\nabla\cdotp\mathbf{B}||_{2}$, exceed $5.26\times 10^{-19}$ for this test. Only for very strong magnetic fields ($\beta<10^{-3}$ for $V_{0}=0.2$) or high amplitudes ($V_{0}>0.5$ for $\beta<10^{-1}$) are no stable solutions reached. This range of parameter space is very similar to what is achieved with the HLLE Riemann solver in Athena [40].

4.1.4 Field-loop Advection

The field-loop advection test [40] is another multi-dimensional MHD test that has been used to demonstrate failure modes of some MHD codes [41]. In this case, some schemes designed to preserve a divergence-free field yield unphysical results owing to an incorrect treatment of the Lorentz force. The failure manifests itself as an anomalous build up of fluid velocity (acceleration) in the direction normal to the field loop propagation.

We initialize the relativistic version of this problem [40], setting the mass density $\rho=1$ and the gas pressure $P=3$. The problem is conducted on a $2N\times{N}$ grid with $N=128$, $-1\leq x\leq 1$, and $-0.5\leq y\leq 0.5$. The velocity in the $x$-direction is set to $V^{x}=0.2/\sqrt{6}$, while the velocity in the $y$-direction is set to $V^{y}=0.1/\sqrt{6}$. The velocity in the $z$-direction is also set to $V^{z}=0.1/\sqrt{6}$ and should remain constant. The magnetic field is initialized as

[TABLE]

with $\mathcal{A}_{0}=10^{-3}$, $R=0.3$, and $r=\sqrt{x^{2}+y^{2}}$.

The field loop propagates in a periodic box, so should return to its starting location undisturbed. Fig. 4.3 displays the magnetic field lines after zero, one, and two grid crossings. The lines are slightly distorted as they propagate across the grid, but the general shape and structure of the field loop is preserved. Further, the total magnetic energy on the grid drops by $<13$%. More importantly, we see no significant increase in $V^{z}$; after two grid crossings, its fractional error is $\delta V^{z}/V^{z}\leq 10^{-7}$, and the $L$-2 norm error of divergence, $||\nabla\cdotp\mathbf{B}||_{2}$, remains below $6.40\times 10^{-16}$.

4.1.5 Orszag-Tang

The Orszag-Tang problem [42] is another MHD vortex problem that has been adopted as one of the standard test problems for MHD. Unlike the smooth MHD vortex above, the Orszag-Tang problem produces turbulence and shocks. To initialize it, we set the mass density to $\rho=\Gamma^{2}/4\pi$, the gas pressure to $P=\Gamma/4\pi$, and the adiabatic index to $\Gamma=5/3$. The initial velocity and magnetic fields are functions of the position, which creates the vortex: $\mathbf{v}=\{-\sin(2\pi y\},\sin(2\pi x))$, and $\mathbf{B}=\{-\sin(2\pi y),\sin(4\pi x)\}$. The evolution of the mass density and magnetic field using the VP scheme at $512\times 512$ resolution is shown in Fig. 4.4. These results are consistent with previously published solutions [e.g. 5, 6, 43, 7, 10] and have very low divergence error, with the $L$-2 norm error, $||\nabla\cdotp\mathbf{B}||_{2}$, never exceeding $3.11\times 10^{-20}$. Notably, the solution remains symmetric, which is significant given the face-centering of our conserved magnetic field components and edge centering of our vector potential components.

4.1.6 Rotor

The next test looks at torsional Alfvén waves given off by a rotating cylinder (the rotor) of fluid threaded by a magnetic field. The magnetic field transfers angular momentum from the rotor to the background fluid. As the field wraps around the rotor, it also distorts its shape.

This problem is set up on a two-dimensional unit square, resolved with $400\times 400$ zones. A static ($\mathbf{v}=\mathbf{0}$) background is set with $\rho=1$, $P=1$, and $\Gamma=1.4$. The rotor is initialized in the center of the grid with a density of $\rho=10$ inside $r<0.1$, tapering off to the background value over a distance of $\Delta r=0.015$. Similarly, the rotor has a uniform angular frequency of 20 inside $r<0.1$, tapering off to zero over the same $\Delta r=0.015$. A uniform magnetic field in the $x$-direction of magnitude $B^{x}=5/\sqrt{4\pi}$ threads the entire domain. This is initialized from a vector potential of the form ${\mathcal{A}}_{z}=B^{x}y$. The problem is run to a time $t=0.15$. Figure 4.5 shows contours of the final density, gas pressure, Mach number, and magnetic pressure. The results look quite comparable to previously published results [3, 7], and again the $L$-2 norm error of divergence, $||\nabla\cdotp\mathbf{B}||_{2}$, remains low ($<1.68\times 10^{-19}$).

4.1.7 Blast

Our final test problem is the explosion of a centrally over pressurized region into a low pressure, low $\beta$ ambient medium. Like the Orszag-Tang problem, this test combines shocked flows and strong magnetic fields. Like the rotor problem, there is also a smooth flow region. Similar to those tests, the results are not particularly quantitative in their measure of the accuracy, but this test adds to our confidence in the robustness of the method. Variants on this problem have been presented by a number of authors [e.g., 3, 43, 41].

The computational domain extends over $-0.5\leq x\leq 0.5$ and $-0.5\leq y\leq 0.5$. A uniform, static ($\mathbf{v}=\mathbf{0}$) background of $\rho=1$, $P=0.1$, and $\Gamma=1.4$ is set everywhere. A magnetic field of $\{B^{x},B^{y}\}=\{100/\sqrt{4\pi},0\}$ threads the entire domain. As in the previous problem, this is initialized with a vector potential of the form ${\mathcal{A}}_{z}=B^{x}y$. This gives a background $\beta$ of $2.5\times 10^{-4}$. Within a sphere of radius $r=0.1$, centered at the origin, we set the pressure to 1000 to initiate the blast. The problem is run to a stop time of $t=0.01$. Figure 4.6 shows contours of the final density, gas pressure, velocity magnitude, and magnetic pressure. The results look quite comparable to previously published results [3, 35], and once again the $L$-2 norm error of divergence, $||\nabla\cdotp\mathbf{B}||_{2}$, remains low ($<3.10\times 10^{-19}$).

4.2 Code Tests with Mesh Refinement

Now that we have demonstrated that our VP method performs well on uniform, static meshes, we move on to consider problems involving refined (static and active) and moving (uniform translation) meshes.

4.2.1 MHD Vortex Revisited

First, we repeat the MHD vortex problem with mesh motion. We set the problem up exactly as in Section 4.1.2, except we add a uniform grid velocity equal to the propagation speed of the vortex, $\{V^{x}_{g},V^{y}_{g}\}=\{1,1\}$. This causes the vortex to remain stationary within the simulation domain, allowing us to quantitatively confirm that our convergence order is preserved whenever the grid velocity in non-zero. This is confirmed in Table 5, where we report the $L$-1 errors for this test. These are comparable to the errors reported in Table 4 for the static mesh, and more importantly, show the same second-order convergence. Also important, the $L$-2 norm error of divergence, $||\nabla\cdotp\mathbf{B}||_{2}$, is still below roundoff ($3.41\times 10^{-18}$ for the $N=50$ test).

4.2.2 Field-loop Advection Revisited

Next we redo the field-loop advection test using independently both adaptive and moving mesh refinement. In Fig. 4.7, we reproduce our results from Fig. 4.3, but instead of using a single-level $256\times 128$ mesh, we use a $64\times 32$ base mesh with two levels of adaptive refinement. Hence, both simulations have the same peak resolution. We begin with the entire grid at the highest level of refinement to ensure proper initialization, but then allow the grid to de-refine and refine as the simulation evolves and the field loop moves. The refinement variable in this case is the magnetic pressure, with refinement set whenever $P_{B}>10^{-9}$ and de-refinement whenever $P_{B}<10^{-10}$. We use the strict refinement setting (Sec. 2.3) and test for refinement every tenth cycle. The point of this test is to demonstrate that we can achieve similar results on an adaptive, multi-layer mesh as we do on a single, static mesh of equivalent peak resolution. We note specifically that the drop in total magnetic energy is now 17% and the error in $V^{z}$ is comparable with and without refinement.

Even more so than adaptive mesh refinement, the field loop advection test lends itself to the application of moving mesh. We, therefore, set up another case where we introduce a grid velocity equal to the fixed-frame advective velocity of the field loop, i.e. $V^{x}_{g}=0.2/\sqrt{6}$ and $V^{y}_{g}=0.1/\sqrt{6}$. In this way, we keep the loop centered on the moving grid. We, again, let the simulation evolve for what would be two crossing periods. We find that the results are quite comparable to the fixed mesh results in Sec. 4.1.4 with the same level of magnetic energy decay, $V^{z}$ errors, and $L$-2 norm error of divergence.

4.2.3 Orszag-Tang Revisited

The next test we try with mesh refinement is the Orszag-Tang problem from Sec. 4.1.5. Here we set up static refinement over two quadrants of the grid, leaving the other two quadrants at a lower resolution ($128\times 128$ base resolution with two levels of refinement added in the lower left and upper right quadrants). The goal is to confirm that no artificial features or mesh imprinting are seen whenever a very dynamical, magnetized fluid moves back and forth across refinement boundaries. Fig. 4.8a shows that the final solution shows the same characteristic features as our single layer mesh case (Fig. LABEL:figure:ot) though not as sharp in the unrefined quadrants. Crucially, we find no evidence of mesh imprinting, while maintaining a divergence error very close to machine precision ($||\nabla\cdotp\mathbf{B}||_{2}<1.10\times 10^{-16}$; Fig. 4.8b).

4.2.4 Magnetized Black Hole Torus

Next we consider the astrophysically more interesting problem of a magnetized torus of gas orbiting a rotating black hole. Although there is no known analytic solution for this problem, it is sufficiently well documented in the literature [e.g. 44, 45, 46, 47] to serve as a useful test of the code. It also represents one class of problems for which the new VP method is intended to be used. The particular equilibrium configuration we start from is the constant specific angular momentum torus solution presented in [48]. Specifically, the spin of the black hole is $a/M=0.9375$, the inner radius of the torus is set to $r_{\rm in}=6GM/c^{2}$, the density maximum is located at $r_{\rm max}=12GM/c^{2}$, and we use an ideal gas EOS with an adiabatic index of $\Gamma=4/3$. To the hydrodynamic solution we add an initially weak poloidal magnetic field loop defined by the vector potential

[TABLE]

which produces poloidal field loops coincident with the isodensity contours of the torus. The field strength is set such that $P_{\rm max}/P_{B,{\rm max}}=100$, where the global maxima of the gas $P_{\rm max}$ and magnetic $P_{B,{\rm max}}$ pressures do not necessarily coincide. In order to spur growth of the magneto-rotational instability inside the torus, the gas pressure is perturbed as

[TABLE]

where $X_{p}$ is a uniformly distributed random variable between $-0.02$ and $0.02$.

In the background region not specified by the torus solution, we set a floor model with $\rho_{\mathrm{fl}}=10^{-5}r^{-3/2}$, $e_{\mathrm{fl}}=10^{-7}r^{-5/2}$, and $V^{i}=0$. Obviously once the evolution begins, the background is no longer static. At first it falls toward the hole, but ultimately most of the background region is filled with a magnetized “corona” and a low-density, high-magnetization outflow launched by pressure forces and the opening of magnetic field lines.

Part of our motivation for running this particular test is that we already have results available from the CT version of Cosmos++ with which to compare. As mentioned throughout this paper, one of the advantages of the VP method over CT is that we are able to utilize mesh refinement. For black hole accretion problems done on spherical polar grids, one of the difficulties one faces is the very small zone sizes near the polar axis. On a single layer mesh, the only way to avoid this problem are to cut out a small region around the pole, as in the left panel of Fig. 4.9. However, with (static) refinement, one can extend the base level grid all the way to the pole without suffering prohibitive timestep constraints because the zones do not necessarily need to be small, while the desired level of refinement can be added away from the pole, as in the right panel of Fig. 4.9. Here, both simulations are run with comparable peak resolutions, yet the VP one runs with almost an order of magnitude larger timestep, while also including the full range of $\theta$.

Both simulations use a logarithmic radial coordinate of the form $x_{1}=1+\ln(r/r_{\mathrm{BH}})$ with radial boundaries at $r_{\mathrm{min}}=0.98r_{\mathrm{BH}}$ and $r_{\mathrm{max}}=354GM/c^{2}$, where $r_{\mathrm{BH}}$ is the radius of the black hole event horizon, and a concentrated latitude coordinate of the form $\theta=x_{2}+0.25\sin(2x_{2})$ with the $x_{1}$ and $x_{2}$ coordinates being uniformly distributed. For the single-layer ($96\times 96$) mesh (CT) simulation, the grid covers the angular range $0.02\pi\leq\theta\leq 0.98\pi$, while the base-level ($24\times 24$) mesh for VP covers the full range $0\leq\theta\leq\pi$. Two levels of refinement are added on top of this base mesh in the region $r_{\mathrm{min}}\leq r\leq 50GM/c^{2}$ and $0.02\pi\leq\theta\leq 0.98\pi$. Both simulations cover the full range of $0\leq\phi\leq 2\pi$, with periodic boundary conditions at $\phi=\{0,2\pi\}$. Outflow boundary conditions are applied on the inner and outer radial boundaries, as well as the cutout region in the CT simulation. Special reflecting boundaries are used along the pole in the VP case.

Fig. 4.10 gives a more quantitative comparison between the two simulations. It shows time-averaged radial profiles of various disk quantities for both the CT and VP simulations. Spatial averages are calculated as

[TABLE]

which are then time averaged over approximately 100 dumps in the interval $5000\leq t\leq 10000GM/c^{3}$. Both simulations show nearly identical profiles, which also compare favorably to the collection of results published in [47], confirming that the new VP method works quite well for this target problem.

5 Discussion & Conclusions

In this paper, we have reported our implementation of the vector potential method for divergence-free evolution of magnetic fields in both Newtonian and relativistic MHD. The key to this approach is to evolve the magnetic vector potential ($\vec{\mathcal{A}}$) directly, rather than the magnetic field ($\mathbf{B}$). The magnetic field components can then be recovered as needed from $\mathbf{B}=\nabla\times\vec{\mathcal{A}}$, a procedure that is guaranteed to satisfy $\nabla\cdot\mathbf{B}=0$ to machine precision. The advantage in our case is that this approach makes it much easier to handle adaptive and moving mesh refinement (and hopefully adaptive order refinement in the future), while still preserving a divergence-free field.

We have shown through a series of multi-dimensional test problems that the method performs as anticipated. These tests were chosen to facilitate comparison with previously published results and to provide quantitative measures of the code’s accuracy and convergence. We found the code maintained stable, accurate solutions in a number of challenging setups, performing as well or better than other published codes on fixed meshes. Further, the method showed no sign of any of the worrisome pathologies that were feared for vector potential evolution [4] nor the critical failures on select tests identified in other MHD approaches [7, 41, 35, 40]. Most importantly, the method was able to do all of this on adaptively refined and uniformly translating meshes, with no signs of mesh imprinting or other issues.

Acknowledgements

We thank the anonymous referees for their helpful suggestions in improving this manuscript. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. P.C.F. acknowledges support from National Science Foundation grant AST-1616185. Work by P.A. was performed in part under the auspices of the U. S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. P.L.S. acknowledges support from the College of Charleston Undergraduate Research and Creative Activities Board, through SURF grant SU2018-24.