Recovery schemes for
primitive variables in
General-relativistic magnetohydrodynamic (GRMHD) simulations are an important tool to study a variety of astrophysical systems such as neutron star mergers, core-collapse supernovae, and accretion onto compact objects. A conservative GRMHD scheme numerically evolves a set of conservation equations for ’conserved’ quantities. In order to advance the conservatives from one time step to the next, a set of primitive variables needs to be computed at every time step, which is typically achieved by solving a set of non-linear equations. This recovery procedure constitutes a core part of any conservative GRMHD scheme and it is closely tied to the equation of state (EOS) of the fluid. In the quest to include nuclear physics, weak interactions, and neutrino physics, state-of-the art GRMHD simulations employ finite-temperature, composition-dependent EOS. Different schemes have individually been proposed for such EOS; however, the recovery problem still remains a major source of error, failure and inefficiency in GRMHD simulations with advanced microphysics. The strengths and weaknesses of the different schemes when compared to each other remain unclear. Here we present the first systematic comparison of various recovery schemes used in different dynamical spacetime GRMHD codes for both analytic and tabulated microphysical EOS. We assess the schemes in terms of speed, accuracy, and robustness. We find large variations among the different schemes and that there is not a single ideal scheme. While the computationally most efficient schemes are less robust, the most robust schemes are computationally less efficient. More robust schemes may require an order of magnitude more calls to the EOS, which are computationally expensive. We propose an optimal strategy of an efficient three-dimensional Newton-Raphson scheme and a slower but more robust one-dimensional scheme as a fall-back.
NASA Einstein Fellow \altaffiliationNASA Einstein Fellow
Magnetized plasmas in gravitating systems where general relativity is important play a pivotal role in understanding a variety of phenomena across astrophysics, including accreting black holes in (X-ray) binaries and active galactic nuclei, short and long gamma-ray bursts, core-collapse supernovae, and compact-object mergers involving neutron stars. Various general-relativistic magnetohydrodynamic (GRMHD) codes have been developed to study these systems from first principles, both in stationary background spacetime (e.g., HARM (Gammie et al., 2003), Komissarov 2005, Antón et al. 2006, ECHO (Del Zanna et al., 2007), Athena++ (White et al., 2016)) and in fully dynamical spacetime (e.g., Duez et al. 2005, WhiskyMHD (Giacomazzo & Rezzolla, 2007), Cerdá-Durán et al. 2008, Anderson et al. 2008, Kiuchi et al. 2012, GRHydro (Mösta et al., 2014a), SpEC (Muhlberger et al., 2014), IllinoisGRMHD (Etienne et al., 2015)).
Conservative GRMHD schemes, such as the ones employed by the aforementioned codes, solve the equations of GRMHD as a set of conservation equations of the form (see Sec. 2)
where and are (known) quantities derived from the spacetime metric. Here, represents a vector of so-called conserved quantities, which are analytic functions of the so-called primitive variables , . Since the fluxes and source terms are functions of the primitives, at every time step of the conservative scheme, the given vector of conserved variables needs to be inverted to obtain the primitive variables in order to evolve to the next time step. This recovery procedure therefore represents a core part of every conservative GRMHD scheme.
Whereas the inversion from conservative to primitive quantities is known in terms of analytical relations in Newtonian MHD (and thus trivial), there is no known such inversion in closed form in GRMHD, and usually a set of non-linear equations needs to be solved to obtain the primitives . Complication arises from the equation of state (EOS), which captures the thermodynamic properties of the fluid and which is closely tied to the recovery problem. Most GRMHD codes have used simple analytic EOS so far, such as (piecewise) polytropic EOS or ideal gas EOS, in which case specific recovery schemes exist (Noble et al., 2006). Even in the context of simple analytic and well-behaved EOS, this recovery problem is one of the most error-prone key parts of any GRMHD evolution code.
Advancing realism in numerical simulations, particularly in fully dynamical spacetime applications like core-collapse supernovae and compact binary mergers, requires GRMHD codes to support composition-dependent, finite-temperature EOS, typically formulated in terms of rest-mass density, temperature, and electron fraction (3-parameter EOS). For example, the inclusion of weak interactions and (approximate) neutrino transport into GRMHD require such EOS. Furthermore, capabilities for composition-dependent, finite-temperature EOS are required for exploring the influence of different EOS on neutron star mergers; this may provide an important tool for inferring the unknown EOS of nuclear matter at high densities from observations of NS mergers by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and electromagnetic facilities. To date, only very few GRMHD codes and studies using composition-dependent EOS exist (e.g., Neilsen et al. 2014; Palenzuela et al. 2015; Mösta et al. 2014b, 2015; Siegel & Metzger 2017a, b).
Composition-dependent EOS are usually supplied to the code partially or entirely in form of tables. EOS calls by the code therefore involve table lookups with interpolation operations, which can be computationally expensive. Widely used recovery schemes that, in principle, support 3-parameter EOS, such as the 2D scheme in Noble et al. (2006) or the recovery scheme in Antón et al. (2006), require additional inversions using the EOS at every iteration step of the scheme from, e.g., specific enthalpy or pressure to temperature, and thus introduce many additional EOS table lookups; this potentially makes the recovery process computationally expensive and introduces many more operations and additional sources of error and/or failure that make the GRMHD scheme less robust.
In this paper, we discuss and compare several recovery schemes that are suitable for composition-dependent, finite-temperature EOS. Most of the schemes presented here have already been used in GRMHD simulations, in different codes and by different groups. Here, we provide the first comparison of these schemes in a well-defined test setting, i.e., independent of a specific application, and assess these schemes in terms of three cardinal criteria: speed (in terms of the number of iterations and the number of EOS calls needed to converge), accuracy, and robustness. Section 2 defines the recovery problem and Sec. 3 introduces the recovery schemes considered in this analysis. Results of our comparison are presented in Sec. 4 and Sec. 5 summarizes our conclusions.
2 The recovery problem
The equations of ideal GRMHD include energy and momentum conservation, baryon number conservation, lepton number conservation, and Maxwell’s equations (e.g., Antón et al. 2006),
where is the energy-momentum tensor, the 4-velocity, the baryon number density, the electron number density, and the dual of the Faraday electromagnetic tensor. The energy-momentum tensor in ideal GRMHD can be written as
where denotes pressure, specific enthalpy, with being the specific internal energy, and the magnetic field vector in the frame comoving with the fluid, , and the space-time metric.111In this paper, Greek indices take space-time values 0,1,2,3, whereas Roman indices represent the spatial components 1,2,3 only. Repeated indices are summed over. We assume that the thermodynamic properties of matter can be described by a composition-dependent, finite-temperature EOS, formulated as a function of density , where denotes the baryon mass, temperature , and electron fraction . Accordingly, we have added an evolution equation for the electron fraction (Eq. (4)) to the original set of ideal GRMHD evolution equations (2), (3), and (5). The terms on the right-hand side of Eqs. (2) and (4) represent source terms that reflect the evolution of the electron fraction due to weak interactions, which lead to the emission of neutrinos and anti-neutrinos that carry away energy and momentum from the system.
For numerical evolution, we adopt a 3+1 split of spacetime into non-intersecting spacelike hypersurfaces of constant coordinate time (Lichnerowicz, 1944; Arnowitt et al., 2008), writing the line element as
where denotes the lapse function, the shift vector, and the metric induced on every spatial hypersurface. The hypersurfaces are characterized by the timelike unit normal , where and . We adopt an Eulerian formulation of the equations of GRMHD in terms of the Eulerian observer, defined as the observer in the 3+1 foliation of spacetime moving with 4-velocity perpendicular to the hypersurfaces of constant coordinate time .
suitable for numerical evolution. Here, is the determinant of the spatial metric and
denotes the vector of conserved variables, where
and are the 3-vector components of the magnetic field as measured by the Eulerian observer; is called the conserved density, represent the conserved momenta, and is the conserved energy. For further reference, we note that
The Eulerian observer measures the fluid 4-velocity with corresponding 3-velocity
is the relative Lorentz factor between and , with . We note that the comoving and Eulerian magnetic field components are related by
where . Furthermore, we shall refer to
as the vector of primitive variables.
Conservative numerical GRMHD schemes evolve the state vector from one time step to the next using Eq. (8). This involves computing the flux terms and source terms for given , for which one needs to obtain the primitive variables from the conserved ones. While the conservative variables as a function of primitive variables, are given in analytic form by Eqs. (10)–(20), the inverse relation , i.e., the recovery of primitive variables from conservative ones, is not known in closed form; this rather requires numerical inversion of the aforementioned set of non-linear equations.
3 Recovery schemes
The recovery problem is, in principle, five-dimensional (5D), i.e., obtaining the primitives , , in Eq. (21) from Eqs. (10)–(20) involves numerical root-finding in a 5D space, given that the recovery of the magnetic field components and the electron fraction from (cf. Eq. (9)) is trivial. While 5D schemes as originally implemented in, e.g., the HARM code (Gammie et al., 2003), were later shown to be slow and inaccurate, reducing the dimensionality of the problem to 1D or 2D with the help of certain scalar quantities proved more efficient, with the 2D schemes being less pathological and therefore being recommended (Noble et al., 2006).
In this section, we first discuss the widely used 2D scheme by Noble et al. (2006) and extend it to 3-parameter EOS, which will serve as a standard reference (hereafter referred to as “2D NR Noble”). We then introduce a variant of the 2D scheme used in Antón et al. (2006), Aasi et al. (2014), Giacomazzo & Rezzolla (2007),Cerdá-Durán et al. (2008), suitable for the use of 3-parameter EOS (henceforth “2D NR”). Furthermore, we discuss several recovery schemes (3D and effective-1D) that have been recently proposed and that are used in current GRMHD codes with support for composition-dependent, finite-temperature EOS.
3.1 2D NR Noble et al. scheme
This scheme reduces the dimensionality of the recovery problem by making use of certain scalar quantities that can be computed from the conservatives; it is based on Eqs. (27) and (29) of Noble et al. (2006), which are two equations in the unknowns and . In our notation, they read
The pressure can be directly obtained in terms of the unknowns and only for simple analytic EOS such as the ideal gas EOS. In order to extend this scheme to 3-parameter EOS, at every iteration of the scheme we invert the specific enthalpy for given , , and with the help of the EOS to find the temperature and then compute using the EOS. This involves inversion with a 1D Newton-Raphson scheme, which falls back to bisection if the inversion fails.
Once converged, the final primitives , , and can be obtained from ,
and, using the EOS,
3.2 2D NR scheme
In the following, we present a 2D scheme starting from the same equations as the 2D schemes in Antón et al. (2006),Giacomazzo & Rezzolla (2007), Cerdá-Durán et al. (2008) and similar to Noble et al. (2006), but formulated in terms of unknowns that are more suitable for the use of modern microphysical (tabulated) EOS provided in terms of , , and . The advantage of our formulation is that additional inversion steps from specific enthalpy or specific internal energy to temperature are not necessary. Such additional inversions at every iteration of the main root-finding process are not required in the case of the ideal gas EOS, but introduce additional operations and sources of error or failure for 3-parameter EOS, and are computationally expensive.
Similar to the 2D Noble scheme (Sec. 3.1), this scheme reduces the dimensionality of the recovery problem to 2D by making use of the two scalar quantities (cf. Eq. (13)) and (cf. Eq. (14)). Setting as before and using the identity
Employing the EOS to compute the pressure from , we take Eqs. (27) and (28) as two equations in the unknowns and by setting , , with being calculated from Eq. (10), and solve using a 2D Newton-Raphson (NR) algorithm. Once converged, the final primitives are obtained from Eqs. (10), (24), and (25).
3.3 3D NR scheme
Here, is computed as
and in Eqs. (29) and (32) as well as in Eq. (31) are computed using the EOS with from Eq. (10). As this scheme also employs the temperature directly as an unknown, it does not require any inversions with the EOS. Once the system of Eqs. (29)–(31) has been solved with a 3D Newton-Raphson scheme, one recovers the final primitives from Eqs. (10), (24), and (25).
3.4 effective 1D schemes
3.4.1 Method of Palenzuela et al.
The scheme recently presented by Neilsen et al. (2014) and Palenzuela et al. (2015) reduces the dimensionality of the recovery problem to 1D, but requires an additional inversion step using the EOS at every iteration of the main scheme if the EOS is not provided in terms of . The scheme utilizes Brent’s method to solve
in the unknown
which is bounded by
Here, following Palenzuela et al. (2015) we have defined
Quantities with a hat in Eq. (33) are computed at every iteration step from and the conservatives by the following procedure:
Using Eq. (10), one obtains .
Using the EOS, one inverts to find the corresponding temperature and thus .
3.4.2 Method of Newman & Hamlin
The recovery scheme presented in Newman & Hamlin (2014) is an effective 1D method, as is the scheme of Palenzuela et al. (2015) and Neilsen et al. (2014) discussed above. It iterates over the fluid pressure to find a solution and as in the method of Palenzuela et al. (2015) requires and additional inversion step using the EOS for every main scheme iteration. It is independent of an initial guess as it starts with as a guess for the pressure. For tabulated EOS with a finite minimum pressure, we employ the pressure obtained from the minimum temperature of the EOS at given and , where is an initial guess for the present Lorentz factor.
Starting from the conserved quantities
(45) (46) (47) (48)
where , the scheme solves a cubic polynomial
where as before and
Evaluating the positive extremum of this equation, one can find a sufficient condition for the existence of a positive root to Eq. (49) as
one can show that for , with
the first root of
corresponds to the correct and physical solution of Eq. (49), from which one obtains .
One can now compute primitive quantities as
(57) (58) (59) (60)
Inverting with the help of the EOS for of Eq. (60) and given , one obtains the corresponding temperature and a new pressure guess .
One iterates over steps (i) and (ii) until the pressure has converged. When three successive pressure estimates are known, the Aitken acceleration scheme as described in Newman & Hamlin (2014) is used. If convergence is reached under the Aitken acceleration, we recalculate the variables as in step (ii). If not, the solver restarts from the beginning of step (i).
Once a converged value for has been acquired, including a consistent set of , , and from (iii), the remaining primitives can be obtained from Eq. (24).
4 Comparison of recovery schemes
4.1 Comparison setup and criteria
Comparing the different recovery schemes described in Sec. 3 proceeds by generating a parameter space of primitive variables , computing the corresponding conservative variables using Eqs. (10)–(13), and applying the recovery schemes to obtain sets of primitives that can be compared to the original ones. We perform this comparison of schemes for three different EOS, the analytic ideal gas EOS (where we arbitrarily set ), the partially analytic, partially tabulated Helmholtz EOS (Timmes & Arnett, 1999; Timmes & Swesty, 2000), and the fully tabulated Lattimer-Swesty EOS (Lattimer & Swesty, 1991) with incompressibility modulus (hereafter LS220), which serves as a widely used standard reference EOS in GRMHD simulations of neutron star mergers and core-collapse supernovae. We reduce the parameter space by varying the Lorentz factor instead of , and then choosing the orientation of the 3-velocity randomly at every point in the parameter space. The magnetic field vector is set by a specified absolute value and oriented relative to the 3-velocity vector by a specified angle. For the results presented below, we shall assume that and are aligned for simplicity; the results for the comparison of the schemes does not strongly depend on this relative orientation. For those recovery schemes that require initial guesses explicitly (such as 2D Noble, 2D NR, 3D NR) or implicitly (such as 1D Brent), we apply a random perturbation to the initial primitives and use the perturbed set as initial guess; this deviation typically represents a conservative bracket on the change from one time step to the next in an actual GRMHD evolution and should thus lead to a conservative upper limit on, e.g., the number of iterations or EOS calls required to converge to the solution. We have checked that the results presented here are essentially unchanged if an even higher perturbation of 10% is applied. We define convergence of a scheme by having surpassed a (maximum) relative error in the iteration variable(s). The resulting recovered primitives are compared to the original set and to the results from the other recovery schemes in terms of the following criteria:
Speed. The speed of a given recovery scheme is assessed by counting the number of iterations needed to converge to a solution. Additionally, we have implemented counters to evaluate the total number of EOS calls required until convergence. As table lookups can dominate the computational cost of the recovery scheme and with some schemes requiring additional EOS inversion steps at each iteration, monitoring the number of EOS calls represents a better proxy for the overall speed of a scheme than just counting the number of iterations.
Accuracy. We assess the accuracy of recovery by comparing the recovered primitives with the original ones and computing the mean relative deviation/error,
where only , , and are included in the summation, as and are recovered trivially. We define successful recovery of the primitives when all primitives have been recovered with a relative error smaller than .
Robustness. Robustness is assessed qualitatively in terms of the following criteria: independence of an initial guess, guaranteed convergence, and independence of thermodynamic derivatives.
Figures 1–3 show comparisons of the five different recovery schemes discussed in this paper. We present results for three different EOS: the ideal gas EOS, the LS220 tabulated EOS, and the partially-tabulated Helmholtz EOS. White spaces correspond to parameter sets for which the recovery process either failed or did not achieve the required accuracy (Sec. 4.1) within a certain number of iterations.
The top panels of Fig. 1 compare the number of iterations required for each scheme to converge to the correct primitive variables within the required accuracy (Sec. 4.1). The parameter range in density and temperature is chosen to cover the validity ranges of the chosen EOS and we assume representative values of and , where is the magnetic pressure.
For the ideal gas EOS, all schemes recover the entire parameter space in density and temperature with only a few iterations (). The 2D NR Noble scheme overall requires the smallest number of iterations.
For the LS220 EOS there is more variation among the different schemes. The 2D NR scheme recovers most of the parameter space with iterations but requires 15–30 iterations in the upper left corner as well as for densities above and temperatures below . The 3D NR scheme needs 15–25 iterations for regions at lower temperature, whereas the 1D Brent method performs well and recovers the entire parameters space with iterations. The 2D NR Noble and 1D Newman methods perform similarly with the exception that there are small regions where the schemes do not recover the primitives within the given accuracy at very high density and low temperature (white pixels in the bottom right corners). These points account for less than 1.5% of the total points. These failures arise because the additional inversion step of finding a temperature for the given relativistic enthalpy fails in the EOS table. We have investigated this in some detail and find that for all of these points the initial first step in the solver leads to relativistic enthalpy values outside of the EOS table. We have implemented a step size limitation for these cases but did not find a straightforward way to successfully apply it consistently to all points without additional side effects. While there are a number of ways in which a more sophisticated treatment of step sizes can be implemented, we have chosen to not finetune inversions for individual schemes too much with respect to others in order to ensure a fair comparison. The requirement of this additional inversion for temperature given the relativistic enthalpy is a feature/weakness of these schemes that limits their robustness.
For the Helmholtz EOS, the 2D and 3D NR solver recover most of the parameter space with less than 10 iterations. Both solvers however need iterations for very low temperatures. The 2D NR solver additionally needs iterations in the upper left corner. The extent of this region is similar to the results of the 2D NR solver for the LS220 EOS. The 1D Brent, 2D NR Noble and 1D Newman methods recover the correct solution with less than 10–12 iterations in the entire parameter space.
Based on the number of iterations to converge the 1D Brent method would be one obvious first choice as it recovers the entire parameter space for all EOS well with iterations.
In the bottom panels of Fig. 1, we show the number of EOS calls for each recovery scheme and the same three EOS as in the top panels of Fig. 1. While there is no drastic variation in the number of iterations to convergence for the different schemes, the number of EOS calls required by the different schemes varies dramatically. All schemes need a small number () of EOS calls to converge for the ideal gas EOS. For the LS220 EOS, the 2D and 3D NR schemes need less than EOS calls to converge in the entire parameter space. The 1D Brent, 2D NR Noble, and 1D Newman schemes, however, require typically an order of magnitude more EOS calls to converge throughout the entire parameter space (in the range of 500–1000 EOS calls). The color bar is adjusted to enhance the dynamic range for the more efficient schemes. The average number of EOS calls of all test points is summarized in Table 1. These large numbers of EOS calls are due to the fact that the 1D Brent, 2D NR Noble, and 1D Newman schemes each require an additional inversion in the EOS table per solver iteration. Each of these inversions, in turn, typically requires tens of steps to converge and a few EOS calls per step. For the Helmholtz EOS the results are similar to LS220, albeit the number of required EOS calls are lower in absolute numbers for the 1D Brent, 2D NR Noble, and 1D Newman schemes.
Fig. 2 shows the relative error between original and recovered primitive variables after the schemes have converged for the five different schemes and three EOS presented in Fig. 1. The NR schemes generally perform better in this metric as one or two additional iterations of the solver improve the accuracy by orders of magnitude.
Table 1 summarizes the different metrics (number of iterations, number of EOS calls, and relative accuracy) averaged over all points in the test parameter space for the five recovery schemes and three EOS tested here. The averaged values well capture the trends evident from Figs. 1 and 2. Judged by the number of iterations the 1D Brent, 2D Noble and Newman methods are the more efficient schemes. However, the opposite is true for the number of EOS calls, in which case these three schemes need up to one or two orders of magnitude more EOS calls on average than the 2D and 3D NR schemes. For the Helmholtz EOS the situation is not as dire as for LS220, and the number of EOS calls doubles relative to the 2D and 3D NR schemes with the 1D Brent, 2D Noble, and Newman methods. The NR and Newman schemes reach a higher overall relative accuracy than the 1D Brent method.
A comparison of the five different recovery schemes covering the parameter space in Lorentz factor and the ratio of magnetic over fluid pressure is shown in Fig. 3 (analogous to the top panels of Fig. 1). Here we choose a representative temperature of () and a density of . In general, all schemes fail for very strongly magnetized fluid and/or high Lorentz factors. For the ideal gas EOS, the methods of 1D Brent, 2D NR Noble, and 1D Newman recover most of the parameter space well and with few iterations (). However, the 2D and 3D NR methods do not recover the correct solution beyond Lorentz factors of 10. For the LS220 EOS, the situation is similar to the ideal gas EOS with the exception of the 3D NR method recovering points up to high Lorentz factors of . For the Helmholtz EOS, similar results as for the LS220 are obtained, with the 3D NR solver now even recovering up to higher Lorentz factors of . All metrics are summarized in Tab. LABEL:tab:summary_Wp, which shows similar trends as Tab. 1.
Robustness of the schemes involving a Newton-Raphson algorithm (2D NR, 3D NR, 2D NR Noble) is severely limited by the dependence on good initial guesses. These schemes are particularly sensitive to a good initial guess for the velocity components in the regime of high Lorentz factors. Perturbations of the level of 5% in the primitives as considered here cause definite failure for the 2D NR scheme above Lorentz factors of for all EOS, partial failure in this regime for the 2D Noble scheme and the LS220 and Helmholtz EOS, and failure for the 3D NR scheme and the ideal gas EOS as well as drastically increased number of iterations in the case of the LS220 and the Helmholtz EOS (cf. Fig. 3). Additionally, all schemes centered around a Newton-Raphson algorithm directly involve thermodynamic derivatives, which often must be numerically computed for tabulated EOS. These derivatives may be thermodynamically inconsistent (i.e., the Maxwell relations may not be satisfied to machine precision) and they may be noisy, which adds additional sources of error or failure that cannot be controlled and may appear quasi-randomly.
Robustness is significantly improved in the 1D Brent and Newman schemes, which do neither explicitly depend on initial guesses nor on thermodynamic derivatives. In particular, the fact that there are known brackets for the 1D Brent scheme guarantees the existence of a root and convergence of the scheme. However, for 3-parameter EOS, additional inversions from specific internal energy or enthalpy to temperature are required at every iteration step of the respective main scheme, which is typically accomplished via a 1D Newton-Raphson algorithm. This then implicitly introduces dependencies on initial guesses and thermodynamic derivatives. As we find here, inversions of specific internal energy (as in the 1D Brent method) can typically be more robust than inversions of the specific enthalpy (as in the 1D Newman scheme; see the above discussion related to Fig. 1). None of the schemes under consideration here is thus entirely independent of initial guesses. However, in contrast to the 1D Brent method, the 1D Newman scheme starts with a well defined initial guess for the temperature (the minimum temperature for the EOS at given and , see Sec. 3.4.2).
5 Discussion and Conclusion
We performed a comparison of a number of conservative to primitive variable conversion schemes used in different GRMHD codes (Sec. 3), employing analytic, partially tabulated and fully tabulated microphysical EOS. We have analyzed these methods using a well-defined uniform standalone testbed (Sec. 4.1), covering 2D parameter spaces in either density and temperature or Lorentz factor and ratio of magnetic to fluid pressure typically encountered in astrophysical simulations. The schemes are assessed in terms of the three criteria speed (according to number of iterations and number of EOS calls required for convergence), accuracy, and robustness (Sec. 4.1).
While for the ideal gas EOS there is not too much variation among the different schemes in terms of speed and accuracy (typically iterations and EOS calls), there are significant differences for the microphysical EOS (LS220 and Helmholtz EOS; see Tables 1 and LABEL:tab:summary_Wp for a summary). The 2D and 3D NR schemes require 10–20 iterations more to converge in some parts of the parameter space for microphysical EOS relative to other schemes, but require tens to hundreds of EOS calls less and are thus overall much faster. The additional EOS calls in the 1D Brent, 2D NR Noble and 1D Newman schemes arise because of an additional temperature inversion at every solver iteration that requires roughly 10–100 EOS calls by itself (depending on the inversion scheme and parameters). For tabulated EOS like the LS220 and Helmholtz EOS, the computational cost of a single EOS call can be significant. Therefore, employing the most efficient scheme (3D NR) can yield a speed-up by a factor of at least a few to ten, potentially up to a factor of hundred compared to one of the least efficient schemes. Since the cost of the primitive variable recovery in an evolution code can be up to of the total simulation cost, this can have significant performance impacts on the evolution code, and may decide on whether or not a certain simulation is computationally feasible.
The 3D NR scheme clearly emerges as the fastest and most accurate one; however, it is not the most robust scheme. All schemes based on a Newton-Raphson algorithm (2D NR, 3D NR, 2D NR Noble) require good initial guesses and thus do not guarantee convergence. In the case of the 2D NR Noble scheme this can limit the ability to recover at high Lorentz factors (see Fig. 3). When used in an evolution code, these schemes thus depend on the previous time step, which represents a severe limitation that can lead to failure, as data from the previous time step can be corrupted due to numerical error. Furthermore, these schemes explicitly depend on thermodynamic derivatives, which often need to be computed numerically and can thus be noisy and thermodynamically inconsistent. This may lead to failure. The 1D Brent and 1D Newman scheme, however, are both guaranteed to converge and do not explicitly depend on initial guesses and thermodynamic derivatives. However, for microphysical EOS, additional inversions involving the EOS are required at every solver step, which implicitly introduces dependence on initial guesses (and thermodynamic derivatives, depending on the inversion method). Only the 1D Newman method starts from an independently defined initial guess and is thus entirely independent of the previous time step when used in an evolution code. Furthermore, this additional inversion step may fail and thus limit the robustness of the main scheme. Enthalpy inversions as in the 1D Newman scheme tend to be less robust and more prone to failure (typically on the percent level in the density-temperature parameter space) than inversions of the specific internal energy as in the 1D Brent scheme.
In conclusion, the optimal recovery scheme for a GRMHD evolution code depends on the EOS. For an analytic EOS, we suggest to employ a more robust scheme such as the 1D Brent or 1D Newman scheme, with the 1D Newman scheme being preferred as it tends to achieve higher accuracy. For microphysical (3-parameter) EOS, we envision a combination of different schemes. In a first attempt, we suggest to employ the 3D NR scheme, which is both the most efficient and accurate scheme among those tested here. For points that fail in this initial pass, one can then switch to a more robust scheme (1D Brent, 1D Newman) in a second step and exploit the higher robustness in the more extreme regions of the parameter space.
- Aasi et al. (2014) Aasi, J., Abbott, B. P., Abbott, R., et al. 2014, Phys. Rev. Lett., 113, 011102
- Anderson et al. (2008) Anderson, M., Hirschmann, E. W., Lehner, L., et al. 2008, Phys. Rev. Lett., 100, 191101
- Antón et al. (2006) Antón, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296
- Arnowitt et al. (2008) Arnowitt, R., Deser, S., & Misner, C. W. 2008, Gen. Rel. Grav., 40, 1997
- Cerdá-Durán et al. (2008) Cerdá-Durán, P., Font, J. A., Antón, L., & Müller, E. 2008, A&A, 492, 937
- Del Zanna et al. (2007) Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11
- Duez et al. (2005) Duez, M. D., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2005, Phys. Rev. D, 72, 024028, astro-ph/0503420
- Etienne et al. (2015) Etienne, Z. B., Paschalidis, V., Haas, R., Mösta, P., & Shapiro, S. L. 2015, Class. Quant. Grav., 32, 175009
- Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
- Giacomazzo & Rezzolla (2007) Giacomazzo, B., & Rezzolla, L. 2007, Class. Quant. Grav., 24, 235
- Kiuchi et al. (2012) Kiuchi, K., Kyutoku, K., & Shibata, M. 2012, Phys. Rev. D, 86, 064008
- Komissarov (2005) Komissarov, S. S. 2005, MNRAS, 359, 801
- Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys. A, 535, 331
- Lichnerowicz (1944) Lichnerowicz, A. 1944, J. Math. Pures et Appl., 23, 37
- Mösta et al. (2015) Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376
- Mösta et al. (2014a) Mösta, P., Mundim, B. C., Faber, J. A., et al. 2014a, Class. Quant. Grav., 31, 015005
- Mösta et al. (2014b) Mösta, P., Richers, S., Ott, C. D., et al. 2014b, ApJ, 785, L29
- Muhlberger et al. (2014) Muhlberger, C. D., Nouri, F. H., Duez, M. D., et al. 2014, Phys. Rev. D, 90, 104014
- Neilsen et al. (2014) Neilsen, D., Liebling, S. L., Anderson, M., et al. 2014, Phys. Rev. D, 89, 104029
- Newman & Hamlin (2014) Newman, W. I., & Hamlin, N. D. 2014, SIAM J. Sci. Comput., 36
- Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
- Palenzuela et al. (2015) Palenzuela, C., Liebling, S. L., Neilsen, D., et al. 2015, Phys. Rev. D, 92, 044045
- Siegel & Metzger (2017a) Siegel, D. M., & Metzger, B. D. 2017a, Phys. Rev. Lett., 119, 231102
- Siegel & Metzger (2017b) —. 2017b, ApJ submitted, arXiv:1711.00868
- Timmes & Arnett (1999) Timmes, F. X., & Arnett, D. 1999, ApJS, 125, 277
- Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501
- White et al. (2016) White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22