Free Access
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/0004-6361/201935559
Published online 05 September 2019

© ESO 2019

1. Introduction

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2. Mathematical formulation

2.1. Equations of GRMHD

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

(1)

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

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

(2)

and are related as where ϵμναβ are the components of the Levi-Civita tensor. Therefore, in a general frame moving with 4-velocity Uμ

(3)

and

(4)

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

(5)

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

(6)

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

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

(7)

where htot := h + b2/ρ is the total specific enthalpy and ptot = p + b2/2 is the total pressure, both in the fluid frame.

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

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

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

(8)

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 3-momentum Si := nμγνiTμν, the total energy U := nμnνTμν and the spatial stress tensor Wij := γμiγνjTμν.

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

(9)

in addition to the solenoidal constraint for the magnetic field . Here, γ is the determinant of the induced three-metric, and the vectors of conserved quantities U, the fluxes Fi, and the sources S are given by

(10)

where is the 3-velocity, Γ := −uμnμ is the Lorentz factor2, 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 Fi, knowledge of the primitive variables P = [ρνi,p,Bi] 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 horizon-scale VLBI images of Sgr A*, as well as the study of quasi-periodic oscillations in accretion disks around neutron stars (de Avellar et al. 2018).

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

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

Table 1.

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 Kerr-Schild coordinates (MKS). These coordinates, from now on labeled as xμ = (tMKS, s, ϑ, ϕMKS), were introduced by McKinney & Gammie (2004) in order to simulate large domains and, at the same time, concentrate the resolution near the black hole and on the equatorial plane. They are related to the standard Kerr-Schild coordinates xμ  =  (tKS, rKS, θKS, ϕKS) by the transformation

(11)

where R0 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 Kerr-Schild coordinates.

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

3. Numerical methods

3.1. Finite-volume scheme

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

To obtain the discretized equations, the simulation domain is divided into control volumes and the system is integrated over each of them. This leads to evolution equations for the average of the conserved quantities inside each cell

(12)

Quantities indicated as, say, F1ΔS1|i + 1/2, j, k represent integrals of the fluxes over the surfaces ΔS1|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 Fn (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 Δxi is the grid spacing in each direction and is Kronecker’s delta. Fn is obtained through the approximate solution of a Riemann problem at the interface, and static integrals such as cell volumes, interface areas and barycenter positions are calculated at initialization using fourth-order Simpson’s rule and stored in memory.

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

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

3.2. Constrained transport

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

(13)

is equivalent to ∇ ⋅ B = 0 in the continuous limit. When the limit is not taken, that is, for the finite-volume case, it follows from the divergence theorem that the average value of the divergence 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

(14)

thumbnail Fig. 1.

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

Open with DEXTER

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

(15)

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

(16)

Since each of the line integrals of the electric field is shared by two faces, but appears with opposite sign in the time update formula for each of them, the rate of change of , 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 divergence-free 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 divergence-free evolution to machine precision was first proposed by Yee (1966). However, although his method was widely known in engineering, a staggered grid was not applied in GRMHD until the work of Evans & Hawley (1988).

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

3.3. Arithmetic averaging (BS)

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

(17)

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

(18)

where is the mean value of the integral from Eq. (15) over the face at j + 1/2. To second-order accuracy, this integral takes the value 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

(19)

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

3.4. Upwind constrained transport (UCT)

Upwind Constrained Transport (UCT) is a method to evolve the magnetic field, first proposed in Londrillo & Del Zanna (2004). As a CT method, like BS, it maintains the divergence of the magnetic field to machine precision; however, it furthermore aims to accurately reproduce the magnetic-field continuity and transport properties by using limited reconstructions. Also, in contrast to arithmetic averaging, UCT reduces to the correct 1-dimensional limit when symmetry in the other two directions is assumed. Two variants of UCT are implemented in BHAC. These are presented in Londrillo & Del Zanna (2004) and in Del Zanna et al. (2007), and will be referred, respectively, as UCT1 and UCT2. Now we describe the procedure to obtain ℰ in each of these algorithms. To simplify the notation, we rename (x1, x2, x3) 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:

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

Open with DEXTER

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

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

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

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

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

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

    (20)

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

The implementation of UCT2 is as follows:

  1. At the time of calculating the numerical fluxes at each interface, store the characteristic speeds 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 x-interface.

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

    (21)

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

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

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

    (22)

    Again, we approximate ℰz = EzΔ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 Bi are reconstructed at the edges. We instead, factor the square root of the metric determinant and multiply the whole expression for Ez by it at the end. This makes the scheme more consistent when evolving a radial magnetic field in the polar regions, and simplifies the AMR operations at the poles, as will be seen in Sect. 4.5.

4. CT-adapted AMR

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

4.1. Notation

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

thumbnail Fig. 3.

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

Open with DEXTER

The coarse coordinate increment in the direction i is denoted ΔXi and the correspondent fine increment Δxi = ΔXi/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 x1 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 ΔS2, 0, 0 is subdivided are labeled Δs2, ±1, ±1. The same applies to ΔS−2, 0, 0 and in the other directions. The twelve fine faces Δs0, ±1, ±1, Δs±1, 0, ±1 and Δs±1, ±1, 0 are not part of any coarse face. The magnetic flux across Δs2, 1, 1 is

(23)

and that across ΔS2, 0, 0 is

(24)

The same applies to quantities defined at cell edges: ℰ2, 2, 0 denotes line integral of the electric field along the coarse edge centered at and ℰ2, 2, 1 that along the fine edge centered at .

In a 2D problem, where symmetry in the x3 direction is assumed, the integrals with respect to x3 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

(25)

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

4.2. Prolongation and restriction operators for face-allocated variables

4.2.1. Restriction formulas

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

(26)

The divergence in the newly created coarse cell will be

(27)

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

4.2.2. Prolongation formulas

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

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

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

(28)

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

(29)

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

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

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

(30)

The first requisite can be expressed as

(31)

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

(32)

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 nlmn, the sign of nlmn must reverse. This results in another set of relations for the coefficients,

(33)

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

(34)

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

(35)

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

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

(36)

where

(37)

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

(38)

where we rename Δx1 → x, Δx2 → y and Δx3 → 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

(39)

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

(40)

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

thumbnail Fig. 4.

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

Open with DEXTER

4.3. Restriction of edge-allocated variables

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

4.4. Parallelisation and ghost-cell exchange

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

thumbnail Fig. 5.

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

Open with DEXTER

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

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

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

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

4.5. Poles in spherical and cylindrical coordinates

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

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

thumbnail Fig. 6.

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

Open with DEXTER

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

4.6. Boundary conditions

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

5. Numerical tests

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

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

5.1. Loop advection

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

On a uniform background with ρ = 1, p = 1, νx = 0.2, and νy = 0.1, the problem consists on advecting a magnetic loop with radius R = 0.3, described by the potential Az = A0(R − r) for r ≤ R and Az = 0 for r >  R, where r2 = x2 + y2. A0 is chosen 10−3, so that β = pfluid/pmag = 106 and the loop is advected passively. The domain is x ∈ [ − 1, 1] and y ∈ [ − 0.5, 0.5] with periodic boundary conditions. The methods employed are the MP5 limiter, HLLE fluxes and the RK3 time-integrator.

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

The magnetic pressure at the initial condition and some snapshots of the evolution are depicted in Fig. 7. The upper-left panel of the same figure shows the divergence of the magnetic field for the UCT2 run at the final time t = 10.0, when the magnetic loop, and thus the refined region, has travelled a complete cycle across the domain. Divergence remains always zero at machine precision (∼10−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 flux-CT.

thumbnail Fig. 7.

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

Open with DEXTER

5.2. Cylindrical explosion

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

The domain extends over [ − 6, 6]×[−6, 6]. Plasma is initially at rest, with a constant magnetic field Bx = 0.1, By = Bz = 0. Defining r2 = x2 + y2, the domain is divided in three regions: a high density and pressure (ρin = 10−2, pin  =  1) region for r <  rin  =  0.8, a low density and pressure (ρout  =  10−4, pout  =  5 × 10−4) region for r >  rout  =  1.0, and a transition region for rin ≤ r ≤ rout where ρ = ρin(r/rin)α1 and p = pin(r/rin)α2, and the exponents α1 and α2 are such that p(r) and ρ(r) are continuous at rin and rout. The boundary conditions are periodic, although this is not relevant in practice, since the region of interest cannot be affected by signals coming form the boundaries during the simulation time.

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

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

thumbnail Fig. 8.

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

Open with DEXTER

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

thumbnail Fig. 9.

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

Open with DEXTER

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ϕ/ut = l0, and adopting a specific relation of the fluid and magnetic pressure with respect to the fluid enthalpy, the configuration is described by the equation

(41)

where 𝒲 = ln|ut| 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/pmag)c = βc, respectively. The parameters of the torus chosen here as a test problem are κ = 4/3, η = 4/3, l0 = 2.8, rc = 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 Kerr-Schild coordinates (i.e., MKS coordinates with ϑ0 = 0 and R0 = 0). The simulation domain spans over θ ∈ [π/5, 4π/5], rKS ∈ [0.95 rh,  50 M] and ϕ ∈ [0, π], where rKS is the radial Kerr-Schild coordinate and rh 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 Nr × Nθ × Nϕ = 100 × 60 × 40. The rotation period of the center of the torus is 68 M, and the simulation is carried out until t = 100 M, thus spanning nearly one and a half orbital periods. Although this kind of solutions are known to be unstable to the MRI and the Papaloizou-Pringle instability in 3D (Wielgus et al. 2015; Bugli et al. 2018), the time elapsed by our simulations still corresponds to the initial phase of the linear growth, and no chaotic behavior can be observed. In fact, the torus here presented is identical to Case A simulated by Wielgus et al. (2015), who reported perturbations at t = 100 M still at the ∼10−4 level.

To emulate vacuum regions outside the torus, a tenuous atmosphere with density and pressure given by , with ρmin = 10−5 and pmin = 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 two-step time integration with CFL factor 0.5.

Figure 10 shows the density distribution, the plasma beta and the divergence of the magnetic field at the initial state and at t = 100 M for the run with resolution Nr = 400. As can be seen, the distribution of the displayed quantities is very well preserved by the scheme. The L1 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.

thumbnail Fig. 10.

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

Open with DEXTER

thumbnail Fig. 11.

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

Open with DEXTER

6. Magnetized black hole accretion

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

For all the simulations in this section, the initial configuration is a torus in equilibrium around a black hole (Fishbone & Moncrief 1976) with spin parameter a = 0.9375, and thus with the outer event horizon located at r = 1.348 M, where r is the radial Kerr-Schild coordinate. The inner radius of the torus is at rin = 6 M, and its density maximum at rmax = 12 M (orbital period of 247 M at the density maximum). The density in the torus is normalized so that it takes 1.0 as its a maximum value. A single-loop poloidal magnetic field is built from the vector potential Aϕ ∝ max(ρ/ρmax − 0.2, 0) and is normalized in such a way that the highest magnetic pressure and the highest thermal pressure are related by the ratio β = pfluid, max/pmag, 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−4r−3/2 and fluid pressure pfl = 1/3 × 10−6r−5/2. We reset the density or the pressure whenever they fall below these floor values. Simulations were evolved using a two-step time-integrator.

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

Table 2.

Summary of the simulations referred in Sect. 6.

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 Nr  ×  Nθ = 200  ×  100 and an effective resolution of Nr  ×  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 high-resolution VLBI images of the SMBH candidates in M 87 and Sgr A*, since radiative models of these sources which are able to reproduce the properties of their radio spectrum consider that most of the radio emission originates at the funnel wall (Mościbrodzka et al. 2016).

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

Open with DEXTER

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

6.2. 3D MKS simulations

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

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

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

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

(42)

where b(i) and Δx(i) are, respectively, the magnetic field and the displacement in the ith direction in the orthonormal frame co-moving 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 MKS192-UCT (see Table 2), is performed on a static grid that has been de-refined 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 time-step, 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 Nr × Nθ × Nϕ = 128 × 48 × 48, and the simulation contains 3 AMR-levels, giving an effective resolution of 512 × 192 × 192 cells. As an additional measure to prevent the time-step penalization from the poles, this time the MKS ϑ0-parameter is set to 0.25.

Figure 13 shows a snapshot of simulation MKS192-UCT 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 MRI-quality 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 low-density outflow and where MRI-driven 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 magnetic-field changes sign (see Eq. (42)).

thumbnail Fig. 13.

Logarithmic rest-frame density (left panel) and MRI quality factor (right panel) for simulation MKS192-UCT, showing the AMR blocks of Nr × 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 de-refining the polar regions.

Open with DEXTER

As can be seen in Fig. 14, the mass-accretion 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 MKS192-fCT, see Table 2), which used the same initial condition and the numerical methods validated in that work including the flux-CT scheme. In both cases, the behavior consists of a transient growth (t <  5000 M) followed by a quasi-stationary state after the saturation of MRI (t >  5000 M).

thumbnail Fig. 14.

Mass-accretion 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 highest-resolution simulation in MKS coordinates shown in Porth et al. (2017), using flux-CT.

Open with DEXTER

As an additional comparison tool, we compute time- and disk-averaged 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, xi), these averages are defined as

(43)

were θmin = π/3, θmin = 2π/3, tmin = 5000 M and tmin = 10 000 M.

As shown in Fig. 15, the averaged profiles of MKS192-UCT show a reasonable quantitative agreement with those of MKS192-fCT, 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 CKS8-UCT, 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 MKS192-fCT and MKS192-UCT are likely more a consequence of the method employed to evolve the magnetic field rather than of the slightly different radial resolution.

thumbnail Fig. 15.

Disk- and time-averaged profiles of rest mass density (a), ϕ-component of the 4-velocity (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 rKS <  10 M for simulation CKS8-UCT probably caused by a poor resolution of the MRI in that region (see also Fig. 16).

Open with DEXTER

6.3. 3D CKS simulations

Despite the strategies to avoid slow time steps from polar regions described in Sect. 6.2 and the self-consistent 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 machine-precision conservation of angular momentum3 and are very efficient at resolving small-scale 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 (∼103rg) simulations impossible, unless some form of mesh refinement is used.

In this section we present a simulation performed in CKS coordinates, labeled as CKS8-UCT 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 σ = b2/ρ and the density ρ, using the Löhner scheme. In order to limit the overhead by alternating refinement and de-refinement of the same regions, re-gridding is performed only every 1000 iterations (since hierarchical time-stepping has not yet been implemented in BHAC, these correspond to the time-step of the finest grid). The simulation employs 8 AMR levels, with a base resolution of Nx × Ny × Nz = 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 cut-off values for the density and pressure in the region r <  0.5(rH− + rH+), where rH− and rH+ are the location of the inner and outer event horizons. In particular, we set ρcut = 10−2 and pcut = 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 cross-section 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 σ = b2/ρ, 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.

thumbnail Fig. 16.

Cuts of the simulation CKS8-UCT (see Table 2) at t = 10 000 M, showing logarithmic plasma magnetization σ = b2/ρ and MRI quality factor in the θ-direction, . Panel a: cross section of the jet at z = 100 M. Panel b: MRI-quality factors in the disk at z = 0. A magenta circumference marks the black hole outer horizon. Panel c: AMR blocks, of 163 cells, following the propagation of the jet.

Open with DEXTER

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 rKS <  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 mass-accretion rate and magnetic flux through the horizon as well as time- and disk-averaged profiles in the same intervals as mentioned in Sect. 6.2.

In Fig. 14, the mass-accretion 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 magnetic-field 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 MKS192-UCT, which employs the same algorithm to evolve the magnetic field. The agreement is practically perfect for the ϕ-component of the 4-velocity in all three simulations and it remains mostly within one standard deviation between simulations CKS8-UCT and MKS192-UCT. The most important deviations in density and plasma β at rKS <  10 M with respect to MKS192-UCT 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 CKS8-UCT and MKS192-UCT 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, time-step 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 large-scale effects on the jet produced by finer features arising from Kelvin Helmholtz or kink instabilities, as those visible in Fig. 16.

Table 3.

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 time-step 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 AMR-codes (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 large-scale 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 jet-launching 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 high-order reconstructions.

In order for the divergence of the magnetic field to be zero to machine precision across coarse-fine boundaries and during the creation and destruction of blocks caused by AMR, special divergence-preserving 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 ghost-cell exchange which was re-designed for staggered variables, and the treatment of the poles in spherical and cylindrical coordinates, as well as on divergence-preserving 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 flux-CT 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 3-dimensional 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 high-order 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 time-step 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 non-AMR codes, which can be performed using the new capabilities added to the BHAC infrastructure, we simulated black hole accretion in Cartesian Kerr-Schild 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 set-ups extremely useful for self-consistently 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 state-of-the-art 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 strong-field 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).


2

This quantity is often also indicated as W (Antón et al. 2006; Rezzolla & Zanotti 2013).

3

By adding to one cell the numerical flux that was subtracted from its neighbors (see Sect. 3.1), finite-volume methods achieve machine-precision 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 LOEWE-Program in HIC for FAIR, and the European Union’s Horizon 2020 Research and Innovation Programme (Grant 671698, call FETHPC-1-2014, project ExaHyPE). During the completion of this work, HO was supported in part by a CONACYT-DAAD 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

All Tables

Table 1.

Coordinates available in BHAC.

Table 2.

Summary of the simulations referred in Sect. 6.

Table 3.

Comparison between Cartesian and spherical runs for the same accretion problem.

All Figures

thumbnail Fig. 1.

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

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 3.

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

Open with DEXTER
In the text
thumbnail Fig. 4.

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

Open with DEXTER
In the text
thumbnail Fig. 5.

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

Open with DEXTER
In the text
thumbnail Fig. 6.

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

Open with DEXTER
In the text
thumbnail Fig. 7.

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

Open with DEXTER
In the text
thumbnail Fig. 8.

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

Open with DEXTER
In the text
thumbnail Fig. 9.

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

Open with DEXTER
In the text
thumbnail Fig. 10.

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

Open with DEXTER
In the text
thumbnail Fig. 11.

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

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 13.

Logarithmic rest-frame density (left panel) and MRI quality factor (right panel) for simulation MKS192-UCT, showing the AMR blocks of Nr × 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 de-refining the polar regions.

Open with DEXTER
In the text
thumbnail Fig. 14.

Mass-accretion 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 highest-resolution simulation in MKS coordinates shown in Porth et al. (2017), using flux-CT.

Open with DEXTER
In the text
thumbnail Fig. 15.

Disk- and time-averaged profiles of rest mass density (a), ϕ-component of the 4-velocity (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 rKS <  10 M for simulation CKS8-UCT probably caused by a poor resolution of the MRI in that region (see also Fig. 16).

Open with DEXTER
In the text
thumbnail Fig. 16.

Cuts of the simulation CKS8-UCT (see Table 2) at t = 10 000 M, showing logarithmic plasma magnetization σ = b2/ρ and MRI quality factor in the θ-direction, . Panel a: cross section of the jet at z = 100 M. Panel b: MRI-quality factors in the disk at z = 0. A magenta circumference marks the black hole outer horizon. Panel c: AMR blocks, of 163 cells, following the propagation of the jet.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.