© ESO 2019

1. Introduction

The prospect of horizon-scale images of the two nearest supermassive black hole (SMBH) candidates Sgr A* and M 87, soon to be obtained by very long baseline interferometry (VLBI) arrays (Doeleman et al. 2008; Akiyama et al. 2015; Fish et al. 2016; Goddi et al. 2017; Issaoun et al. 2019; Event Horizon Telescope Collaboration 2019a,b,c,d,e,f), opens the possibility to study in both fundamental and astrophysical aspects of these objects in great detail. The most exciting possibilities include direct observations of the black hole shadow (Abdujabbarov et al. 2015; Younsi et al. 2016; Psaltis et al. 2016, 2015; Mizuno et al. 2018) and the movement of hot spots in the accretion flow (Broderick & Loeb 2006), as well as deciphering the cause of flares and non-thermal emission mechanisms (Özel et al. 2000; Dexter et al. 2012; Davelaar et al. 2018, 2019). In addition to Sgr A* and M 87, VLBI observations will provide high-resolution images of other sources of great interest. For instance, observations of the two-sided jet in the nearby active galactic nucleus (AGN) Cen A (Kim et al. 2018) offer unique possibilities to study the dynamics of jet formation and propagation in these objects. VLBI images could also be crucial to discriminate between models for the periodic variability in the BL Lac object OJ287, namely a secondary SMBH plowing through the accretion disk (Valtonen et al. 2008), or a precessing disk (Katz 1997) or jet (Britzen et al. 2018).

General relativistic magnetohydrodynamical (GRMHD) simulations are an invaluable tool and possibly the only available one to assess the validity of theoretical models with respect to astronomical observations. However, an important challenge for codes performing such simulations is the interplay between different scales relevant for the accretion process. Such a large excursion of lengthscales can easily make some problems prohibitive even for massively parallel numerical simulations. In fact, to aid with the interpretation of astronomical observations, simulations of accretion flows must reproduce global quantities of the system, such as accretion rates, spectra and light curves, or large-scale features such as the position of re-collimation shocks in the jet. However, many of these observables are determined by turbulent phenomena occurring at much smaller scales, which must be resolved to a reasonable degree to properly reproduce the physics. For instance, the mechanism of angular momentum transport that permits accretion (Shakura & Sunyaev 1973) is now believed to be magnetic turbulence driven by the magnetorotational instability (MRI, Balbus & Hawley 1991), and the processes re-collimating the jet and leading to equipartition (Porth & Komissarov 2015), are driven by several magnetohydrodynamical instabilities (see e.g., Mizuno et al. 2012) which interact with the acceleration processes (see e.g., Aloy & Rezzolla 2006).

By concentrating resolution only at places where it is most needed, adaptive mesh refinement (AMR) offers a flexible solution to this problem. In addition, a great advantage of AMR is the possibility of the grid to follow moving features while applying sufficient resolution. This can be of special importance, for instance, for systems with complex geometries, such as precessing jets (Britzen et al. 2018; Liska et al. 2018) and disks (Fragile et al. 2007; McKinney et al. 2013) or tidal disruption events (Tchekhovskoy et al. 2014).

Among the variety of GRMHD codes reported on the literature (Hawley et al. 1984; Kudoh 2000; De Villiers & Hawley 2003; Gammie et al. 2003; Baiotti et al. 2005; Duez et al. 2005; Anninos et al. 2005, 2017; Antón et al. 2006; Mizuno et al. 2006; Del Zanna et al. 2007; Giacomazzo & Rezzolla 2007; Radice & Rezzolla 2012; Radice et al. 2014; McKinney et al. 2014; Etienne et al. 2015; Zanotti & Dumbser 2015; White et al. 2016; Meliani et al. 2017; Liska et al. 2018; Fambri et al. 2018) an increasing number is making use of AMR techniques, (see e.g., Anninos et al. 2005; Zanotti et al. 2015; White et al. 2016; Liska et al. 2018). One of them is the publicly available BHAC (the Black Hole Accretion Code¹), which was described in Porth et al. (2017) and is the main focus of this work. The consistency of results obtained by several of the codes in this list is thoroughly analyzed in the forthcoming work of Porth et al. (2019).

An important challenge to the application of AMR in GRMHD simulations comes from the solenoidal constraint of the magnetic field, ∇ ⋅ B = 0, which must be fulfilled in order for the simulation to represent a physical state. The numerous schemes devised to enforce this constraint present different degrees of complication to be coupled with AMR techniques. Divergence-cleaning schemes, which require only the solution of an additional equation damping and propagating away violations to the constraint, can be coupled straightforwardly to AMR grids, as done by Anninos et al. (2005), Zanotti et al. (2015), Porth et al. (2017).

Unfortunately, comparisons among several schemes in non-relativistic magnetohydrodynamics (Toth 2000; Balsara & Kim 2004; Mocz et al. 2016) show that all tested variants of the divergence-cleaning method produce important spurious oscillations in the magnetic-field energy. Balsara & Kim (2004) have attributed these oscillations to the non-locality introduced by the constraint-damping equation, thus suggesting that all divergence-cleaning methods may suffer from the same problem. In a previous work (Porth et al. 2017), we have observed the same behavior in GRMHD when comparing simulations on uniform grids performed using the divergence-cleaning technique known as Generalized Lagrange Multipliers (GLM, Dedner et al. 2002) and the method known as flux-interpolated Constrained Transport (flux-CT, Toth 2000).

By adopting a discretization of Faraday’s law consistent with that of the constraint, Constrained-Transport (CT) methods maintain the solenoidal condition to machine precision throughout the solution of the GRMHD equations. In general, this requires to define the components of the magnetic field at cell faces, in a grid that is staggered with respect to that of the hydrodynamic variables.

A drawback of these methods is that the fulfillment of the constraint at every next step depends on its fulfillment at the current step; therefore, the initial condition must be divergence-free and care must be taken in order not to generate magnetic-field divergence by any other means. This can happen, in particular, at boundary cells, at coarse-fine interfaces in AMR simulations, and when creating and destroying blocks at prolongation (refining) and restriction (coarsening). Although the generation divergence at resolution jumps can by avoided by advancing in time the magnetic vector potential and computing the magnetic field as its curl (Etienne et al. 2010), the method is sensitive to the gauge condition employed, giving especially bad results for gauges with zero-speed modes (Etienne et al. 2012).

Furthermore, even though cell-centered versions of these methods exist (Toth 2000) and have been applied in GRMHD codes (see e.g., Gammie et al. 2003; Mizuno et al. 2006; Porth et al. 2017), these are incompatible with AMR, since the problem of finding divergence-preserving prolongation and restriction operation for cell-centered magnetic fields is under-determined, that is, information is lost when abandoning the face representation in favor of the cell-centered representation (Tóth & Roe 2002).

The approach used in this work consists in applying divergence-preserving restriction and prolongation operators for face-allocated magnetic fields. Examples of these operators were derived independently by Tóth & Roe (2002) and Balsara (2001), and have been employed in codes such as BATSRUS (Tóth et al. 2012) and AstroBEAR (Cunningham et al. 2009).

In addition, staggered magnetic fields allow the use of upwind magnetic field evolution schemes. In fact, as shown by Flock et al. (2010), not-upwind algorithms such as the widely used flux-CT and the Balsara & Spicer (1999) method (BS) are prone to numerical instabilities. Upwind CT schemes include the second-order algorithm by Gardiner & Stone (2005), used for example in Athena++ (White et al. 2016), and the two algorithms presented in Londrillo & Del Zanna (2004) and in Del Zanna et al. (2007), which allow the use of high-order schemes for the integration of Faraday’s law.

This work largely focuses on the AMR-compatible implementation of staggered grids and the upwind CT methods by Londrillo & Del Zanna (2004) and Del Zanna et al. (2007) in BHAC. Special attention is given to a new derivation of divergence-preserving prolongation operators, which constitute a coordinate-independent generalization of those found in the literature for the special case of Cartesian geometry (Tóth & Roe 2002).

Though these methods are fully general, technical details on their implementation in BHAC are provided, with the purpose of documenting how they can be incorporated in a GRMHD code. These include the ghost-cell exchange, the data structures and the special treatments employed at coarse-fine interfaces, and at the polar axis. The new methods implemented in the code are validated through several tests, which are also used to highlight the differences between the newly implemented schemes and those already present and validated. In particular, we present simulations of magnetised accretion onto a Kerr black hole, which take advantage of the newly implemented methods. This problem is of central importance for the interpretations of future images to be obtained by the EHT, and constitutes the main scientific application of the code. As shown later, we find that the divergence control technique employed significantly affects the dynamics of this system.

Moreover, the use of AMR makes it possible to simulate this system in Cartesian coordinates while still resolving the growth of the MRI with sufficient accuracy to reproduce the same mass-accretion rate and fluxes through the horizon found in spherical coordinates. This might be useful to avoid and possibly quantify directional biases introduced by grids in spherical coordinates, such as those commonly used in simulations of the accretion process, or eliminate the need of special boundary conditions in the vicinity of the polar axis.

This work is organised as follows: Sect. 2 presents the formulation of the equations of GRMHD employed in BHAC, as well as a brief description of the coordinates and equations of state available in the code. Section 3 summarizes the numerical methods used to solve these equations, namely finite-volume (Sect. 3.1) and CT methods (Sect. 3.2). The mathematical and technical aspects of BHAC’s AMR-compatible implementation of staggered grids is contained in Sect. 4, with Sect. 4.2 focusing on the new prolongation operators. Section 5 presents the problems used to validate the new methods, while the simulations of accretion onto black holes using the newly implemented methods are described in Sect. 6. Finally, in Sect. 7 we summarize and discuss the prospects of future research. Hereafter, we use Greek symbols for spacetime indices that run from zero to three and Latin symbols for for hypersurface indices, that run form one to three.

2. Mathematical formulation

2.1. Equations of GRMHD

BHAC was initiated as an extension of the MPI-AMRVAC framework toward GRMHD simulations in one, two, and three dimensions using finite-volume methods and a variety of modern numerical methods, described more in detail in Porth et al. (2017). It exploits much of MPI-AMRVAC’s infrastructure for parallelization and block-based automated adaptive mesh refinement (AMR), resulting in a potentially significant saving in computational time and computing resources. BHAC solves the equations of ideal magnetohydrodynamics (MHD) in general relativity but also in more generic spacetimes (Mizuno et al. 2018). The GRMHD equations are expressed as conservation equations for the rest-mass density, the energy and momentum, and the homogeneous Maxwell equations

$$\begin{array}{cc}& {\mathrm{\nabla }}_{\mu }(\rho u^{\mu })=0,\hfill \\ \hfill & {\mathrm{\nabla }}_{\mu }{T}^{\mu \nu }=0,\hfill \\ \hfill & {\mathrm{\nabla }}_{\mu }{}^{\ast }{F}^{\mu \nu }=0,\hfill \end{array}$$(1)

where ρ is the particle number density in the fluid frame, u^μ is the fluid 4-velocity in the observer frame, T^μν is the energy-momentum tensor and ^*​F^μν is the dual of the Faraday tensor.

The Faraday tensor F^μν and its dual ^*​F^μν are such that the electric and magnetic fields measured by an observer moving at 4-velocity U^μ are

$$\begin{array}{c}\hfill E^{\mu }:=F^{\mu \nu }{U}_{\nu }\mathrm{and}{B}^{\mu }:{=}^{\ast }{F}^{\mu \nu }{U}_{\nu },\end{array}$$(2)

and are related as ${}^{\ast }{F}^{\mu \nu }:=\frac{1}{2}{ϵ}^{\mu \nu \alpha \beta }{F}_{\alpha \beta }$ where ϵ^μναβ are the components of the Levi-Civita tensor. Therefore, in a general frame moving with 4-velocity U^μ

$$\begin{array}{c}\hfill F^{\mu \nu }:=U^{\mu }{E}^{\nu }-U^{\nu }{E}^{\mu }-{ϵ}^{\mu \nu \lambda \delta }{U}_{\lambda }{B}_{\delta },\end{array}$$(3)

and

$$\begin{array}{c}\hfill {}^{\ast }{F}^{\mu \nu }:=U^{\mu }{B}^{\nu }-U^{\nu }{B}^{\mu }+{ϵ}^{\mu \nu \lambda \delta }{U}_{\lambda }{E}_{\delta }.\end{array}$$(4)

Instead of solving the inhomogeneous Maxwell equations, we close the system by enforcing the ideal-MHD condition

$$\begin{array}{c}\hfill e^{\mu }=F^{\mu \nu }{u}_{\nu }=0,\end{array}$$(5)

where e^μ denotes the electric field in the fluid frame. Consequently, the magnetic field in the fluid frame is b^μ=^*​F^μνu_ν and

$$\begin{array}{c}\hfill {}^{\ast }{F}^{\mu \nu }=b^{\mu }{u}^{\nu }-b^{\nu }{u}^{\mu }.\end{array}$$(6)

Physically, the ideal-MHD condition corresponds to a plasma with an infinite conductivity and has the important consequence that in this case the magnetic field is simply advected with the fluid (frozen-flux theorem).

The energy-momentum tensor T^μν includes both fluid and electromagnetic contributions. Using Eq. (6) to write its electromagnetic part in terms of b^μ only (Anile 1990), it reads

$$\begin{array}{c}\hfill T^{\mu \nu }=\rho h_{\mathrm{tot}}{u}^{\mu }{u}^{\nu }+p_{\mathrm{tot}}{g}^{\mu \nu }-b^{\mu }{b}^{\nu },\end{array}$$(7)

where h_tot := h + b²/ρ is the total specific enthalpy and p_tot = p + b²/2 is the total pressure, both in the fluid frame.

The system is closed once an equation of state relating the specific enthalpy of the fluid with the density and the pressure h = h(ρ, p) is specified. The equations of state for an ideal gas, a Synge gas and a polytropic fluid are currently implemented in BHAC (see Porth et al. 2017, for more details), although the implementation of any additional one is straightforward as long as it is analytic. Additionally, BHAC can be used together with the microphysics routines developed in Most et al. (2019) that can handle temperature-dependent equations of state and provide a neutrino leakage scheme (Ruffert et al. 1996; O’Connor & Ott 2010; Galeazzi et al. 2013).

In order to formulate system (1) as a set of evolution equations, we employ the 3+1 decomposition of the spacetime (see e.g., Alcubierre 2008; Rezzolla & Zanotti 2013). The metric is decomposed as g_μν = γ_μν − n_μn_ν, where n_μ and γ_μν are such that γ_μ ^νn_ν = 0 and n_νn^ν = −1. Thus, γ_μ ^ν defines a projection operator onto a hypersurface normal to the vector field n^ν and contains the 3-metric induced in the hypersurface.

The explicit form of the line element in a 3+1 split spacetime is

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\gamma }_{\mathit{ij}}(\mathrm{d}{x}^i+{\beta }^i\mathrm{d}t)(\mathrm{d}{x}^j+{\beta }^j\mathrm{d}t),\end{array}$$(8)

where α and βⁱ are the lapse and the shift vector, respectively. Using γ_μ ^ν and n_ν, it is possible to define variables suitable to be evolved. These are the conserved variables in the Eulerian frame: the number density D := −ρu^νn_ν, the covariant 3-momentum S​_i := n_μγ_νiT^μν, the total energy U := n_μn_νT^μν and the spatial stress tensor W^ij := γ_μiγ_νjT^μν.

In terms of these variables, system (1) can be written as a set of conservation equations

$$\begin{array}{c}\hfill {\partial }_t(\sqrt{\gamma }\mathit{U})+{\partial }_i(\sqrt{\gamma }{\mathit{F}}^i)=\sqrt{\gamma }\mathit{S},\end{array}$$(9)

in addition to the solenoidal constraint for the magnetic field ${\partial }_i\sqrt{\gamma }{B}^i=0$. Here, γ is the determinant of the induced three-metric, and the vectors of conserved quantities U, the fluxes Fⁱ, and the sources S are given by

$$\begin{array}{cc}\hfill \mathit{U}& =\left[\begin{array}{c}D\\ S_j\\ \tau \\ B^j\end{array}\right],{\mathit{F}}^i=\left[\begin{array}{c}{\mathcal{V}}^{i}D\\ \alpha W_j^i-{\beta }^i{S}_j\\ \alpha (S^i-{\nu }^{i}D)-{\beta }^i\tau \\ {\mathcal{V}}^i{B}^j-B^i{\mathcal{V}}^j\end{array}\right],\hfill \\ \hfill \mathit{S}& =\left[\begin{array}{c}0\\ \frac{1}{2}\alpha W^{\mathit{ik}}{\partial }_j{\gamma }_{\mathit{ik}}+S_i{\partial }_j{\beta }^i-U{\partial }_j\alpha \\ \frac{1}{2}{W}^{\mathit{ik}}{\beta }^j{\partial }_j{\gamma }_{\mathit{ik}}+W_i^j{\partial }_j{\beta }^i-S^j{\partial }_j\alpha \\ 0\end{array}\right],\hfill \end{array}$$(10)

where ${\nu }^i:={\gamma }_{\mu }^i{u}^{\mu }/\mathrm{\Gamma }$ is the 3-velocity, Γ := −u^μn_μ is the Lorentz factor², and 𝒱ⁱ := ανⁱ − βⁱ is the transport velocity. Evolving τ := U − D instead of U makes the evolution more accurate in regions of low energy and allows to recover the Newtonian limit.

To calculate the fluxes Fⁱ, knowledge of the primitive variables P = [ρ,Γνⁱ,p,Bⁱ] is required. While it is straightforward to find U(P), P(U) requires numerical inversion. The inversion process then consists of finding the auxiliary variables A = [Γ,ξ], where ξ := Γ²ρh and h is the specific enthalpy compatible with U and P, which is in turn used to find the primitives. Once A is found, BHAC stores it in order to facilitate new inversions, thus extending the array U(P), as detailed in Porth et al. (2017).

2.2. Coordinates in BHAC

The main target application of BHAC is the simulation of accretion onto compact objects in arbitrary spacetimes. This has allowed to simulate accretion onto Kerr black holes to build theoretical models of M 87 (Davelaar et al. 2019), but also to explore the consequences of alternative theories of gravity (Mizuno et al. 2018) or of the presence of a boson star (Olivares et al. 2018) at the Galactic Center in the forthcoming horizon-scale VLBI images of Sgr A*, as well as the study of quasi-periodic oscillations in accretion disks around neutron stars (de Avellar et al. 2018).

As a result, BHAC is designed with a great flexibility to adopt new spacetime metrics. Its modular structure not only allows us to add straightforwardly new analytic coordinates, but another interesting feature is that it is not compulsory to provide derivatives of the metric, which are used to calculate the sources. The user needs only to provide the metric functions, and derivatives can be calculated internally using the complex-step derivative (Squire & Trapp 1998), a highly accurate numerical method which is capable of achieving machine precision for derivatives of algebraic functions (Martins et al. 2003).

The code can handle spacetime metrics given either in 3+1 form (i.e., expressed in terms of γ_ij, β_i and α) or as the full 4-metric g_μν, as the necessary conversions to 3+1 form are done internally. In addition, it is able to read tabulated numerical metrics when the spacetime is known only numerically. A list of metrics currently available in the code is reported in Table 1, which complements that presented in Porth et al. (2017).

The need to store the metric components both at barycenter positions and at cell interfaces is handled efficiently by adopting the metric data structure described in Porth et al. (2017), which takes advantage of the symmetry of the metric components and the fact that the metric and its derivatives are usually sparsely populated, that is, many of the tensor components are zero.

Due to its use in later parts of this work, we next describe the less well known modified Kerr-Schild coordinates (MKS). These coordinates, from now on labeled as x^μ = (t_MKS, s, ϑ, ϕ_MKS), were introduced by McKinney & Gammie (2004) in order to simulate large domains and, at the same time, concentrate the resolution near the black hole and on the equatorial plane. They are related to the standard Kerr-Schild coordinates x^μ′  =  (t_KS, r_KS, θ_KS, ϕ_KS) by the transformation

$$\begin{array}{cc}& t_{\mathrm{KS}}=t_{\mathrm{MKS}},\hfill \\ \hfill & r_{\mathrm{KS}}=R_0+e^{\mathrm{s}},\hfill \\ \hfill & {\theta }_{\mathrm{KS}}=\vartheta +\frac{{\vartheta }_0}{2}sin(2\vartheta ),\hfill \\ \hfill & {\varphi }_{\mathrm{KS}}={\varphi }_{\mathrm{MKS}},\hfill \end{array}$$(11)

where R₀ and ϑ₀ are constant parameters. The convention used here for the parameter ϑ₀ is different from that in McKinney & Gammie (2004). The one employed here is chosen so that ϑ reduces to θ_KS when ϑ₀ = 0, thus retrieving logarithmic Kerr-Schild coordinates.

It is worth to mention that, in contrast to Porth et al. (2017), in this work the inverse relation ϑ(θ_KS) is not approximated using a cubic polynomial, but the resulting transcendental equation is solved numerically whenever necessary.

3. Numerical methods

3.1. Finite-volume scheme

To introduce the notation used in this work, we next summarize the finite-volume scheme used by BHAC to evolve the hydrodynamic variables in system (9). This scheme was already described in Porth et al. (2017), to which we refer the reader for more details.

To obtain the discretized equations, the simulation domain is divided into control volumes and the system is integrated over each of them. This leads to evolution equations for the average ${\overline{\mathit{U}}}_{i,j,k}$ of the conserved quantities inside each cell

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\overline{\mathit{U}}}_{i,j,k}}{\mathrm{d}t}=& -\frac{1}{\mathrm{\Delta }{V}_{i,j,k}}[{\mathit{F}}^{\mathbf{1}}\mathrm{\Delta }{S}^1{|}_{i+1/2,j,k}-{\mathit{F}}^{\mathbf{1}}\mathrm{\Delta }{S}^1{|}_{i-1/2,j,k}\hfill \\ \hfill & +{\mathit{F}}^{\mathbf{2}}\mathrm{\Delta }{S}^2{|}_{i,j+1/2,k}-{\mathit{F}}^{\mathbf{2}}\mathrm{\Delta }{S}^2{|}_{i,j-1/2,k}\hfill \\ \hfill & +{\mathit{F}}^{\mathbf{3}}\mathrm{\Delta }{S}^3{|}_{i,j,k+1/2}-{\mathit{F}}^{\mathbf{3}}\mathrm{\Delta }{S}^3{|}_{i,j,k-1/2}]+{\overline{\mathit{S}}}_{i,j,k}.\hfill \end{array}$$(12)

Quantities indicated as, say, F¹ΔS​¹|_{i + 1/2, j, k} represent integrals of the fluxes over the surfaces ΔS​¹|_{i + 1/2, j, k} bounding the control volume, and ${\overline{\mathit{S}}}_{i,j,k}$ is the volume average of the sources. Both kinds of integrals are approximated to second order, by assigning to Fⁿ (n = 1, 2, 3) the point value of the flux at the face center and to ${\overline{\mathit{S}}}_{i,j,k}$ the point value at the cell barycenter, that is, at

$$\begin{array}{cc}& x_{(\mathrm{face}k)}^i=x^i+{\delta }_k^i\frac{\mathrm{\Delta }{x}^i}{2},\hfill \\ \hfill & x_{(\mathrm{barycenter})}^i=\frac{{\int }_{\mathrm{cell}}\sqrt{\gamma }{x}^i\mathrm{d}{x}^1\mathrm{d}{x}^2\mathrm{d}{x}^3}{{\int }_{\mathrm{cell}}\sqrt{\gamma }\mathrm{d}{x}^1\mathrm{d}{x}^2\mathrm{d}{x}^3},\hfill \end{array}$$

respectively, where Δxⁱ is the grid spacing in each direction and ${\delta }_k^i$ is Kronecker’s delta. Fⁿ is obtained through the approximate solution of a Riemann problem at the interface, and static integrals such as cell volumes, interface areas and barycenter positions are calculated at initialization using fourth-order Simpson’s rule and stored in memory.

The Riemann solvers currently available in BHAC are the Rusanov method, also known as total variation diminishing Lax-Friedrichs scheme, and the HLL solver of Harten et al. (1983). The left and right states for the Riemann problem are obtained using limited reconstructions from the cell centers (see also Toro 2009; Rezzolla & Zanotti 2013).

BHAC features a variety of reconstruction schemes, some of which are listed in Keppens et al. (2012), Porth et al. (2017). They include both second-order, total variation diminishing (e.g., minmod, vanLeer), and high-order such as PPM (Colella & Woodward 1984), LIMO3 (Čada & Torrilhon 2009) and MP5 (Suresh & Huynh 1997). Recent additions to that list are the high-order schemes WENO5 and WENOZ+ (Acker et al. 2016). Equation (12) can be integrated using any of the methods present in the MPI-AMRVAC toolkit. These include the simple half-step modified Euler, the third order Runge-Kutta RK3 (Gottlieb & Shu 1998) and the strong-stability preserving s-step, pth-order RK schemes SSPRK (s,p) schemes: SSPRK (4,3), SSPRK (5,4) due to Spiteri & Ruuth (2002) (see Porth et al. 2014, for implementation details).

3.2. Constrained transport

Constrained transport (CT) is a divergence-control method first proposed by Evans & Hawley (1988). In essence, instead of eliminating the divergence of the magnetic field once it is created, it modifies the way in which magnetic-field transport is evolved so as to prevent the creation of divergence. CT is able to keep a discretization of ∇ ⋅ B = 0 to a precision close to that of floating point operations by ensuring that the sum of the magnetic fluxes through the surfaces bounding a cell is zero to machine precision. Recalling the definition of the divergence of a vector field B

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{B}=\underset{\mathrm{\Delta }V\to 0}{lim}\frac{1}{\mathrm{\Delta }V}{{\displaystyle \oint }}_{S=\partial \mathrm{\Delta }V}\mathit{B}·\mathit{n}\mathrm{d}S,\end{array}$$(13)

is equivalent to ∇ ⋅ B = 0 in the continuous limit. When the limit is not taken, that is, for the finite-volume case, it follows from the divergence theorem that the average value of the divergence ${\left(\overline{\mathrm{\nabla }·\mathit{B}}\right)}_{\mathrm{cell}}$ is zero within the cell.

The central idea from constrained transport is to give the electromagnetic variables a special spatial location, as depicted in Fig. 1. In particular, to each face of the cell corresponds a magnetic flux calculated, for example, as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{i+1/2,j,k}={{\displaystyle \int }}_{\partial V(x_{i+1/2}^1)}{\gamma }^{1/2}{B}^1\mathrm{d}{x}^2\mathrm{d}{x}^3,\end{array}$$(14)

Fig. 1. Spatial location of variables for a cell with indices (i, j, k), correspondent to directions (x, y, z), respectively. Line integrals of the electric field ℰ are located at its edges, and magnetic and numerical fluxes Φ and Fi (the latter used for the BS algorithm) are located at its faces. The conserved hydrodynamic variables belonging to the array U, are located at cell centers (not shown).

and on each edge is associated a line integral of the electric field, similar to

$$\begin{array}{c}\hfill {\mathcal{E}}_{i+1/2,j+1/2,k}=-{{\displaystyle \int }}_{x_{k-1/2}^3}^{x_{k+1/2}^3}{E_3|}_{x_{i+1/2}^1,x_{j+1/2}^2}\mathrm{d}{x}^3.\end{array}$$(15)

The magnetic flux at each face is updated from the circulation of the electric field, using the integral form of Faraday’s law

$$\begin{array}{cc}\hfill \frac{d}{\mathrm{d}t}{\mathrm{\Phi }}_{i+1/2,j,k}=& {\mathcal{E}}_{i+1/2,j+1/2,k}-{\mathcal{E}}_{i+1/2,j-1/2,k}\hfill \\ \hfill & -{\mathcal{E}}_{i+1/2,j,k+1/2}+{\mathcal{E}}_{i+1/2,j,k-1/2}.\hfill \end{array}$$(16)

Since each of the line integrals of the electric field is shared by two faces, but appears with opposite sign in the time update formula for each of them, the rate of change of ${\left(\overline{\mathrm{\nabla }·\mathit{B}}\right)}_{\mathrm{cell}}$, which can be calculated as the sum of the rate of change of the outgoing flux through all faces, vanishes. Therefore, the CT time update ensures that ${\left(\overline{\mathrm{\nabla }·\mathit{B}}\right)}_{\mathrm{cell}}$ is kept constant at machine precision from one iteration to the next.

In order for B to be divergence-free during the whole simulation, ${\left(\overline{\mathrm{\nabla }·\mathit{B}}\right)}_{\mathrm{cell}}$ must be zero at the initial condition. This can be accomplished by initialising the line integrals of magnetic vector potential along cell edges and calculating the initial magnetic fluxes as the circulation of the vector potential around the cell faces, in the same way as the rate of change of the magnetic flux is calculated from the circulation of the electric field.

After the time update, the magnetic field is interpolated to the cell center in order to use it for the inversion from conserved to primitive variables.

The idea of staggering the magnetic components of the electromagnetic field to achieve a divergence-free evolution to machine precision was first proposed by Yee (1966). However, although his method was widely known in engineering, a staggered grid was not applied in GRMHD until the work of Evans & Hawley (1988).

In the formulas written so far, no approximations have been made. Approximations come into play when deciding how to calculate the line integrals of the electric field. The way these approximations are done distinguishes the different variants of CT. Three of these variants are described in the following sections.

3.3. Arithmetic averaging (BS)

This variant was introduced by Balsara & Spicer (1999). In the comparison work of Toth (2000), it is referred to as Flux-interpolated Constrained Transport (flux-CT), although in most of the literature this name is used exclusively for the cell-centered version of this method, proposed also in Toth (2000). In this work, we will use the abbreviation flux-CT to refer to the cell-centered version of the method and BS to refer to the staggered one. Both variants are particularly suitable for finite-volume schemes, since the electric field at cell edges is estimated as the arithmetic average of the numerical fluxes returned by the Riemann solver. In fact, the numerical fluxes corresponding to the magnetic-field components are surface integrals of the electric field, for example, the flux in the x²-direction for B¹ is

$$\begin{array}{c}\hfill {\mathrm{\Delta }{S}^2{\overline{F}}^2|}_{i,j+1/2,k}={{\displaystyle \int }}_{x_{i-1/2}^1}^{x_{i+1/2}^1}{{\displaystyle \int }}_{x_{k-1/2}^3}^{x_{k+1/2}^3}{E_{x^3}|}_{j+1/2}\mathrm{d}{x}^3\mathrm{d}{x}^1.\end{array}$$(17)

The innermost integral is the same as that of Eq. (15), so the average flux can be interpreted as

$$\begin{array}{c}\hfill {\mathrm{\Delta }{S}^2{\overline{F}}^2|}_{i,j+1/2,k}=-\mathrm{\Delta }{x}_i{\stackrel{\sim }{\mathcal{E}}}_{i,j+1/2,k},\end{array}$$(18)

where ${\stackrel{\sim }{\mathcal{E}}}_{i,j+1/2,k}$ is the mean value of the integral from Eq. (15) over the face at j + 1/2. To second-order accuracy, this integral takes the value ${\stackrel{\sim }{\mathcal{E}}}_{i,j+1/2,k}$ at the middle of the cell. Therefore, ℰ_{i + 1/2, j + 1/2, k} can be found by interpolating the averaged fluxes from the four adjacent cell faces as

$$\begin{array}{cc}\hfill {\mathcal{E}}_{i+1/2,j+1/2,k}=& \frac{1}{4}(\frac{\mathrm{\Delta }{S}^2{\overline{F}}^2{|}_{i,j+1/2,k}}{\mathrm{\Delta }{x}_i}+\frac{\mathrm{\Delta }{S}^2{\overline{F}}^2{|}_{i+1,j+1/2,k}}{\mathrm{\Delta }{x}_{i+1}}\hfill \\ \hfill & -\frac{\mathrm{\Delta }{S}^1{\overline{F}}^1{|}_{i+1/2,j,k}}{\mathrm{\Delta }{y}_j}-\frac{\mathrm{\Delta }{S}^1{\overline{F}}^1{|}_{i+1/2,j+1,k}}{\mathrm{\Delta }{y}_{j+1}})·\hfill \end{array}$$(19)

Although this algorithm, especially in its cell-centered version, is widely used in the community (see e.g., Gammie et al. 2003; Noble et al. 2009), it is known to lack a proper upwinding and to not reduce to the correct limit for one-dimensional flow (Gardiner & Stone 2005), as well as to cause numerical instabilities (Flock et al. 2010). Efforts to overcome this problems motivated the development of methods such as those by Gardiner & Stone (2005), Londrillo & Del Zanna (2004), Del Zanna et al. (2007). In the next section, we will briefly describe the methods by Londrillo & Del Zanna (2004) and Del Zanna et al. (2007), which we have also implemented in BHAC.

3.4. Upwind constrained transport (UCT)

Upwind Constrained Transport (UCT) is a method to evolve the magnetic field, first proposed in Londrillo & Del Zanna (2004). As a CT method, like BS, it maintains the divergence of the magnetic field to machine precision; however, it furthermore aims to accurately reproduce the magnetic-field continuity and transport properties by using limited reconstructions. Also, in contrast to arithmetic averaging, UCT reduces to the correct 1-dimensional limit when symmetry in the other two directions is assumed. Two variants of UCT are implemented in BHAC. These are presented in Londrillo & Del Zanna (2004) and in Del Zanna et al. (2007), and will be referred, respectively, as UCT1 and UCT2. Now we describe the procedure to obtain ℰ in each of these algorithms. To simplify the notation, we rename (x¹, x², x³) as (x, y, z). We focus on the calculation of ℰ_z, but the other cases can be obtained by iteratively replacing x → y, y → z and z → x (see Fig. 2). In BHAC, the implementation of UCT1 is as follows:

Fig. 2. Reconstruction of quantities for UCT1 (left) and UCT2 (right). Quantities in blue are defined at cell faces in x and y and quantities in red at the cell edge in z. The resulting line integral of the electric field ℰz over the edge is shown in green. Continuous arrows represent limited reconstructions (steps 2 and 3 of UCT1, and 3 of UCT2) and dashed lines the choice of the maximum value (step 5 of UCT1 and 4 of UCT2, the labels L and R for the characteristic speeds are omitted in the figure to avoid saturation). The average transverse transport velocities needed in UCT2 (step 2) ${\overline{\mathcal{V}}}^{x,\mathit{y}}$ are computed using the quantities in the boxes.

1. At the time of calculating the numerical fluxes at x-interfaces, store the characteristic speeds $c_i^{\text{min}}$ and $c_i^{\text{max}}$ (i = x, y, z) obtained from the Riemann solver, as well as the transport velocity ${\mathcal{V}}_{\text{L}}^x$, ${\mathcal{V}}_{\mathrm{L}}^{\mathit{y}}$, ${\mathcal{V}}_{\text{R}}^x$, and ${\mathcal{V}}_{\mathrm{R}}^{\mathit{y}}$.

2. Reconstruct transport velocities along direction y toward the edges z and obtain one correspondent to the four states around them ${\mathcal{V}}_{\text{LL}}^i$, ${\mathcal{V}}_{\text{LR}}^i$, ${\mathcal{V}}_{\text{RL}}^i$ and ${\mathcal{V}}_{\text{RR}}^i$ (i = x, y).

3. Reconstruct the magnetic fields on the faces to the same edges and obtain $B_{\text{L}}^x$, $B_{\text{R}}^x$, $B_{\mathrm{L}}^{\mathit{y}}$ and $B_{\mathrm{R}}^{\mathit{y}}$.

4. Calculate approximations to the electric field on the edge from each of the possible states as

   $$\begin{array}{cc}\hfill E_z^{\mathrm{LL}}& ={\mathcal{V}}_{\mathrm{LL}}^y{B}_{\mathrm{L}}^x-{\mathcal{V}}_{\mathrm{LL}}^x{B}_{\mathrm{L}}^y,\hfill \\ \hfill E_z^{\mathrm{LR}}& ={\mathcal{V}}_{\mathrm{LR}}^y{B}_{\mathrm{R}}^x-{\mathcal{V}}_{\mathrm{LR}}^x{B}_{\mathrm{L}}^y,\hfill \\ \hfill E_z^{\mathrm{RL}}& ={\mathcal{V}}_{\mathrm{RL}}^y{B}_{\mathrm{L}}^x-{\mathcal{V}}_{\mathrm{RL}}^x{B}_{\mathrm{R}}^y,\hfill \\ \hfill E_z^{\mathrm{RR}}& ={\mathcal{V}}_{\mathrm{RR}}^y{B}_{\mathrm{R}}^x-{\mathcal{V}}_{\mathrm{RR}}^x{B}_{\mathrm{R}}^y.\hfill \end{array}$$

5. Take the maximum characteristic speeds in each direction from the four states,

   $$\begin{array}{cc}\hfill c_x^{-}& =max(c_{x,\mathrm{L}}^{\min },c_{x,\mathrm{R}}^{\min }),\hfill \\ \hfill c_x^{+}& =max(c_{x,\mathrm{L}}^{\max },c_{x,\mathrm{R}}^{\max }),\hfill \\ \hfill c_y^{-}& =max(c_{y,\mathrm{L}}^{\min },c_{y,\mathrm{R}}^{\min }),\hfill \\ \hfill c_y^{+}& =max(c_{y,\mathrm{L}}^{\max },c_{y,\mathrm{R}}^{\max }),\hfill \end{array}$$

   where $c_{i,\text{L}(\text{R})}^{\text{min}}$ $\left(c_{i,\text{L}(\text{R})}^{\text{max}}\right)$ are defined as positive in the −i (+i) direction, and they are zero otherwise.

6. Finally, approximate the electric field, as an average of those at the four states weighted by the characteristic speeds and corrected by the transport of magnetic field, using the formula by Londrillo & Del Zanna (2004):

   $$\begin{array}{cc}\hfill E_z=& \sqrt{\gamma }[\frac{c_x^{+}{c}_y^{+}{E}_z^{\mathrm{LL}}+c_x^{+}{c}_y^{-}{E}_z^{\mathrm{LR}}+c_x^{-}{c}_y^{+}{E}_z^{\mathrm{RL}}+c_x^{-}{c}_y^{-}{E}_z^{\mathrm{RR}}}{(c_x^{+}+c_x^{-})(c_y^{+}+c_y^{-})}\hfill \\ \hfill & +\frac{c_x^{+}{c}_x^{-}}{c_x^{+}+c_x^{-}}(B_y^{\mathrm{R}}-B_y^{\mathrm{L}})-\frac{c_y^{+}{c}_y^{-}}{c_y^{+}+c_y^{-}}(B_x^{\mathrm{R}}-B_x^{\mathrm{L}})]·\hfill \end{array}$$(20)

   Since this is a point value, we take the second order approximation that ℰ_z = E_zΔz and use it to calculate the circulation.

The implementation of UCT2 is as follows:

1. At the time of calculating the numerical fluxes at each interface, store the characteristic speeds $c_i^{\text{min}}$ and $c_i^{\text{max}}$ (i = x, y, z) obtained from the Riemann solver, and the transverse transport velocities at the left and right states, for example ${\mathcal{V}}_{\mathrm{L}}^{\mathit{y}}$, ${\mathcal{V}}_{\text{L}}^z$, ${\mathcal{V}}_{\mathrm{R}}^{\mathit{y}}$, and ${\mathcal{V}}_{\text{R}}^z$ for the x-interface.

2. Obtain a new transverse transport velocity weighting those of the left and right states by the characteristic speeds (Del Zanna et al. 2007). For the x-interface, this is

   $$\begin{array}{c}\hfill {\overline{\mathcal{V}}}^{y,z}=\frac{c_x^{\max }{\mathcal{V}}_{\mathrm{L}}^{y,z}+c_x^{\min }{\mathcal{V}}_{\mathrm{R}}^{y,z}}{c_x^{\max }+c_x^{\min }}·\end{array}$$(21)

3. Reconstruct the magnetic fields and the transport velocities to the edge where ℰ_z needs to be calculated. This gives $B_{\mathrm{L},\mathrm{R}}^{x,\mathit{y}}$ and ${\overline{\mathcal{V}}}_{\mathrm{R},\mathrm{L}}^{x,\mathit{y}}$.

4. Take the maximum characteristic speeds in the same way as in step (5) of UCT1.

5. Approximate the electric field using the formula by (Del Zanna et al. 2007):

   $$\begin{array}{cc}\hfill E_z=& \sqrt{\gamma }[-\frac{c_x^{+}{\overline{\mathcal{V}}}_{\mathrm{L}}^x{B}_{\mathrm{L}}^y+c_x^{-}{\overline{\mathcal{V}}}_{\mathrm{R}}^x{B}_{\mathrm{R}}^y-c_x^{-}{c}_x^{+}(B_{\mathrm{R}}^y-B_{\mathrm{L}}^y)}{c_x^{+}+c_x^{-}}\hfill \\ \hfill & +\frac{c_y^{+}{\overline{\mathcal{V}}}_{\mathrm{L}}^y{B}_{\mathrm{L}}^x+c_x^{-}{\overline{\mathcal{V}}}_{\mathrm{R}}^y{B}_{\mathrm{R}}^x-c_y^{-}{c}_y^{+}(B_{\mathrm{R}}^x-B_{\mathrm{L}}^x)}{c_y^{+}+c_y^{-}}]·\hfill \end{array}$$(22)

   Again, we approximate ℰ_z = E_zΔz and use it to calculate the circulation.

Since UCT2 appears to be more symmetric than UCT1, we prefer it for the simulations presented in this work. However, it is worth to mention that we did not observe important differences, especially no directional bias, when using with UCT1.

An important difference in our implementation with respect to that of Londrillo & Del Zanna (2004) and Del Zanna et al. (2007) is that in those works the quantities $\sqrt{\gamma }{B}^i$ and not the magnetic fields Bⁱ are reconstructed at the edges. We instead, factor the square root of the metric determinant and multiply the whole expression for E_z by it at the end. This makes the scheme more consistent when evolving a radial magnetic field in the polar regions, and simplifies the AMR operations at the poles, as will be seen in Sect. 4.5.

4. CT-adapted AMR

Constrained transport schemes ensure that, during evolution, no divergence of the magnetic field is created from one iteration to the next, at machine precision. However, by construction, no mechanism to eliminate divergence that has already been created is present. Therefore, in order to exploit the advantages of constrained transport in an AMR code, it is necessary to resort to prolongation and restriction operators that also preserve the constraint ∇ ⋅ B = 0 to machine precision. In addition, care must be taken to synchronise different representations of the same electric and magnetic field components across fine-coarse interfaces, which is done in a step similar to the refluxing of finite-volume methods. After defining the notation in Sect. 4.1, we here describe in detail such operations.

4.1. Notation

To avoid using lengthy subscripts, we employ a notation similar to that used by Tóth & Roe (2002). We refer to Fig. 3 for a schematic view.

Fig. 3. Notation for magnetic fluxes and line integrals of the electric field for Sect. 4.1.

The coarse coordinate increment in the direction i is denoted ΔXⁱ and the correspondent fine increment Δxⁱ = ΔXⁱ/2. The center of the coarse cell is defined at $(x_0^1,x_0^2,x_0^3)$ and those of the eight fine cells that result after refinement at $(x_0^1\pm \mathrm{\Delta }{x}^1/2,x_0^2\pm \mathrm{\Delta }{x}^2/2,x_0^3\pm \mathrm{\Delta }{x}^3/2)$. Quantities defined at the fine cell centers are labeled by the subscripts ±1, ±1, ±1, and those at the coarse cell center by 0, 0, 0, although this label is often omitted to keep the notation simple. Accordingly, ΔV is the volume of the coarse cell.

The coarse faces defined by the coordinate surfaces of x¹ are centered at $(x_0^1+\mathrm{\Delta }{X}^1/2,x_0^2,x_0^3)$. Quantities defined at these faces are labeled by the subscripts ±2, 0, 0. For instance, both these faces and their areas are denoted ΔS​_{±2, 0, 0}. The four fine faces in which ΔS​_2, 0, 0 is subdivided are labeled Δs_{2, ±1, ±1}. The same applies to ΔS​_{−2, 0, 0} and in the other directions. The twelve fine faces Δs_{0, ±1, ±1}, Δs_{±1, 0, ±1} and Δs_{±1, ±1, 0} are not part of any coarse face. The magnetic flux across Δs_2, 1, 1 is

$$\begin{array}{c}\hfill {\varphi }_{2,1,1}={{\displaystyle \int }}_{\mathrm{\Delta }{s}_{2,1,1}}{B}^1\sqrt{\gamma }\mathrm{d}{x}^2\mathrm{d}{x}^3,\end{array}$$(23)

and that across ΔS​_2, 0, 0 is

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{2,0,0}={{\displaystyle \int }}_{\mathrm{\Delta }{S}_{2,0,0}}{B}^1\sqrt{\gamma }\mathrm{d}{x}^2\mathrm{d}{x}^3.\end{array}$$(24)

The same applies to quantities defined at cell edges: ℰ_2, 2, 0 denotes line integral of the electric field along the coarse edge centered at $(x_0^1+\mathrm{\Delta }{X}^1/2,x_0^2+\mathrm{\Delta }{X}^2/2,x_0^3)$ and ℰ_2, 2, 1 that along the fine edge centered at $(x_0^1+\mathrm{\Delta }{X}^1/2,x_0^2+\mathrm{\Delta }{X}^2/2,x_0^3+\mathrm{\Delta }{x}^3/2)$.

In a 2D problem, where symmetry in the x³ direction is assumed, the integrals with respect to x³ can be ignored, since they only result in multiplying all quantities by a constant factor that is cancelled in every equation.

The divergence of the magnetic field in the coarse cell is discretized as

$$\begin{array}{c}\hfill (\mathrm{\nabla }·\mathit{B})=\frac{1}{\mathrm{\Delta }V}({\mathrm{\Phi }}_{2,0,0}-{\mathrm{\Phi }}_{-2,0,0}+{\mathrm{\Phi }}_{0,2,0}-{\mathrm{\Phi }}_{0,-2,0}+{\mathrm{\Phi }}_{0,0,2}-{\mathrm{\Phi }}_{0,0,-2}),\end{array}$$(25)

which reduces to the standard definition of divergence (Eq. (13)) in the limit ΔV → 0, and corresponding definitions are valid for each of the fine cells.

4.2. Prolongation and restriction operators for face-allocated variables

4.2.1. Restriction formulas

As in Tóth & Roe (2002), to obtain a restriction formula that naturally preserves the divergence, we simply require that the flux through one of the coarse faces is the sum of the fluxes through the fine faces which form part of it, that is,

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\pm 200}& ={\displaystyle \sum _{j,k=-1,1}}{\varphi }_{\pm 2jk},\hfill \\ \hfill {\mathrm{\Phi }}_{0\pm 20}& ={\displaystyle \sum _{i,k=-1,1}}{\varphi }_{i\pm 2k},\hfill \\ \hfill {\mathrm{\Phi }}_{00\pm 2}& ={\displaystyle \sum _{i,j=-1,1}}{\varphi }_{ij\pm 2}.\hfill \end{array}$$(26)

The divergence in the newly created coarse cell will be

$$\begin{array}{c}\hfill (\mathrm{\nabla }·\mathit{B})=\frac{1}{\mathrm{\Delta }V}{\displaystyle \sum _{i,j,k=-1,1}}\mathrm{\Delta }{V}_{i,j,k}{(\mathrm{\nabla }·\mathit{B})}_{i,j,k},\end{array}$$(27)

which is consistent with definition (25) since the interior fluxes cancel in pairs. If all (∇⋅B)_i, j, k are numerical zeros, Eq. (27) implies that (∇⋅B) is as well.

4.2.2. Prolongation formulas

While divergence-free restriction operators can be derived straightforwardly as seen above, it is not trivial to find prolongation operators with the same property. To the best of our knowledge, the only such prolongation operators for vector fields collocated at cell faces are those found by Balsara (2001) and Tóth & Roe (2002), who also derived similar operators for continuous interpolation and for prolongation of vertex-collocated fields. Since the divergence constraint alone results in an under-determined system of equations for the fine magnetic fields, Tóth & Roe (2002) close the system by requiring also that, on the fine grid, two different discretizations of the curl in Cartesian coordinates give the same numerical value. The formulas by Balsara (2001) are identical to those of divergence-preserving continuous interpolation from Tóth & Roe (2002), after making consistent the conventions used in both works for the cell dimensions. However, it is worth mentioning that Balsara (2001) obtained them from different requirements, namely by asking the interpolation of the magnetic field to be second-order, total-variation-diminishing (TVD), in addition to have vanishing divergence. The formulas by Tóth & Roe (2002) for face-collocated fields and those for continuous interpolation differ in that the former include additional terms arising from the requirement of constant discretized curl, while for the latter the curl is allowed to vary linearly inside the parent cell. Balsara (2004) provided a generalised version of his expressions for orthonormal cylindrical and spherical coordinates, by means of appropriate substitutions in the Cartesian formula, for instance, of the form B^x  →  rB^r, B^y  →  B^ϕ/r, B^z  →  B^z for cylindrical coordinates. The aforementioned formulas were successfully implemented and tested for non-relativistic MHD in the codes AstroBEAR (Cunningham et al. 2009) and BATSRUS (Tóth et al. 2012).

Due to the relation between the curl of the magnetic field and the electric current, asking the prolonged fields to preserve the value of the curl is well motivated from a physical point of view. However, this becomes more complicated when moving to a curved spacetime. On a general 3-dimensional hypersurface, the expressions for the divergence (∇ ⋅ B = γ^−1/2∂_iγ^1/2Bⁱ) and the curl ((∇×B)ⁱ = ϵ^ijk∂_jB_k) of a 3-vector not only involve additional quantities from the metric, but also operations of raising and lowering indices, or derivatives with respect to the covariant coordinates. For non-orthonormal bases, this requires to compute additional approximations to derivatives or interpolations of the magnetic-field components to faces in which they were not originally defined, thus increasing the computational cost and potentially decreasing the accuracy of the algorithm. For this reason, we revisit the problem of finding divergence preserving prolongation operators for face-allocated quantities using an alternative approach, and derive a general class of such formulas which reduces as a special case to the operator of Tóth & Roe (2002) for face-collocated fields.

Prolongation consists in two steps. First, interpolation is done for the fine fluxes which are part of a flux on the coarse faces (the exterior fluxes), and then the remaining (interior) fine fluxes are computed. Interpolation is done on the magnetic fluxes, not on the fields, and any formula that satisfies Eq. (26) can be used. The expressions employed in this work are

$$\begin{array}{cc}\hfill {\varphi }_{\pm 2jk}& =\omega ({\mathrm{\Phi }}_{\pm 200}+j\frac{\mathrm{\Delta }{x}^2}{2}{({\partial }_2\mathrm{\Phi })}_{\pm 200}+k\frac{\mathrm{\Delta }{x}^3}{2}{({\partial }_3\mathrm{\Phi })}_{\pm 200})\hfill \\ \hfill {\varphi }_{i\pm 2k}& =\omega ({\mathrm{\Phi }}_{0\pm 20}+i\frac{\mathrm{\Delta }{x}^1}{2}{({\partial }_1\mathrm{\Phi })}_{0\pm 20}+k\frac{\mathrm{\Delta }{x}^3}{2}{({\partial }_3\mathrm{\Phi })}_{0\pm 20})\hfill \\ \hfill {\varphi }_{ij\pm 2}& =\omega ({\mathrm{\Phi }}_{00\pm 2}+i\frac{\mathrm{\Delta }{x}^1}{2}{({\partial }_1\mathrm{\Phi })}_{00\pm 2}+j\frac{\mathrm{\Delta }{x}^2}{2}{({\partial }_2\mathrm{\Phi })}_{00\pm 2}),\hfill \end{array}$$(28)

where ω = 1/4 in 3D and 1/2 in 2D and (∂_iΦ) is any approximation to the slope of Φ. As Tóth & Roe (2002), to compute these slopes we choose second order formulas such as

$$\begin{array}{c}\hfill {({\partial }_2\mathrm{\Phi })}_{\pm 200}\approx \frac{1}{8\mathrm{\Delta }{x}^2}({\mathrm{\Phi }}_{\pm 2+40}-{\mathrm{\Phi }}_{\pm 2-40}).\end{array}$$(29)

However, it is also possible to use limited slopes, for example minmod, as suggested by Balsara (2001).

Once the fine magnetic fluxes on the exterior faces have been calculated, the next steps have the purpose of filling the interior fluxes without creating divergence. We looked for the most general linear formula that gave these fine interior fluxes in terms of the exterior ones, and which fulfilled three requisites: first, to keep constant a magnetic flux that was originally constant; second, to be reversible (specular symmetry); and third, to preserve the discretized divergence of the magnetic field in the coarse cell.

First, we write each of interior magnetic fluxes ϕ_l, m, n as a combination of the exterior ones and 24 coefficients n_i, j, k.

$$\begin{array}{cc}\hfill {\varphi }_{l,m,n}=& {\displaystyle \sum _{i=-2,2}}{\displaystyle \sum _{j,k=-1,1}}{n}_{l,m,n}^{i,j,k}{\varphi }_{i,j,k}+{\displaystyle \sum _{j=-2,2}}{\displaystyle \sum _{i,k=-1,1}}{n}_{l,m,n}^{i,j,k}{\varphi }_{i,j,k}\hfill \\ \hfill & +{\displaystyle \sum _{k=-2,2}}{\displaystyle \sum _{i,j=-1,1}}{n}_{l,m,n}^{i,j,k}{\varphi }_{i,j,k}.\hfill \end{array}$$(30)

The first requisite can be expressed as

$$\begin{array}{cc}& {\varphi }_{-2,j,k}={\varphi }_{2,j,k},{\varphi }_{i,-2,k}={\varphi }_{i,2,k},{\varphi }_{i,j,-2}={\varphi }_{i,j,2}\hfill \\ \hfill & \Rightarrow {\varphi }_{0,j,k}={\varphi }_{-2,j,k},{\varphi }_{i,0,k}={\varphi }_{i,-2,k},{\varphi }_{i,j,0}={\varphi }_{i,j,-2}.\hfill \end{array}$$(31)

Equation (31) is more general than simply requiring that a constant magnetic field is kept constant. It means that if the magnetic flux does not change after entering and exiting the cell, it will not change inside. This allows to prolong in a consistent way, for example, a monopolar field in spherical symmetry, for which the magnetic flux is not the same for all the cell faces at r = r₀, but for the same cell must not change between the face at r = r₀ and that at r = r₀ + Δr. Applied to the fluxes in the first direction, requirement Eq. (31) results in 12 equations for the coefficients of each flux ϕ_lmn, namely

$$\begin{array}{cc}\hfill n_{0,m,n}^{-2,j,k}+n_{0,m,n}^{2,j,k}& ={\delta }_m^j{\delta }_n^k,\hfill \\ \hfill n_{0,m,n}^{i,-2,k}+n_{0,m,n}^{i,2,k}& =0,\hfill \\ \hfill n_{0,m,n}^{i,j,-2}+n_{0,m,n}^{i,j,2}& =0,\hfill \end{array}$$(32)

where ${\delta }_j^i$ is the Kronecker delta. The next requirement, which we call reversibility, is the symmetry condition that when all the exterior fluxes are reflected with respect to the surface containing n_lmn, the sign of n_lmn must reverse. This results in another set of relations for the coefficients,

$$\begin{array}{cc}& n_{0,m,n}^{-2,j,k}=n_{0,m,n}^{2,j,k},\hfill \\ \hfill & n_{0,m,n}^{-1,\pm 2,k}=-n_{0,m,n}^{1,\pm 2,k},\hfill \\ \hfill & n_{0,m,n}^{-1,j,\pm 2}=-n_{0,m,n}^{1,j,\pm 2},\hfill \end{array}$$(33)

and analogous expression for the other interior fluxes. Combining with the first requirement (Eq. (32)), and removing the redundant relations, we can already find that

$$\begin{array}{cc}\hfill n_{0,m,n}^{-2,j,k}& =n_{0,m,n}^{2,j,k}=\frac{1}{2}{\delta }_m^j{\delta }_n^k,\hfill \\ \hfill n_{0,m,n}^{1,2,k}& =-n_{0,m,n}^{1,-2,k}=-n_{0,m,n}^{-1,2,k}=n_{0,m,n}^{-1,-2,k},\hfill \\ \hfill n_{0,m,n}^{1,j,2}& =-n_{0,m,n}^{1,j,-2}=-n_{0,m,n}^{-1,j,2}=n_{0,m,n}^{-1,j,-2}.\hfill \end{array}$$(34)

Therefore, only four free parameters remain for determining all the coefficients for ϕ_0, m, n. These could be, for example $n_{0,m,n}^{1,2,\pm 1}$ and $n_{0,m,n}^{1,\pm 1,2}$. The same reasoning leads to similar relations for the coefficients $n_{l,0,n}^{i,j,k}$ and $n_{l,m,0}^{i,j,k}$, which can be obtained for the other internal fluxes after cyclic permutations.

Finally, to fulfill the third requirement we ask that the monopolar magnetic charge (the integral of the divergence) is split evenly between the eight fine cells. This can be expressed as

$$\begin{array}{cc}\hfill (\pm {\varphi }_{\pm 2,j,k}\mp {\varphi }_{0,j,k})& +(\pm {\varphi }_{i,\pm 2,k}\mp {\varphi }_{i,0,k})+(\pm {\varphi }_{i,j,\pm 2}\mp {\varphi }_{i,j,0})\hfill \\ \hfill & =\frac{1}{8}\mathrm{\Delta }V(\mathrm{\nabla }·\mathit{B}).\hfill \end{array}$$(35)

where ∇ ⋅ B follows from the discretization of Eq. (25). After substituting the expressions for the coarse fluxes in terms of the fine ones (Eq. (26)) and grouping the coefficients associated with each of the exterior fine fluxes, it is possible to use the independence of these fluxes to obtain nine equations for the twelve still unknown coefficients in each of the eight fine cells (more precisely, 24 equations are obtained for each fine cell, but 12 of them are redundant due to Eq. (34) and three are not independent from the others, leaving nine equations). Therefore, this constraint leaves three free parameters for each of the eight fine cells.

Since each of the internal faces is shared by two fine cells and its free parameters need to be compatible on both sides, it turns out that all of the coefficients of the internal fluxes can be written in terms of the same three parameters, which we name α₁, α₂ and α₃. For example, the interior magnetic fluxes in direction x¹ can be obtained in a divergence-free way as

$$\begin{array}{cc}& {\varphi }_{0-1-1}=\frac{1}{2}({\varphi }_{-2-1-1}+{\varphi }_{2-1-1})\hfill \\ \hfill & +\frac{1}{16}[(3+{\alpha }_2)F_{12-1}+(1-{\alpha }_2)F_{121}\hfill \\ \hfill & +(3-{\alpha }_3)F_{1-12}+(1+{\alpha }_3)F_{112}]\hfill \\ \hfill & {\varphi }_{0-11}=\frac{1}{2}({\varphi }_{-2-11}+{\varphi }_{2-11})\hfill \\ \hfill & +\frac{1}{16}[(1-{\alpha }_2)F_{12-1}+(3+{\alpha }_2)F_{121}\hfill \\ \hfill & +(3-{\alpha }_3)F_{1-12}+(1+{\alpha }_3)F_{112}]\hfill \\ \hfill & {\varphi }_{01-1}=\frac{1}{2}({\varphi }_{-21-1}+{\varphi }_{21-1})\hfill \\ \hfill & +\frac{1}{16}[(3+{\alpha }_2)F_{12-1}+(1-{\alpha }_2)F_{121}\hfill \\ \hfill & +(1+{\alpha }_3)F_{1-12}+(3-{\alpha }_3)F_{112}]\hfill \\ \hfill & {\varphi }_{011}=\frac{1}{2}({\varphi }_{-211}+{\varphi }_{211})\hfill \\ \hfill & +\frac{1}{16}[(1-{\alpha }_2)F_{12-1}+(3+{\alpha }_2)F_{121}\hfill \\ \hfill & +(1+{\alpha }_3)F_{1-12}+(3-{\alpha }_3)F_{112}],\hfill \end{array}$$(36)

where

$$\begin{array}{cc}& F_{12-1}={\varphi }_{12-1}-{\varphi }_{1-2-1}-{\varphi }_{12-1}+{\varphi }_{-1-2-1}\hfill \\ \hfill & F_{121}={\varphi }_{121}-{\varphi }_{1-21}-{\varphi }_{-121}+{\varphi }_{-1-21}\hfill \\ \hfill & F_{1-12}={\varphi }_{1-12}-{\varphi }_{1-1-2}-{\varphi }_{-1-12}+{\varphi }_{-1-1-2}\hfill \\ \hfill & F_{112}={\varphi }_{112}-{\varphi }_{11-2}-{\varphi }_{-112}+{\varphi }_{-11-2}.\hfill \end{array}$$(37)

Corresponding formulas for the other directions can be obtained by circular permutation of direction of the fluxes and of the indices of α_is. Since the parameters α_i affect only the interior fluxes, they may vary from cell to cell depending, for instance, on the cell geometry. For example, the formulas reduce to those by Tóth & Roe (2002) when

$$\begin{array}{c}\hfill {\alpha }_1=\frac{\mathrm{\Delta }{y}^2-\mathrm{\Delta }{z}^2}{\mathrm{\Delta }{y}^2+\mathrm{\Delta }{z}^2},{\alpha }_2=\frac{\mathrm{\Delta }{z}^2-\mathrm{\Delta }{x}^2}{\mathrm{\Delta }{z}^2+\mathrm{\Delta }{x}^2},{\alpha }_3=\frac{\mathrm{\Delta }{x}^2-\mathrm{\Delta }{y}^2}{\mathrm{\Delta }{x}^2+\mathrm{\Delta }{y}^2},\end{array}$$(38)

where we rename Δx¹ → x, Δx² → y and Δx³ → z to keep the notation simple. This choice has a directional bias which depends on the geometry of the cell, and which disappears when Δx = Δy = Δz and thus α_i = 0, which gives the form of the operator used by Cunningham et al. (2009). Another feature of this choice, and for any other in which α_i do not depend on the exterior magnetic fluxes, is that, prolongation can produce a z component in a field originally having only x and y components. We noticed that this can be avoided at the expense of making the operator nonlinear, by resorting to a heuristic formula for α_i to control the magnetic flux mixing. First, we quantify this mixing in the z direction as

$$\begin{array}{c}\hfill {\sigma }_z=\frac{{\sum }_{\mathrm{up}}|\varphi |-{\sum }_{\mathrm{down}}|\varphi |}{{\sum }_{\mathrm{up}}|\varphi |+{\sum }_{\mathrm{down}}|\varphi |},\end{array}$$(39)

where ∑_up runs on the exterior fluxes in the x and y direction residing in the upper half of the coarse cell, and ∑_down on those residing in the lower half. Corresponding expressions are defined for quantifying the mixing in the x and y directions. Then, we set the α_i-functions for interpolation as

$$\begin{array}{c}\hfill {\alpha }_1={\sigma }_y-{\sigma }_z,{\alpha }_2={\sigma }_z-{\sigma }_x,{\alpha }_3={\sigma }_x-{\sigma }_y.\end{array}$$(40)

By substituting, it is possible to verify that any magnetic-field configuration that originally had no components in a given coordinate direction will not show them when prolonged. A comparison between the interior fields interpolated using this nonlinear operator and that by Tóth & Roe is shown in Fig. 4.

Fig. 4. Comparison between the different projections of the interior fields obtained using Tóth & Roe’s operator (blue) and the nonlinear operator here derived (red), from the same exterior fluxes (black) and for a cubic cell. Each panel corresponds to a slice of a coarse cell, for which the indices of the fine cells are given at the upper left (cf. Fig. 3). The direction of the original flux is better preserved by the nonlinear operator (cf. right panels), where a component in the y-direction that was not originally present appears for the Tóth & Roe operator. The field interpolated using the nonlinear operator is exactly zero for cells with j = −1. For clarity, vectors in the lower left panel are represented with a double size.

4.3. Restriction of edge-allocated variables

On each block, magnetic fluxes are updated using the representation of ℰ correspondent to the level of the block. However, at a fine-coarse interface, two different representations of the magnetic flux exist. On the coarse side we have, for example, the magnetic flux Φ_2, 0, 0, while on the fine side, we have ∑_{i, j = −1, 1}ϕ_2, j, k. Due to the use of different representations of ℰ for the update, the two representations of the magnetic flux differ, although in absence of surface monopolar magnetic charge they should coincide. In order to make them consistent without creating divergence on the coarse side, we replace the coarse representation of ℰ used in the calculation of the circulation with the fine one. To do this, we store both representations of ℰ on the interface at the time they are calculated. After the time update, we subtract the line integrals of the electric field contained in the interface from the fluxes in contact with it. We then communicate the fine representations to the coarse side and construct a new coarse representation. For example, for an interface at constant x¹ between a coarse region at the left and a fine region at the right, we replace

$$\begin{array}{cc}& {\mathcal{E}}_{2,2,-1}+{\mathcal{E}}_{2,2,1}\to {\mathcal{E}}_{2,2,0}\hfill \\ \hfill & {\mathcal{E}}_{2,-1,2}+{\mathcal{E}}_{2,1,2}\to {\mathcal{E}}_{2,0,2}\hfill \\ \hfill & {\mathcal{E}}_{2,-2,-1}+{\mathcal{E}}_{2,-2,1}\to {\mathcal{E}}_{2,-2,0}\hfill \\ \hfill & {\mathcal{E}}_{2,-1,-2}+{\mathcal{E}}_{2,1,-2}\to {\mathcal{E}}_{2,0,-2}\hfill \end{array}$$

for each coarse cell and recalculate the circulation to update all the fluxes except for Φ_{−2, 0, 0}. In 2D, no sum is necessary and the coarse representations are simply replaced by the co-spatial fine ones.

4.4. Parallelisation and ghost-cell exchange

To simplify ghost-cell communications in the code, a block’s neighbors and its children are identified by local grid indices i1, i2, i3 = { − 1, 0, 1} and local child grid indices ic1, ic2, ic3 = {0, 1, 2, 3}, both from left (back, bottom) to right (front, top). For example, if block A is directly at the left of block B both belong to the same resolution level, block A is identified as neighbor ( − 1, 0) of block B, and block B is identified as neighbor (1, 0) of block A. This identification strategy is inherited from MPI-AMRVAC (see e.g., Fig. 4 of Keppens et al. 2012). When exchanging ghost cells, different array ranges are associated to each index combination, corresponding to the extent of the ghost zones that need to be filled from that neighbor. The identification of neighbors and ghost zones with integer indices allows to easily determine communication patterns using integer arithmetic operations. As an example, Fig. 5 shows some of these communication patterns for the case of a two-dimensional grid.

Fig. 5. Ghost-cell exchange for blocks (a) at the same resolution level and for coarse-fine interfaces in (b) restriction and (c) prolongation. The numbers represented are the local grid indices to identify neighbors and ghost regions, as explained in Sect. 4.4. Lines drawn in color denote (magnetic) fluxes on the boundary faces that are communicated. Panel c: dashed lines denote boundaries that contain only tangential fluxes necessary for restriction which reside only on the coarse buffer of a block. The size of fine blocks including ghost regions is marked by a gray line.

The three panels of Fig. 5 show examples of how these local indices are used to identify different ghost regions for (a) same resolution level communications, (b) restriction and (c) prolongation. To ensure that the communicated values are consistent, first same resolution communications are performed, then restrictions and finally prolongations. As in MPI-AMRVAC, to minimise the size of communications at resolution jumps, variables are restricted before being sent and are prolonged after being received. Each block possess a buffer containing a coarse representation of it to facilitate this operations.

In contrast to MPI-AMRVAC, for which inter-processor ghost-cell communications were based on MPI types created using the functions MPI_CREATE_SUBARRAY and MPI_TYPE_COMMIT (see Keppens et al. 2012), BHAC packs all communications, both for staggered and cell centered variables, in one-dimensional arrays using FORTRAN’s reshape function. These buffer arrays are transmitted to other processors using non-blocking MPI communications through the functions MPI_ISEND and MPI_IRECV. When received at the target destination, they are unpacked to fit the shape of the ghost region that they need to fill.

A reason for this change is that the array segments that need to be communicated for each component of the staggered field are different among them and different from those of the cell-centered variables; therefore, using a single buffer instead of a different MPI_TYPE for each of them reduces the number of necessary communications. These and other changes done in the ghost-cell exchange routines are oriented to a future upgrade of the code to a task-based scheduling, which would result in a significant speedup and facilitate the use of hierarchical time-stepping.

Finally, fluxes and electric fields necessary for the consistency steps described in Sect. 4.3 are also written in special communication buffers for sending and receiving. These buffers are likewise one-dimensional arrays filled and extracted using the Fortran reshape function.

4.5. Poles in spherical and cylindrical coordinates

In three-dimensional simulations in spherical and cylindrical geometry, the ghost-cell exchange at the poles requires special attention. Despite being located at the minimum (and maximum) θ for spherical and ρ for cylindrical coordinates, in BHAC the pole is not treated as a physical boundary nor it is excised as is usually done in many codes (see e.g., Shiokawa et al. 2012). Instead, the connectivity of the blocks touching the poles is changed so that each on them considers as its neighbors the blocks at the opposite side of the polar axis and reconstructions can be performed across the latter.

The only difference respect to the usual ghost-cell exchange is that, after being received, the elements of the receive buffer are reversed in the ϕ direction, as is shown in Fig. 6. Likewise, periodic boundary conditions such as those in the ϕ direction, are handled by changing the connectivity of the grid such that the first block in the ϕ direction becomes the neighbor of the last one.

Fig. 6. Ghost-cell exchange at the south and north poles in spherical coordinates for same resolution level (a), restriction (b) and prolongation (c). For clarity, the pole is represented as an expanded black circumference to identify the cell faces touching it, although no there is no excised region. As in Fig. 5, colored lines denote (magnetic) fluxes on the boundary faces that are communicated. For prolongation, dashed lines denote boundaries that contain only tangential fluxes necessary for restriction and the total extent of fine blocks including ghost regions is marked by a gray line. For cylindrical coordinates, the ghost-cell exchange at the z pole happens as for the north pole.

The flux-fixing operation is not necessary at poles, since the areas of cell surfaces touching the pole are exactly zero and no flux should be present. Similarly, the contribution of the circulation of the electric field is exactly zero along the polar axis, by recalling that in a coordinate basis $\sqrt{\gamma }$ vanishes there (see Eqs. (20) and (22)). Therefore, no fixings are performed and both numerical fluxes and electric fields are set to zero at the pole.

4.6. Boundary conditions

When non-periodic boundary conditions are used, it is necessary to ensure that no divergence is created at boundaries. To accomplish this, we first extrapolate the components of the magnetic field tangential to the boundary according to a user-given prescription, for example, symmetric, antisymmetric, flat, etc. Then, the normal component is filled layer by layer outward from the boundary, cancelling the sum of the magnetic fluxes for each cell. For symmetric and antisymmetric boundary conditions, this procedure can be applied without modifications also at resolution jumps; however, this is not the case for flat boundary conditions. The reason is that they consist in copying the value of a variable at the last cell of the physical domain to all of the ghost cells; and in a resolution jump, the sum of two times the last magnetic flux on the fine side is not necessarily equal to the last magnetic flux on the coarse side. Therefore, instead of filling the ghost cells with the values of the last row in the domain, we fill them with an average of the last two rows, weighted by the normal surfaces. In this way, the fine and the coarse representations of the fields at the ghost cells are consistent.

5. Numerical tests

In this section we describe the results of three well known numerical tests performed to validate the new schemes, two in special relativity and one in general relativity. These are the Gardiner & Stone advected magnetic loop (Gardiner & Stone 2005), the cylindrical explosion of Komissarov (1999), and the highly magnetized stationary torus by Komissarov (2006). Since the new additions (CT) and AMR can only be tested in 2D and 3D, and the code has already been validated for 1D problems (see Porth et al. 2017), we omit 1D tests. We also omit tests using both the cell-centered and the staggered variants of arithmetic averaging, except for comparisons with UCT, as the former have already been published in Porth et al. (2017) and Olivares Sánchez et al. (2018b).

We are therefore interested in testing the UCT algorithms and the new AMR features, as well as the recently implemented limiters, WENO Z+ and MP5. All of the tests presented here use the equation of state of an ideal fluid with $\widehat{\gamma }=4/3$. The prescription for triggering AMR is the Löhner scheme (Löhner 1987), which decides to coarsen or refine based on a quantification of oscillations in the numerical solution.

5.1. Loop advection

A well known problem to assess the importance of spurious effects due to the divergence-control technique employed in an MHD code is the advection of a weakly magnetized loop. This test was originally used by Gardiner & Stone (2005) to study a divergence-preserving scheme with the upwind property. In this section, we present a comparison of the results obtained using the cell-centered version of flux-CT and the upwind method UCT2.

On a uniform background with ρ = 1, p = 1, ν^x = 0.2, and ν^y = 0.1, the problem consists on advecting a magnetic loop with radius R = 0.3, described by the potential A_z = A₀(R − r) for r ≤ R and A_z = 0 for r >  R, where r² = x² + y². A₀ is chosen 10⁻³, so that β = p_fluid/p_mag = 10⁶ and the loop is advected passively. The domain is x ∈ [ − 1, 1] and y ∈ [ − 0.5, 0.5] with periodic boundary conditions. The methods employed are the MP5 limiter, HLLE fluxes and the RK3 time-integrator.

Two realizations of the same set-up are evolved: one using the flux-CT method, and the other UCT2. While for the former it is not possible to use AMR, for the latter three levels of AMR are used. The highest AMR-level has a resolution equivalent to that of the flux-CT run, i.e, 256 × 128.

The magnetic pressure at the initial condition and some snapshots of the evolution are depicted in Fig. 7. The upper-left panel of the same figure shows the divergence of the magnetic field for the UCT2 run at the final time t = 10.0, when the magnetic loop, and thus the refined region, has travelled a complete cycle across the domain. Divergence remains always zero at machine precision (∼10⁻¹⁶) despite resolution jumps, and the coarsening and refining of blocks. It can be seen that, in addition to the AMR capacity, the staggered UCT scheme is able to preserve much better the original shape of the loop, without creating as many spurious oscillations as flux-CT.

Fig. 7. Gardiner & Stone’s loop advection test. Top left: divergence of the magnetic field after ten time units, using UCT2. Bottom left: magnetic pressure at the initial configuration. Top right: after ten time units, using UCT2, showing the AMR blocks of 8 × 8 cells. Bottom right: after ten time units, using flux-CT.

5.2. Cylindrical explosion

Another challenging problem for MHD codes is that of a cylindrical blast wave expanding in a homogeneous plasma with an initially constant magnetic field. We present results of such problem, which highlights the importance of upwinding in the calculation of electric fields, and show the different behavior of the BS and the UCT2 schemes. The problem set-up is the same that appears in Del Zanna et al. (2007) and Komissarov (1999).

The domain extends over [ − 6, 6]×[−6, 6]. Plasma is initially at rest, with a constant magnetic field B^x = 0.1, B^y = B^z = 0. Defining r² = x² + y², the domain is divided in three regions: a high density and pressure (ρ_in = 10⁻², p_in  =  1) region for r <  r_in  =  0.8, a low density and pressure (ρ_out  =  10⁻⁴, p_out  =  5 × 10⁻⁴) region for r >  r_out  =  1.0, and a transition region for r_in ≤ r ≤ r_out where ρ = ρ_in(r/r_in)^α₁ and p = p_in(r/r_in)^α₂, and the exponents α₁ and α₂ are such that p(r) and ρ(r) are continuous at r_in and r_out. The boundary conditions are periodic, although this is not relevant in practice, since the region of interest cannot be affected by signals coming form the boundaries during the simulation time.

The system is evolved until t = 4.0 with an RK3 integrator and a CFL factor of 0.35, and HLLE fluxes with WENOZ+ reconstruction are used. The simulation used three AMR levels, with a base resolution of 100 × 100, giving an effective resolution of 400 × 400. To compare the two algorithms for magnetic-field evolution, we run two simulations, one using the BS algorithm and the other using UCT2.

Figure 8 shows the last snapshot of the evolution for both methods. Spurious oscillations behind the reverse shock are more pronounced for the BS method, similar to those already encountered in Sect. 5.1.

Fig. 8. Cylindrical explosion. y-component of the magnetic field at t = 4.0 when using the BS algorithm (left) and UCT2 (right).

The situation greatly improves when adopting the UCT2 algorithm, as can be further confirmed by examining the profiles of the y-component of the magnetic field along the line x = 3 shown in Fig. 9.

Fig. 9. Cylindrical explosion. Profile of the y-component of the magnetic field along a line at x = 3. The Balsara & Spicer algorithm shows a more oscillatory behavior as compared with UCT2.

5.3. Stationary torus

The test presented in this section is the analytic solution found by Komissarov (2006) for an equilibrium torus with a toroidal magnetic field. As in Porth et al. (2017), the implementation of the initial condition in BHAC is based on routines kindly provided by Chris Fragile. Taking advantage of the stationarity of this analytic solution, we use this test to quantify the convergence of the numerical solution. Recalling the results from Komissarov (2006), for a fluid of constant angular momentum distribution l = −u_ϕ/u_t = l₀, and adopting a specific relation of the fluid and magnetic pressure with respect to the fluid enthalpy, the configuration is described by the equation

$$\begin{array}{c}\hfill \mathcal{W}-{\mathcal{W}}_{\mathrm{in}}=-\frac{\kappa }{\kappa -1}\frac{p}{\rho h}-\frac{\eta }{\eta -1}\frac{p_m}{\rho h},\end{array}$$(41)

where 𝒲 = ln|u_t| and κ and η are constants. To completely specify the solution, it is necessary to give two more parameters, namely the fluid enthalpy and the ratio of fluid to magnetic pressure at the center of the torus, (ρh)_c = ω_c and (p/p_mag)_c = β_c, respectively. The parameters of the torus chosen here as a test problem are κ = 4/3, η = 4/3, l₀ = 2.8, r_c = 4.62, 𝒲_in = −0.030, ω_c = 1.0, and β_c = 0.1. The dimensionless spin of parameter of the black hole is set to a = 0.9, and the plasma follows the equation of state of an ideal fluid with $\widehat{\gamma }=4/3$. Simulations are performed in 3D and using logarithmic Kerr-Schild coordinates (i.e., MKS coordinates with ϑ₀ = 0 and R₀ = 0). The simulation domain spans over θ ∈ [π/5, 4π/5], r_KS ∈ [0.95 r_h,  50 M] and ϕ ∈ [0, π], where r_KS is the radial Kerr-Schild coordinate and r_h is the (outer) event horizon radius of the BH. To test convergence, we performed three simulations progressively doubling the resolution from a base grid of N_r × N_θ × N_ϕ = 100 × 60 × 40. The rotation period of the center of the torus is 68 M, and the simulation is carried out until t = 100 M, thus spanning nearly one and a half orbital periods. Although this kind of solutions are known to be unstable to the MRI and the Papaloizou-Pringle instability in 3D (Wielgus et al. 2015; Bugli et al. 2018), the time elapsed by our simulations still corresponds to the initial phase of the linear growth, and no chaotic behavior can be observed. In fact, the torus here presented is identical to Case A simulated by Wielgus et al. (2015), who reported perturbations at t = 100 M still at the ∼10⁻⁴ level.

To emulate vacuum regions outside the torus, a tenuous atmosphere with density and pressure given by $\rho ={\rho }_{\min }{r}_{\mathrm{KS}}^{-3/2}$, $p=p_{\min }{r}_{\mathrm{KS}}^{-5/2}$ with ρ_min = 10⁻⁵ and p_min = 10⁻⁷ was added. Initially, the velocity of the atmosphere is set to zero in the Eulerian frame. However, the atmosphere is free to evolve, and only density and pressure are reset to the atmosphere value whenever they fall below it. The numerical methods chosen for this test were the HLLE Riemann solver, LIMO3 reconstruction and two-step time integration with CFL factor 0.5.

Figure 10 shows the density distribution, the plasma beta and the divergence of the magnetic field at the initial state and at t = 100 M for the run with resolution N_r = 400. As can be seen, the distribution of the displayed quantities is very well preserved by the scheme. The L₁ and L_∞ errors (calculated by comparing with the last state with the initial condition) are shown in Fig. 11. Consistently with the algorithm employed, second order convergence can be observed.

Fig. 10. Rest mass density, plasma beta and divergence of the magnetic field at t = 0 and t = 100 for the run with Nr = 400 of the Komisssarov torus test.

Fig. 11. L∞ and L1 norms of the error in density ρ for the strongly-magnetised Komissarov torus at t = 70 M, i.e., after one orbital period. The scheme shows second order convergence.

6. Magnetized black hole accretion

In this section we present simulations where the newly implemented numerical methods are applied to evolve the accretion of plasma onto black holes, one of the main target applications of the code. In Porth et al. (2017), simulations of this kind, performed at uniform resolution, were used to validate the code by means of qualitative and quantitative comparisons with the widely used GRMHD code HARM (Gammie et al. 2003; Noble et al. 2009). Therefore, this section will be focused in showing the advantages of AMR and the effects of the choice of divergence control technique.

For all the simulations in this section, the initial configuration is a torus in equilibrium around a black hole (Fishbone & Moncrief 1976) with spin parameter a = 0.9375, and thus with the outer event horizon located at r = 1.348 M, where r is the radial Kerr-Schild coordinate. The inner radius of the torus is at r_in = 6 M, and its density maximum at r_max = 12 M (orbital period of 247 M at the density maximum). The density in the torus is normalized so that it takes 1.0 as its a maximum value. A single-loop poloidal magnetic field is built from the vector potential A_ϕ ∝ max(ρ/ρ_max − 0.2, 0) and is normalized in such a way that the highest magnetic pressure and the highest thermal pressure are related by the ratio β = p_fluid, max/p_mag, max = 100. To break this equilibrium state, random perturbations of 4% are added to the pressure. This eventually triggers the MRI, allowing the plasma to accrete. It is worth mentioning that even without explicitly adding a perturbation, numerical errors produced by the discretization would be amplified to produce a turbulent state very similar to that obtained. However, in order to have some control on the initial state, the perturbation we add is significantly larger than numerical errors (which are not easily predictable) but still small respect to the equilibrium pressure. In addition, since the initial growth rate of the perturbation depends on its amplitude, the saturated state can be reached faster with larger perturbations. For this reason, starting with seed perturbations is computationally less expensive than waiting for the initial discretization errors in the equilibrium state to grow. Moreover, the value of 4% is the same used in Porth et al. (2017), and is therefore useful for making comparisons.

To avoid vacuum regions, the rest of the simulation is filled with a tenuous atmosphere with density ρ_fl = 10⁻⁴ r^−3/2 and fluid pressure p_fl = 1/3 × 10⁻⁶ r^−5/2. We reset the density or the pressure whenever they fall below these floor values. Simulations were evolved using a two-step time-integrator.

As mentioned in Sect. 2.2, several coordinate systems for black hole spacetimes are available in BHAC. The simulations shown here are performed in MKS, as well as in CKS coordinates. Both coordinate systems possess advantages and weaknesses that must be taken into account depending on the aspect of the physical system under study, as will be explained below. The simulations whose results are reported in this section are summarized in Table 2.

6.1. 2D MKS simulations

The first pair of simulations is performed in 2D on the meridional plane. The domain covers θ ∈ [0, π] and r ∈ [1.2, 2500], and is resolved using three AMR levels triggered by the Löhner scheme, with a base resolution of N_r  ×  N_θ = 200  ×  100 and an effective resolution of N_r  ×  N_θ = 800  ×  400. Since the time step for these 2D simulations is not penalised by the small width of cells in contact with the polar axis, the parameter ϑ₀ of MKS coordinates is set to zero. Reconstruction is performed using PPM. The purpose of these two simulations is to investigate the effect of the method utilised for the evolution of the magnetic field; therefore, the only difference between them was the use of either the Balsara & Spicer algorithm or UCT2.

A comparison at t  =  5000 (Fig. 12) displays what we observe to be systematic differences between the two simulations. Namely, a higher magnetization inside the funnel for the UCT2 run and a different morphology of the jet, which acquires a parabolic shape for the Balsara & Spicer run and is more conical for the UCT2 run. Another difference lays in the transition between the funnel and the mildly magnetized wind surrounding the jet, which is sharper for the UCT2 case. These differences could have important observational consequences, for instance, in future high-resolution VLBI images of the SMBH candidates in M 87 and Sgr A*, since radiative models of these sources which are able to reproduce the properties of their radio spectrum consider that most of the radio emission originates at the funnel wall (Mościbrodzka et al. 2016).

Fig. 12. Comparison between 2D simulations of black hole accretion using the Balsara & Spicer method (left) and the UCT2 algorithm (right) at t = 5000 M. For each panel, rest mass density ρ is shown at the left and plasma β at the right.

From these simulations, it is possible to see that the choice of method for evolving the magnetic field can have a visible effect on the flow behavior. Among these two methods, we would recommend the use of UCT2 due to its upwind properties and to the possibility of using high-order reconstructions. Although some asymmetry respect to the equator can be seen in both the BS and the UCT2 simulations, this is only a consequence of the chaotic evolution of the initial random perturbations driven by turbulence, and not of the numerical schemes. Even though not included in the tests presented in this work, in order to verify that the scheme does not possess any directional bias, we have simulated the development of the magnetic Kelvin-Helmholtz instability with a controlled perturbation (Tóth & Roe 2002) and employing the same methods (PPM and BS or UCT2). When comparing with a simulation starting from an initial condition obtained by specular reflection, we obtained identical results. Finally, it is worth mentioning that the simulation using the Balsara & Spicer method was tested previously against another run performed in an equivalent uniform resolution and using the widespread flux-CT method. The AMR run was able to achieve a speedup factor of 7.1, while obtaining quantitatively similar results in the mass accretion rate and the magnetic flux through the horizon, as discussed in Olivares Sánchez et al. (2018b).

6.2. 3D MKS simulations

Since self-sustaining dynamo activity leading to the perpetuation of the MRI cannot occur in strict azimuthal symmetry (Cowling 1933; Balbus & Hawley 1991), 3D simulations are necessary to study the accretion flow in the saturated state. The AMR capabilities of the code become even more important in this case, where the computational cost of simulations rapidly increases due to the larger number of computing cells. In addition, the polar coordinate system naturally leads to a significantly higher resolution close to the polar axis, which is not always justified, and the narrow cells in these regions usually cause a penalization in time-step size.

A possibility to alleviate this limitation is to de-refine cells close to the polar axis, as done by White et al. (2016). In typical accretion scenarios, the polar region is occupied by an evacuated, smooth funnel, thus no especially high resolution is required. In contrast, equatorial regions are populated by turbulent structures that need to be properly resolved.

Newtonian shearing box simulations have shown that an insufficient resolution can suppress the growth of the MRI when the wavelength of its fastest growing mode λ_MRI is resolved with less than six cells, and move the simulation away from the ratio of magnetic to gas pressure expected at the saturated state (Sano et al. 2004). A relativistic version of this criterion has been applied to estimate whether global accretion simulations are properly resolved (see e.g., Noble et al. 2010; McKinney et al. 2012).

The MRI quality factor is defined as the ratio between the grid spacing and λ_MRI, which approximately gives the number of cells employed to resolve it. Following Noble et al. (2010), we define the relativistic MRI quality factor as

$$\begin{array}{c}\hfill Q_{\mathrm{MRI}}^{(i)}=\frac{{\lambda }_{\mathrm{MRI}}^{(i)}}{\mathrm{\Delta }{x}^{(i)}}=\frac{2\pi }{\mathrm{\Omega }}\frac{|b^{(i)}|}{\mathrm{\Delta }{x}^{(i)}\sqrt{\rho h+b^2}},\end{array}$$(42)