Issue 
A&A
Volume 629, September 2019



Article Number  A61  
Number of page(s)  21  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201935559  
Published online  05 September 2019 
Constrained transport and adaptive mesh refinement in the Black Hole Accretion Code
^{1}
Institute for Theoretical Physics, GoetheUniversity, 60438 Frankfurt am Main, Germany
email: olivares@th.physik.unifrankfurt.de
^{2}
Astronomical Institute Anton Pannekoek, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
^{3}
Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{4}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{5}
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
Received:
27
March
2019
Accepted:
16
June
2019
Context. Worldwide very long baseline radio interferometry (VLBI) arrays are expected to obtain horizonscale images of supermassive black hole candidates and of relativistic jets in several nearby active galactic nuclei. This, together with the expected detection of electromagnetic counterparts of gravitationalwave signals, motivates the development of models for magnetohydrodynamic flows in strong gravitational fields.
Aims. The Black Hole Accretion Code (BHAC) is a publicliy available code intended to aid with the modeling of such sources by performing general relativistic magnetohydrodynamical simulations in arbitrary stationary spacetimes. New additions to the code are required in order to guarantee an accurate evolution of the magnetic field when small and large scales are captured simultaneously.
Methods. We discuss the adaptive mesh refinement (AMR) techniques employed in BHAC, which are essential to keep several problems computationally tractable, as well as staggeredmeshbased constrained transport (CT) algorithms to preserve the divergencefree constraint of the magnetic field. We also present a general class of prolongation operators for faceallocated variables compatible with them.
Results. After presenting several standard tests for the new implementation, we show that the choice of the divergencecontrol method can produce qualitative differences in the simulation results for scientifically relevant accretion problems. We demonstrate the ability of AMR to decrease the computational costs of black hole accretion simulations while sufficiently resolving turbulence arising from the magnetorotational instability. In particular, we describe a simulation of an accreting Kerr black hole in Cartesian coordinates using AMR to follow the propagation of a relativistic jet while selfconsistently including the jet engine, a problem set up for which the new AMR implementation is particularly advantageous.
Conclusions. The CT methods and AMR strategies discussed here are currently being used in the simulations performed with BHAC for the generation of theoretical models for the Event Horizon Telescope collaboration.
Key words: magnetohydrodynamics (MHD) / relativistic processes / methods: numerical / accretion / accretion disks / black hole physics
© ESO 2019
1. Introduction
The prospect of horizonscale 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 nonthermal 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 highresolution images of other sources of great interest. For instance, observations of the twosided 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 largescale features such as the position of recollimation 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 recollimating 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^{1}), 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. Divergencecleaning 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 nonrelativistic magnetohydrodynamics (Toth 2000; Balsara & Kim 2004; Mocz et al. 2016) show that all tested variants of the divergencecleaning method produce important spurious oscillations in the magneticfield energy. Balsara & Kim (2004) have attributed these oscillations to the nonlocality introduced by the constraintdamping equation, thus suggesting that all divergencecleaning 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 divergencecleaning technique known as Generalized Lagrange Multipliers (GLM, Dedner et al. 2002) and the method known as fluxinterpolated Constrained Transport (fluxCT, Toth 2000).
By adopting a discretization of Faraday’s law consistent with that of the constraint, ConstrainedTransport (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 divergencefree and care must be taken in order not to generate magneticfield divergence by any other means. This can happen, in particular, at boundary cells, at coarsefine 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 zerospeed modes (Etienne et al. 2012).
Furthermore, even though cellcentered 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 divergencepreserving prolongation and restriction operation for cellcentered magnetic fields is underdetermined, that is, information is lost when abandoning the face representation in favor of the cellcentered representation (Tóth & Roe 2002).
The approach used in this work consists in applying divergencepreserving restriction and prolongation operators for faceallocated 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), notupwind algorithms such as the widely used fluxCT and the Balsara & Spicer (1999) method (BS) are prone to numerical instabilities. Upwind CT schemes include the secondorder 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 highorder schemes for the integration of Faraday’s law.
This work largely focuses on the AMRcompatible 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 divergencepreserving prolongation operators, which constitute a coordinateindependent 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 ghostcell exchange, the data structures and the special treatments employed at coarsefine 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 massaccretion 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 finitevolume (Sect. 3.1) and CT methods (Sect. 3.2). The mathematical and technical aspects of BHAC’s AMRcompatible 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 MPIAMRVAC framework toward GRMHD simulations in one, two, and three dimensions using finitevolume methods and a variety of modern numerical methods, described more in detail in Porth et al. (2017). It exploits much of MPIAMRVAC’s infrastructure for parallelization and blockbased 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 restmass density, the energy and momentum, and the homogeneous Maxwell equations
where ρ is the particle number density in the fluid frame, u^{μ} is the fluid 4velocity in the observer frame, T^{μν} is the energymomentum 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 4velocity U^{μ} are
and are related as where ϵ^{μναβ} are the components of the LeviCivita tensor. Therefore, in a general frame moving with 4velocity U^{μ}
and
Instead of solving the inhomogeneous Maxwell equations, we close the system by enforcing the idealMHD condition
where e^{μ} denotes the electric field in the fluid frame. Consequently, the magnetic field in the fluid frame is b^{μ}=^{*}F^{μν}u_{ν} and
Physically, the idealMHD 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 (frozenflux theorem).
The energymomentum 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
where h_{tot} := h + b^{2}/ρ is the total specific enthalpy and p_{tot} = p + b^{2}/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 temperaturedependent 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 3metric induced in the hypersurface.
The explicit form of the line element in a 3+1 split spacetime is
where α and β^{i} 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 3momentum S_{i} := n_{μ}γ_{νi}T^{μν}, the total energy U := n_{μ}n_{ν}T^{μν} and the spatial stress tensor W^{ij} := γ_{μi}γ_{νj}T^{μν}.
In terms of these variables, system (1) can be written as a set of conservation equations
in addition to the solenoidal constraint for the magnetic field . Here, γ is the determinant of the induced threemetric, and the vectors of conserved quantities U, the fluxes F^{i}, and the sources S are given by
where is the 3velocity, Γ := −u^{μ}n_{μ} is the Lorentz factor^{2}, and 𝒱^{i} := αν^{i} − β^{i} 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^{i}, knowledge of the primitive variables P = [ρ,Γν^{i},p,B^{i}] 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 ξ := Γ^{2}ρ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 horizonscale VLBI images of Sgr A*, as well as the study of quasiperiodic 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 complexstep 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 4metric 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).
Coordinates available in BHAC.
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 KerrSchild 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 KerrSchild coordinates x^{μ′} = (t_{KS}, r_{KS}, θ_{KS}, ϕ_{KS}) by the transformation
where R_{0} and ϑ_{0} are constant parameters. The convention used here for the parameter ϑ_{0} is different from that in McKinney & Gammie (2004). The one employed here is chosen so that ϑ reduces to θ_{KS} when ϑ_{0} = 0, thus retrieving logarithmic KerrSchild 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. Finitevolume scheme
To introduce the notation used in this work, we next summarize the finitevolume 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 of the conserved quantities inside each cell
Quantities indicated as, say, F^{1}ΔS^{1}_{i + 1/2, j, k} represent integrals of the fluxes over the surfaces ΔS^{1}_{i + 1/2, j, k} bounding the control volume, and is the volume average of the sources. Both kinds of integrals are approximated to second order, by assigning to F^{n} (n = 1, 2, 3) the point value of the flux at the face center and to the point value at the cell barycenter, that is, at
respectively, where Δx^{i} is the grid spacing in each direction and is Kronecker’s delta. F^{n} 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 fourthorder Simpson’s rule and stored in memory.
The Riemann solvers currently available in BHAC are the Rusanov method, also known as total variation diminishing LaxFriedrichs 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 secondorder, total variation diminishing (e.g., minmod, vanLeer), and highorder such as PPM (Colella & Woodward 1984), LIMO3 (Čada & Torrilhon 2009) and MP5 (Suresh & Huynh 1997). Recent additions to that list are the highorder schemes WENO5 and WENOZ+ (Acker et al. 2016). Equation (12) can be integrated using any of the methods present in the MPIAMRVAC toolkit. These include the simple halfstep modified Euler, the third order RungeKutta RK3 (Gottlieb & Shu 1998) and the strongstability preserving sstep, pthorder 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 divergencecontrol 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 magneticfield 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
is equivalent to ∇ ⋅ B = 0 in the continuous limit. When the limit is not taken, that is, for the finitevolume case, it follows from the divergence theorem that the average value of the divergence 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
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 F^{i} (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
The magnetic flux at each face is updated from the circulation of the electric field, using the integral form of Faraday’s law
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 , 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 is kept constant at machine precision from one iteration to the next.
In order for B to be divergencefree during the whole simulation, 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 divergencefree 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 Fluxinterpolated Constrained Transport (fluxCT), although in most of the literature this name is used exclusively for the cellcentered version of this method, proposed also in Toth (2000). In this work, we will use the abbreviation fluxCT to refer to the cellcentered version of the method and BS to refer to the staggered one. Both variants are particularly suitable for finitevolume 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 magneticfield components are surface integrals of the electric field, for example, the flux in the x^{2}direction for B^{1} is
The innermost integral is the same as that of Eq. (15), so the average flux can be interpreted as
where is the mean value of the integral from Eq. (15) over the face at j + 1/2. To secondorder accuracy, this integral takes the value 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
Although this algorithm, especially in its cellcentered 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 onedimensional 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 magneticfield continuity and transport properties by using limited reconstructions. Also, in contrast to arithmetic averaging, UCT reduces to the correct 1dimensional 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^{1}, x^{2}, x^{3}) 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) are computed using the quantities in the boxes. 

At the time of calculating the numerical fluxes at xinterfaces, store the characteristic speeds and (i = x, y, z) obtained from the Riemann solver, as well as the transport velocity , , , and .

Reconstruct transport velocities along direction y toward the edges z and obtain one correspondent to the four states around them , , and (i = x, y).

Reconstruct the magnetic fields on the faces to the same edges and obtain , , and .

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

Take the maximum characteristic speeds in each direction from the four states,
where are defined as positive in the −i (+i) direction, and they are zero otherwise.

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):
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:

At the time of calculating the numerical fluxes at each interface, store the characteristic speeds and (i = x, y, z) obtained from the Riemann solver, and the transverse transport velocities at the left and right states, for example , , , and for the xinterface.

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 xinterface, this is

Reconstruct the magnetic fields and the transport velocities to the edge where ℰ_{z} needs to be calculated. This gives and .

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

Approximate the electric field using the formula by (Del Zanna et al. 2007):
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 and not the magnetic fields B^{i} 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. CTadapted 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 finecoarse interfaces, which is done in a step similar to the refluxing of finitevolume 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.
The coarse coordinate increment in the direction i is denoted ΔX^{i} and the correspondent fine increment Δx^{i} = ΔX^{i}/2. The center of the coarse cell is defined at and those of the eight fine cells that result after refinement at . 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^{1} are centered at . 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
and that across ΔS_{2, 0, 0} is
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 and ℰ_{2, 2, 1} that along the fine edge centered at .
In a 2D problem, where symmetry in the x^{3} direction is assumed, the integrals with respect to x^{3} 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
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 faceallocated 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,
The divergence in the newly created coarse cell will be
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 divergencefree 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 vertexcollocated fields. Since the divergence constraint alone results in an underdetermined 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 divergencepreserving 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 secondorder, totalvariationdiminishing (TVD), in addition to have vanishing divergence. The formulas by Tóth & Roe (2002) for facecollocated 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 nonrelativistic 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 3dimensional hypersurface, the expressions for the divergence (∇ ⋅ B = γ^{−1/2}∂_{i}γ^{1/2}B^{i}) and the curl ((∇×B)^{i} = ϵ^{ijk}∂_{j}B_{k}) of a 3vector 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 nonorthonormal bases, this requires to compute additional approximations to derivatives or interpolations of the magneticfield 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 faceallocated 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 facecollocated 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
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
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}.
The first requisite can be expressed as
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_{0}, but for the same cell must not change between the face at r = r_{0} and that at r = r_{0} + Δr. Applied to the fluxes in the first direction, requirement Eq. (31) results in 12 equations for the coefficients of each flux ϕ_{lmn}, namely
where 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,
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
Therefore, only four free parameters remain for determining all the coefficients for ϕ_{0, m, n}. These could be, for example and . The same reasoning leads to similar relations for the coefficients and , 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
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 α_{1}, α_{2} and α_{3}. For example, the interior magnetic fluxes in direction x^{1} can be obtained in a divergencefree way as
where
Corresponding formulas for the other directions can be obtained by circular permutation of direction of the fluxes and of the indices of α_{i}s. 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
where we rename Δx^{1} → x, Δx^{2} → y and Δx^{3} → 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
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
By substituting, it is possible to verify that any magneticfield 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 ydirection 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 edgeallocated variables
On each block, magnetic fluxes are updated using the representation of ℰ correspondent to the level of the block. However, at a finecoarse 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^{1} between a coarse region at the left and a fine region at the right, we replace
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 cospatial fine ones.
4.4. Parallelisation and ghostcell exchange
To simplify ghostcell 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 MPIAMRVAC (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 twodimensional grid.
Fig. 5. Ghostcell exchange for blocks (a) at the same resolution level and for coarsefine 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 MPIAMRVAC, 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 MPIAMRVAC, for which interprocessor ghostcell 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 onedimensional arrays using FORTRAN’s reshape function. These buffer arrays are transmitted to other processors using nonblocking 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 cellcentered 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 ghostcell exchange routines are oriented to a future upgrade of the code to a taskbased scheduling, which would result in a significant speedup and facilitate the use of hierarchical timestepping.
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 onedimensional arrays filled and extracted using the Fortran reshape function.
4.5. Poles in spherical and cylindrical coordinates
In threedimensional simulations in spherical and cylindrical geometry, the ghostcell 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 ghostcell 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. Ghostcell 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 ghostcell exchange at the z pole happens as for the north pole. 
The fluxfixing 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 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 nonperiodic 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 usergiven 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 cellcentered 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 . 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 divergencecontrol 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 divergencepreserving scheme with the upwind property. In this section, we present a comparison of the results obtained using the cellcentered version of fluxCT 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_{0}(R − r) for r ≤ R and A_{z} = 0 for r > R, where r^{2} = x^{2} + y^{2}. A_{0} is chosen 10^{−3}, so that β = p_{fluid}/p_{mag} = 10^{6} 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 timeintegrator.
Two realizations of the same setup are evolved: one using the fluxCT 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 AMRlevel has a resolution equivalent to that of the fluxCT run, i.e, 256 × 128.
The magnetic pressure at the initial condition and some snapshots of the evolution are depicted in Fig. 7. The upperleft 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^{−16}) 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 fluxCT.
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 fluxCT. 
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 setup 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^{2} = x^{2} + y^{2}, the domain is divided in three regions: a high density and pressure (ρ_{in} = 10^{−2}, p_{in} = 1) region for r < r_{in} = 0.8, a low density and pressure (ρ_{out} = 10^{−4}, p_{out} = 5 × 10^{−4}) region for r > r_{out} = 1.0, and a transition region for r_{in} ≤ r ≤ r_{out} where ρ = ρ_{in}(r/r_{in})^{α1} and p = p_{in}(r/r_{in})^{α2}, and the exponents α_{1} and α_{2} 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 magneticfield 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. ycomponent 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 ycomponent of the magnetic field along the line x = 3 shown in Fig. 9.
Fig. 9. Cylindrical explosion. Profile of the ycomponent 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_{0}, and adopting a specific relation of the fluid and magnetic pressure with respect to the fluid enthalpy, the configuration is described by the equation
where 𝒲 = lnu_{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_{0} = 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 . Simulations are performed in 3D and using logarithmic KerrSchild coordinates (i.e., MKS coordinates with ϑ_{0} = 0 and R_{0} = 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 KerrSchild 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 PapaloizouPringle 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^{−4} level.
To emulate vacuum regions outside the torus, a tenuous atmosphere with density and pressure given by , with ρ_{min} = 10^{−5} and p_{min} = 10^{−7} 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 twostep 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_{1} 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 N_{r} = 400 of the Komisssarov torus test. 
Fig. 11. L_{∞} and L_{1} norms of the error in density ρ for the stronglymagnetised 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 KerrSchild 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 singleloop 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^{−4} r^{−3/2} and fluid pressure p_{fl} = 1/3 × 10^{−6} r^{−5/2}. We reset the density or the pressure whenever they fall below these floor values. Simulations were evolved using a twostep timeintegrator.
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 ϑ_{0} 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 highresolution 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 highorder 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 KelvinHelmholtz 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 fluxCT 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 selfsustaining 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 timestep size.
A possibility to alleviate this limitation is to derefine 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
where b^{(i)} and Δx^{(i)} are, respectively, the magnetic field and the displacement in the ith direction in the orthonormal frame comoving with the fluid, and Ω is the orbital angular velocity. Since the value of b^{(i)} depends on the amplification of the magnetic field experienced during the simulation, can only be determined a posteriori. For a fluid rotating in the ϕ direction, is the relevant MRI quality factor.
The simulation presented in this section, labeled as MKS192UCT (see Table 2), is performed on a static grid that has been derefined at the poles. More specifically, the AMR blocks touching both the polar axis and the inner radial boundary, which contain the smallest cells and therefore determine the global timestep, were forced to the coarsest refinement level. In addition to the settings described at the beginning of this section, the simulation domain spans r = [1.18 M, 10 000/3 M] and the extent in θ and ϕ subtends the whole 4π steradians solid angle. The resolution at the base level is N_{r} × N_{θ} × N_{ϕ} = 128 × 48 × 48, and the simulation contains 3 AMRlevels, giving an effective resolution of 512 × 192 × 192 cells. As an additional measure to prevent the timestep penalization from the poles, this time the MKS ϑ_{0}parameter is set to 0.25.
Figure 13 shows a snapshot of simulation MKS192UCT at t = 10 000 M. The AMR blocks, of 32 × 8 × 8 cells each, are drawn in white lines. The left panel of that figure shows the logarithm of density, while the right panel shows the MRIquality factor. It can be observed that the simulation grid is able to maintain in the disk without requiring a similarly high resolution at the polar regions, which are occupied by a smooth lowdensity outflow and where MRIdriven angular momentum transport does not play an important role. It is necessary to keep in mind that places where within the disk correspond to regions where the magneticfield changes sign (see Eq. (42)).
Fig. 13. Logarithmic restframe density (left panel) and MRI quality factor (right panel) for simulation MKS192UCT, showing the AMR blocks of N_{r} × N_{θ} × N_{ϕ} = 32 × 8 × 8 cells each. The grid is able to concentrate resolution in the disk in order to achieve in the torus, as recommended in Sano et al. (2004), while saving computational costs by derefining the polar regions. 
As can be seen in Fig. 14, the massaccretion rate Ṁ and magnetic flux through the horizon Φ_{BH} (Tchekhovskoy et al. 2011) are consistent with those from the highest resolution of Porth et al. (2017) (labeled here as MKS192fCT, see Table 2), which used the same initial condition and the numerical methods validated in that work including the fluxCT scheme. In both cases, the behavior consists of a transient growth (t < 5000 M) followed by a quasistationary state after the saturation of MRI (t > 5000 M).
Fig. 14. Massaccretion rate (top) and magnetic flux through the horizon (bottom) for the 3D simulations shown in this work in MKS and CKS coordinates using UCT2 (see Table 2), and for the highestresolution simulation in MKS coordinates shown in Porth et al. (2017), using fluxCT. 
As an additional comparison tool, we compute time and diskaveraged profiles similar to those shown in Beckwith et al. (2008), Shiokawa et al. (2012), White et al. (2016), Porth et al. (2017). For a quantity q(t, x^{i}), these averages are defined as
were θ_{min} = π/3, θ_{min} = 2π/3, t_{min} = 5000 M and t_{min} = 10 000 M.
As shown in Fig. 15, the averaged profiles of MKS192UCT show a reasonable quantitative agreement with those of MKS192fCT, despite the use of different methods and resolution. It should be noted, however, that the angular resolution in theta, which is essential to capture the MRI, is effectively the same for the two simulations performed in MKS coordinates, due to the use of three AMR levels in the former. As will be shown in Sect. 6.3, a better agreement is found with simulation CKS8UCT, which employs CKS coordinates, with differences arising from well understood reasons. The fact that such a better agreement can be found between simulations using different coordinate systems, and thus completely different spatial discretizations, shows that the differences between MKS192fCT and MKS192UCT are likely more a consequence of the method employed to evolve the magnetic field rather than of the slightly different radial resolution.
Fig. 15. Disk and timeaveraged profiles of rest mass density (a), ϕcomponent of the 4velocity (b), and inverse plasma β (c), in the interval 5000−10 000 M, for the 3D simulations of magnetized accretion referred to in this work (see Table 2). Shaded regions show the standard deviation from the average value for each simulation. (These are, however, hardly noticeable for u^{ϕ} due the small standard deviation and the excellent agreement between the three curves.) After taking into account the different numerical methods and coordinates employed, the three simulations show reasonably consistent profiles, with the deviations in density and plasma β at r_{KS} < 10 M for simulation CKS8UCT probably caused by a poor resolution of the MRI in that region (see also Fig. 16). 
6.3. 3D CKS simulations
Despite the strategies to avoid slow time steps from polar regions described in Sect. 6.2 and the selfconsistent treatment of poles by rearranging the grid topology described in Sect. 4.5, the poles of a coordinate system still represent an inhomogeneity in the numerical domain that could in principle introduce artifacts or a directional bias in the simulation. In addition, although logarithmic grids in polar coordinates ensure machineprecision conservation of angular momentum^{3} and are very efficient at resolving smallscale features when the interesting dynamics occur near the origin, they loose resolution far away from it.
Therefore, in order to study systems that are extended in space and for which an artificial directional bias could be misinterpreted as a physical property of the system, it might be useful to resort to coordinate systems that are more spatially homogeneous. An obvious choice are coordinate systems based on Cartesian coordinates, which are often used for GRMHD simulations in full general relativity (see e.g., Giacomazzo & Rezzolla 2007; Etienne et al. 2015). However, the necessity to properly resolve the black hole horizon and the MRI renders large scale (∼10^{3} r_{g}) simulations impossible, unless some form of mesh refinement is used.
In this section we present a simulation performed in CKS coordinates, labeled as CKS8UCT in Table 2. A combination of AMR and static mesh refinement is used to ensure sufficient resolution at the black hole horizon and the MRI turbulence while following the jet with dynamic mesh refinement.
The domain is structured as a set of nested boxes for which a different maximum refinement level is specified. The highest AMR level can be reached only by the innermost box, which contains the black hole horizon. In order to follow the jet dynamics, refinement is triggered by variations in the plasma magnetization σ = b^{2}/ρ and the density ρ, using the Löhner scheme. In order to limit the overhead by alternating refinement and derefinement of the same regions, regridding is performed only every 1000 iterations (since hierarchical timestepping has not yet been implemented in BHAC, these correspond to the timestep of the finest grid). The simulation employs 8 AMR levels, with a base resolution of N_{x} × N_{y} × N_{z} = 96 × 96 × 192, and extending over x, y ∈ [ − 500 M, 500 M], and z ∈ [ − 1000 M, 1000 M]. The magnetic field is evolved using the UCT2 algorithm. In order to prevent an unphysical outflow from the black hole interior, we apply cutoff values for the density and pressure in the region r < 0.5(r_{H−} + r_{H+}), where r_{H−} and r_{H+} are the location of the inner and outer event horizons. In particular, we set ρ_{cut} = 10^{−2} and p_{cut} = 10^{−4}.
Figure 16 displays different cuts of the simulation at t = 10 000 M. The two left panels are horizontal cuts at z = 100 M (panel a) and z = 0 (panel b), and the right panel is a vertical cut at y = 0 (panel c). In panel a it is possible to appreciate a crosssection of the jet, which is now completely free to wobble and to change shape independently of the coordinate surfaces. Panel c shows the automated mesh refinement following the evolution of the jet, indicated by a high plasma magnetization σ = b^{2}/ρ, as it propagates in the ±z direction. Also in this panel it can be noticed that the shape of the jet does not appear constrained by the coordinate surfaces. In a future work we will provide a more detailed comparison of the jet dynamics in simulations using Cartesian and spherical grids.
Fig. 16. Cuts of the simulation CKS8UCT (see Table 2) at t = 10 000 M, showing logarithmic plasma magnetization σ = b^{2}/ρ and MRI quality factor in the θdirection, . Panel a: cross section of the jet at z = 100 M. Panel b: MRIquality factors in the disk at z = 0. A magenta circumference marks the black hole outer horizon. Panel c: AMR blocks, of 16^{3} cells, following the propagation of the jet. 
At the same time as the jet is resolved with such detail, panel b shows that in most of the disk MRI is resolved with high quality factors, being the only exception the region for which r_{KS} < 10 M.
To quantitatively compare the results of this Cartesian simulation with those presented in the previous sections, also in this case we compute the massaccretion rate and magnetic flux through the horizon as well as time and diskaveraged profiles in the same intervals as mentioned in Sect. 6.2.
In Fig. 14, the massaccretion rate shows a remarkably consistent behavior for all of the three simulations, with variations being a consequence of the turbulent nature of the process. On the other hand, the absolute magnetic flux through the horizon, though overall consistent in magnitude, shows significantly less variations for the Cartesian case. This is accompanied by a smaller maximum in the initial transient growth. This could likely be attributed to a poorer resolution of the disk region which hinders magneticfield amplification due to MRI. In fact, although overall high quality factors are obtained within the disk, the decrease in due to the density increase close to the black hole is not accompanied by a corresponding increase in resolution as is the case for spherical polar coordinates (cf. panel b of Fig. 16 and right panel of Fig. 13). A more adequate resolution could be achieved by allowing higher AMR levels in this region.
The averaged profiles shown in Fig. 15 are as well in good agreement with the other 3D simulations, especially with MKS192UCT, which employs the same algorithm to evolve the magnetic field. The agreement is practically perfect for the ϕcomponent of the 4velocity in all three simulations and it remains mostly within one standard deviation between simulations CKS8UCT and MKS192UCT. The most important deviations in density and plasma β at r_{KS} < 10 M with respect to MKS192UCT could likely be attributed to the poorer resolution of MRI in that region, which hinders the amplification of the magnetic field and the angular momentum transport leading to an accumulation of mass in that region. In fact, panel c of Fig. 15 is consistent with a result from Sano et al. (2004), namely that a poor resolution leads to values of the magnetic pressure (and thus of β^{−1}) smaller than those expected at the saturation of MRI.
Finally, Table 3 lists several properties of runs CKS8UCT and MKS192UCT at the end of the simulation, related to how well they can resolve physics. These are: number of cells resolving the horizon, total cell population, timestep and average volume of occupied cells, that is, the average proper volume of cells which contain matter coming from the disk, which is identified by a passively advected tracer. Unsurprisingly, it can be seen that, although MKS simulations are able to resolve better the horizon and reach higher MRI quality factors in the disk region, overall the domain is more resolved in the CKS case. As a consequence, simulations in Cartesian coordinates could be useful to study the largescale effects on the jet produced by finer features arising from Kelvin Helmholtz or kink instabilities, as those visible in Fig. 16.
Comparison between Cartesian and spherical runs for the same accretion problem.
Furthermore the small sizes of cells close to the polar axis or the outer event horizon produces a penalization in timestep for MKS simulations which is absent for the Cartesian case, which is thus able to advance more physical time per iteration. However, the price to pay for a more resolved domain in the latter is a much larger cell population which currently is updated simultaneously, significantly increasing the computational cost of these simulations. In the future, this limitation will be overcome to some extent by adopting a hierarchical time step as is done in several AMRcodes (see e.g., Cunningham et al. 2009; Liska et al. 2018).
In summary, due to the advantages mentioned above, Cartesian adaptive meshes appear as a very interesting resource to study largescale jet propagation in simulations which self consistently include the jet engine, as well as other systems for which no symmetry is to be assumed a priori, as tidal disruption events or precessing disks. Currently, the Cartesian simulation described above is being used to model the appearance of the jetlaunching region in M 87 (Davelaar et al. 2019).
7. Conclusion
We described in detail new additions to the GRMHD code BHAC, namely three CT algorithms based on staggered meshes, as well as AMR strategies compatible with them. The variants of CT implemented in BHAC are the arithmetic average of Balsara & Spicer (1999, BS) and the two upwind schemes by Londrillo & Del Zanna (2004) and Del Zanna et al. (2007) (UCT), which allow the use of highorder reconstructions.
In order for the divergence of the magnetic field to be zero to machine precision across coarsefine boundaries and during the creation and destruction of blocks caused by AMR, special divergencepreserving restriction and prolongation operators are required. We derived and employ a class of such operators that generalises those obtained for Cartesian geometry by Tóth & Roe (2002). In addition, we presented technical information on the data structures used, the ghostcell exchange which was redesigned for staggered variables, and the treatment of the poles in spherical and cylindrical coordinates, as well as on divergencepreserving boundary conditions.
We validated these new additions by showing the results of tests commonly used in the community, specifically, the loop advection test by Gardiner & Stone (2005) and the cylindrical explosion by Komissarov (1999) in special relativity, and the magnetised stationary torus by Komissarov (2006) in general relativity. We observed that, in agreement with analogous tests present in the literature, UCT methods generate less spurious oscillations than the fluxCT and BS methods, both still widely used in GRMHD and MHD codes, and that the algorithm converges at the expected rate.
By performing 2 and 3dimensional simulations of magnetized accretion onto a black hole, we could notice that the BS method can give different results in physically relevant problems as compared to UCT methods. Due to its upwind properties and to the possibility of employing highorder reconstructions in UCT methods, we decided to use the latter schemes for future simulations. In addition, we showed that the code’s AMR capabilities can be exploited in black hole accretion simulations in order to eliminate the penalization in timestep caused by the small width of the cells around the pole in spherical coordinates, while maintaining an optimal resolution at the turbulent equatorial regions.
As an example of a problem inaccessible, or at least extremely expensive for nonAMR codes, which can be performed using the new capabilities added to the BHAC infrastructure, we simulated black hole accretion in Cartesian KerrSchild coordinates. The use of a Cartesian mesh could permit the study of jet dynamics including self consistently the black hole engine, while avoiding any possible directional bias introduced in the mesh by the presence of a polar axis. Furthermore, AMR can be used to accurately simulate magnetohydrodynamical instabilities between the disk wind and the jet occurring during its propagation, making similar setups extremely useful for selfconsistently modeling sources such as M 87 (Davelaar et al. 2019) and Cen A (Fromm et al., in prep.), as will be shown in future work. As mentioned before, Cartesian coordinates and AMR could be very useful also in accretion scenarios without symmetries, such as tidal disruption events or precessing disks.
Recent developments in the code that will be presented as well in forthcoming publications include an accurate modeling of electron thermodynamics (Mizuno, in prep.) and the addition of a module for resistive GRMHD (Ripperda et al., in prep.) which also employs the staggered grid infrastructure developed in the present work. A set of comprehensive tests showing consistency between results obtained with several stateoftheart GRMHD codes including BHAC for the same accretion problem will be presented in Porth et al. (2019).
All of these tools, together with the current capabilities of BHAC, are meant to contribute to a detailed modeling of strongfield field phenomena in astrophysics, which is becoming increasingly relevant for international efforts such as the EHT and BlackHoleCam (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f).
This quantity is often also indicated as W (Antón et al. 2006; Rezzolla & Zanotti 2013).
By adding to one cell the numerical flux that was subtracted from its neighbors (see Sect. 3.1), finitevolume methods achieve machineprecision conservation of the conserved variables (see e.g., Leveque 2002; Rezzolla & Zanotti 2013). When solving the GRMHD equations in polar coordinates, these include the covariant momentum in the ϕdirection.
Acknowledgments
We would like to thank Alejandro Cruz Osorio, Lukas Weih, David Kling, Jonas Köhler, and Mariafelicia de Laurentis for useful discussions. This research is supported by the ERC synergy grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant No. 610058), the “NewCompStar” COST Action MP1304, the LOEWEProgram in HIC for FAIR, and the European Union’s Horizon 2020 Research and Innovation Programme (Grant 671698, call FETHPC12014, project ExaHyPE). During the completion of this work, HO was supported in part by a CONACYTDAAD scholarship. The simulations were performed in part on the SuperMUC cluster at the LRZ in Garching, the LOEWE cluster at the CSC in Frankfurt, the Iboga cluster at the ITP Frankfurt, the HazelHen cluster at the HLRS in Stuttgart, as well as on the Dutch National Supercomputing cluster Cartesius, funded by the NWO computing grant 164. We acknowledge technical support from Thomas Coelho.
References
 Abdujabbarov, A. A., Rezzolla, L., & Ahmedov, B. J. 2015, MNRAS, 454, 2423 [NASA ADS] [CrossRef] [Google Scholar]
 Acker, F., de Borges, R. B. R., & Costa, B. 2016, J. Comput. Phys., 313, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Akiyama, K., Lu, R.S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Alcubierre, M. 2008, Introduction to 3+1 Numerical Relativity (Oxford, UK: Oxford University Press) [Google Scholar]
 Aloy, M. A., & Rezzolla, L. 2006, ApJ, 640, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Anile, A. M. 1990, Relativistic Fluids and Magnetofluids (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
 Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723 [NASA ADS] [CrossRef] [Google Scholar]
 Anninos, P., Bryant, C., Fragile, P., et al. 2017, ApJS, 231, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Antón, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Baiotti, L., Hawke, I., Montero, P. J., et al. 2005, Phys. Rev. D, 71, 024035 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S. 2001, J. Comput. Phys., 174, 614 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S. 2004, ApJS, 151, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., & Kim, J. 2004, ApJ, 602, 1079 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [Google Scholar]
 Beckwith, K., Hawley, J. F., & Krolik, J. H. 2008, ApJ, 678, 1180 [NASA ADS] [CrossRef] [Google Scholar]
 Britzen, S., Fendt, C., Witzel, G., et al. 2018, MNRAS, 478, 3199 [NASA ADS] [Google Scholar]
 Broderick, A. E., & Loeb, A. 2006, MNRAS, 367, 905 [NASA ADS] [CrossRef] [Google Scholar]
 Bugli, M., Guilet, J., Müller, E., et al. 2018, MNRAS, 475, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cowling, T. G. 1933, MNRAS, 94, 39 [NASA ADS] [Google Scholar]
 Cunningham, A. J., Frank, A., Varnière, P., Mitran, S., & Jones, T. W. 2009, ApJS, 182, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Davelaar, J., Mościbrodzka, M., Bronzwaer, T., & Falcke, H. 2018, A&A, 612, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davelaar, J., Olivares, H., Porth, O., et al. 2019, ArXiv eprints [arXiv:1906.10065] [Google Scholar]
 de Avellar, M. G. B., Porth, O., Younsi, Z., & Rezzolla, L. 2018, MNRAS, 474, 3967 [NASA ADS] [CrossRef] [Google Scholar]
 De Villiers, J.P., & Hawley, J. F. 2003, ApJ, 589, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Doeleman, S. S., Weintroub, J., Rogers, A. E. E., et al. 2008, Nature, 455, 78 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Duez, M. D., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2005, Phys. Rev. D, 72, 024028 [NASA ADS] [CrossRef] [Google Scholar]
 Etienne, Z. B., Liu, Y. T., & Shapiro, S. L. 2010, Phys. Rev. D, 82, 084031 [NASA ADS] [CrossRef] [Google Scholar]
 Etienne, Z. B., Paschalidis, V., Liu, Y. T., & Shapiro, S. L. 2012, Phys. Rev. D, 85, 024013 [NASA ADS] [CrossRef] [Google Scholar]
 Etienne, Z. B., Paschalidis, V., Haas, R., Mösta, P., & Shapiro, S. L. 2015, Class. Quantum Grav., 32, 175009 [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019c, ApJ, 875, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019d, ApJ, 875, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019e, ApJ, 875, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019f, ApJ, 875, L6 [Google Scholar]
 Fambri, F., Dumbser, M., Köppel, S., Rezzolla, L., & Zanotti, O. 2018, MNRAS, 477, 4543 [NASA ADS] [Google Scholar]
 Fish, V. L., Johnson, M. D., Doeleman, S. S., et al. 2016, ApJ, 820, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Galeazzi, F., Kastaun, W., Rezzolla, L., & Font, J. A. 2013, Phys. Rev. D, 88, 064009 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 458 [NASA ADS] [CrossRef] [Google Scholar]
 García, A., Galtsov, D., & Kechkin, O. 1995, Phys. Rev. Lett., 74, 1276 [NASA ADS] [CrossRef] [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [Google Scholar]
 Giacomazzo, B., & Rezzolla, L. 2007, Class. Quantum Grav., 24, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Goddi, C., Falcke, H., Kramer, M., et al. 2017, Int. J. Mod. Phys. D, 26, 1730001 [NASA ADS] [CrossRef] [Google Scholar]
 Gottlieb, S., & Shu, C. 1998, Math. Comp., 67, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [Google Scholar]
 Hartle, J. B., & Thorne, K. S. 1968, ApJ, 153, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984, ApJ, 277, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Issaoun, S., Johnson, M. D., Blackburn, L., et al. 2019, ApJ, 871, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, J. I. 1997, ApJ, 478, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Keppens, R., Meliani, Z., van Marle, A. J., et al. 2012, J. Comput. Phys., 231, 718 [Google Scholar]
 Kim, J., Marrone, D. P., Roy, A. L., et al. 2018, ApJ, 861, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S. 1999, MNRAS, 303, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S. 2006, MNRAS, 368, 993 [NASA ADS] [CrossRef] [Google Scholar]
 Kudoh, S. K. D. M. K. S. T. 2000, ApJ, 536, 668 [NASA ADS] [CrossRef] [Google Scholar]
 Leveque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems (New York: Cambridge University Press) [CrossRef] [Google Scholar]
 Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81 [NASA ADS] [CrossRef] [Google Scholar]
 Löhner, R. 1987, Comput. Meth. Appl. Mech. Eng., 61, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Londrillo, P., & Del Zanna, L. 2004, J. Comput. Phys., 195, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Martins, J. R. R. A., Sturdza, P., & Alonso, J. J. 2003, ACM Trans. Math. Softw., 29, 245 [CrossRef] [Google Scholar]
 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2013, Science, 339, 49 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., Sadowski, A., & Narayan, R. 2014, MNRAS, 441, 3177 [NASA ADS] [CrossRef] [Google Scholar]
 Meliani, Z., Mizuno, Y., Olivares, H., et al. 2017, A&A, 598, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mizuno, Y., Nishikawa, K. I., Koide, S., Hardee, P., & Fishman, G. J. 2006, ArXiv eprints [arXiv:astroph/0609004] [Google Scholar]
 Mizuno, Y., Lyubarsky, Y., Nishikawa, K.I., & Hardee, P. E. 2012, ApJ, 757, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, Y., Younsi, Z., Fromm, C. M., et al. 2018, Nat. Astron., 2, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Mocz, P., Pakmor, R., Springel, V., et al. 2016, MNRAS, 463, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Phys. Rev. Lett., 122, 061101 [NASA ADS] [CrossRef] [Google Scholar]
 Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959 [NASA ADS] [CrossRef] [Google Scholar]
 O’Connor, E., & Ott, C. D. 2010, Class. Quantum Grav., 27, 114103 [NASA ADS] [CrossRef] [Google Scholar]
 Olivares, H., Younsi, Z., Fromm, C. M., et al. 2018, ArXiv eprints [arXiv:1809.08682] [Google Scholar]
 Olivares Sánchez, H., Porth, O., & Mizuno, Y. 2018b, J. Phys. Conf. Ser., 1031, 012008 [CrossRef] [Google Scholar]
 Özel, F., Psaltis, D., & Narayan, R. 2000, ApJ, 541, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., & Komissarov, S. S. 2015, MNRAS, 452, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4 [Google Scholar]
 Porth, O., Olivares, H., Mizuno, Y., et al. 2017, Comput. Astrophys. Cosmol., 4, 1 [Google Scholar]
 Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Psaltis, D., Özel, F., Chan, C.K., & Marrone, D. P. 2015, ApJ, 814, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Psaltis, D., Wex, N., & Kramer, M. 2016, ApJ, 818, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., & Rezzolla, L. 2012, A&A, 547, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Radice, D., Rezzolla, L., & Galeazzi, F. 2014, MNRAS, 437, L46 [NASA ADS] [CrossRef] [Google Scholar]
 Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford, UK: Oxford University Press) [CrossRef] [Google Scholar]
 Rezzolla, L., & Zhidenko, A. 2014, Phys. Rev. D, 90, 084009 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ruffert, M., Janka, H.T., & Schaefer, G. 1996, A&A, 311, 532 [NASA ADS] [Google Scholar]
 Sano, T., Inutsuka, S.I., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shiokawa, H., Dolence, J. C., Gammie, C. F., & Noble, S. C. 2012, ApJ, 744, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Spiteri, R. J., & Ruuth, S. J. 2002, SIAM J. Numer. Anal., 40, 469 [CrossRef] [Google Scholar]
 Squire, W., & Trapp, G. 1998, SIAM Rev., 40, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Stergioulas, N., & Friedman, J. L. 1995, ApJ, 444, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Suresh, A., & Huynh, H. T. 1997, J. Comput. Phys., 136, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Metzger, B. D., Giannios, D., & Kelley, L. Z. 2014, MNRAS, 437, 2744 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd edn. (SpringerVerlag) [Google Scholar]
 Toth, G. 2000, J. Comput. Phys., 161, 605 [Google Scholar]
 Tóth, G., & Roe, P. 2002, J. Comput. Phys., 180, 736 [NASA ADS] [CrossRef] [Google Scholar]
 Tóth, G., van der Holst, B., Sokolov, I. V., et al. 2012, J. Comput. Phys., 231, 870 [Google Scholar]
 Čada, M., & Torrilhon, M. 2009, J. Comput. Phys., 228, 4118 [NASA ADS] [CrossRef] [Google Scholar]
 Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Wielgus, M., Fragile, P. C., Wang, Z., & Wilson, J. 2015, MNRAS, 447, 3593 [NASA ADS] [CrossRef] [Google Scholar]
 Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302 [Google Scholar]
 Younsi, Z., Zhidenko, A., Rezzolla, L., Konoplya, R., & Mizuno, Y. 2016, Phys. Rev. D, 94, 084025 [NASA ADS] [CrossRef] [Google Scholar]
 Zanotti, O., & Dumbser, M. 2015, Comput. Phys. Commun., 188, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Zanotti, O., Fambri, F., & Dumbser, M. 2015, MNRAS, 452, 3010 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
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 F^{i} (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). 

In the text 
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) are computed using the quantities in the boxes. 

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

In the text 
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 ydirection 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. 

In the text 
Fig. 5. Ghostcell exchange for blocks (a) at the same resolution level and for coarsefine 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. 

In the text 
Fig. 6. Ghostcell 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 ghostcell exchange at the z pole happens as for the north pole. 

In the text 
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 fluxCT. 

In the text 
Fig. 8. Cylindrical explosion. ycomponent of the magnetic field at t = 4.0 when using the BS algorithm (left) and UCT2 (right). 

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

In the text 
Fig. 10. Rest mass density, plasma beta and divergence of the magnetic field at t = 0 and t = 100 for the run with N_{r} = 400 of the Komisssarov torus test. 

In the text 
Fig. 11. L_{∞} and L_{1} norms of the error in density ρ for the stronglymagnetised Komissarov torus at t = 70 M, i.e., after one orbital period. The scheme shows second order convergence. 

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

In the text 
Fig. 13. Logarithmic restframe density (left panel) and MRI quality factor (right panel) for simulation MKS192UCT, showing the AMR blocks of N_{r} × N_{θ} × N_{ϕ} = 32 × 8 × 8 cells each. The grid is able to concentrate resolution in the disk in order to achieve in the torus, as recommended in Sano et al. (2004), while saving computational costs by derefining the polar regions. 

In the text 
Fig. 14. Massaccretion rate (top) and magnetic flux through the horizon (bottom) for the 3D simulations shown in this work in MKS and CKS coordinates using UCT2 (see Table 2), and for the highestresolution simulation in MKS coordinates shown in Porth et al. (2017), using fluxCT. 

In the text 
Fig. 15. Disk and timeaveraged profiles of rest mass density (a), ϕcomponent of the 4velocity (b), and inverse plasma β (c), in the interval 5000−10 000 M, for the 3D simulations of magnetized accretion referred to in this work (see Table 2). Shaded regions show the standard deviation from the average value for each simulation. (These are, however, hardly noticeable for u^{ϕ} due the small standard deviation and the excellent agreement between the three curves.) After taking into account the different numerical methods and coordinates employed, the three simulations show reasonably consistent profiles, with the deviations in density and plasma β at r_{KS} < 10 M for simulation CKS8UCT probably caused by a poor resolution of the MRI in that region (see also Fig. 16). 

In the text 
Fig. 16. Cuts of the simulation CKS8UCT (see Table 2) at t = 10 000 M, showing logarithmic plasma magnetization σ = b^{2}/ρ and MRI quality factor in the θdirection, . Panel a: cross section of the jet at z = 100 M. Panel b: MRIquality factors in the disk at z = 0. A magenta circumference marks the black hole outer horizon. Panel c: AMR blocks, of 16^{3} cells, following the propagation of the jet. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.