Overview
The modular structure of BHAC allows to add new analytic coordinates or read numerical ones in a straightforward way.
While reading a numerical metric only requires preparing the tabulated metric in a specific format (see below), the implementation of a new analytic metric requires the addition of a new ‘coordinates’ file inside the directory $BHAC_DIR/src/geometry . These files have names similar to mod_coord_bl.t or mod_coord_cks.t, where ‘bl’ stands for Boyer-Lindquist and ‘cks’ stands for Cartesian Kerr-Schild coordinates. The suffix given to the file will be used to indicate the code to use that specific file by running the preprocessor as e.g.
$> $BHAC_DIR/setup.pl -coord=cks
before compilation. For a recent list of the implemented metrics and the meaning of suffixes you may check the list in Table 1 of Olivares et al. 2019, or look at the header information in each of the coordinates files.
Although the code solves the GRMHD equations in the 3+1 form, it can handle spacetime metrics given either in 3+1 form (i.e., expressed in terms of the shift, lapse and spatial metric) or as the full 4-metric, as the necessary conversions to the 3+1 form are done internally. To indicate the code that it should make the conversion, the flag ‘init_from_g4’ must be set to true.
A set of subroutines get_alpha, get_beta, and get_g_component can be filled by the user with the desired expressions for the metric components. Due to the sparsity of many metrics, a set of flags ‘is_zero’ has been defined in order to save memory by making all of the zero components point to a common zero.
The user has the option to provide analytic expression for other needed metric-related quantities, namely the inverse metric, the square root of the spatial metric determinant and the metric derivatives, although these can also be calculated numerically by the code. For the square root of the metric determinant and the inverse metric, this is controlled by the the flag sqrtgamma_is_analytic and by the function pointer get_gammainv_component_analytic. 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 latter requires defining complex functions for the metric components, as it will be seen below. However, the tedious and error-prone process of transcribing long expressions for several metric derivatives makes worth this additional small effort.
Finally, it might be necessary to add coordinate and vector transformation routines. This is needed in order to visualize simulations using software such as VisIt or Paraview, for which data needs to be given in Cartesian coordinates in order to be correctly plotted. It can also be needed in cases where initial or boundary data is given in a different set of coordinates (e.g. when trying to simulate a Fishbone-Moncrief torus in Cartesian Kerr-Schild coordinates when the initial data is given in Boyer-Lindquist coordinates).
Analytic metric
- Create a new coordinates file. The easiest way to do this is by copying one of the already existing ones, ideally one that describe similar coordinates as those we want to implement. For example, we could copy mod_coord_cart.t if we want to implement a new metric in Cartesian-like coordinates, or mod_coord_mks.t if we want to implement a metric in radially logarithmic Kerr-Schild-like coordinates. The name of the new file must be ‘mod_coord_’Indicate at the header global information on the metric, as e.g. parameters or whether it needs to be initialized from the 4-metric.
- Fill the get_alpha, get_beta and get_g_component subroutines. Explain indices and use ‘is_zero’ flags. If the 4-metric is provided, only get_g_component is needed.
- Do you want to provide an analytic expression for sqrt(gamma)?
- Do you want to provide an analytic expression for the inverse metric?
- Do you want to calculate derivatives numerically?
- Provide transformation routines
- Verify correctness of the metric and metric derivatives
Numerical metric
In order to use a numerical metric, we need to set the coordinates to ‘num’, make clean and make:
$> $BHAC_DIR/setup.pl -coord=num $> make clean && make
To run the code, you need to have an ascii file called ‘numerical.met’ inside your working directory. This file should contain the following information:
Metric name
NEW LINE-------
Dimensions (1, 2, or 3)
Type of coordinates (1=spherical, 2=cylindrical, 3=cartesian*)
Coordinate 1 (r) is present (T/F)?
Coordinate 2 (theta) is present (T/F)?
Coordinate 3 (phi) is present (T/F)?
NEW LINE-------
Is the metric in 3+1 form (T/F)**
NEW LINE-------
Are the metric components 'default' (i.e. same as Schwarzschild with M=1 in Boyer-Lindquist)? T/F in the following order:
alpha, beta1, beta2, beta3, g11, g12, g13, g22, g23, g33
NEW LINE-------
Number of points (Total)
No. of points in direction 1
No. of points in direction 2
No. of points in direction 3
NEW LINE-------
List of values in columns
coordinate 1 (coordinate 2) (coordinate 3) metric components that are not 'default' in the same order as above.
* At the moment, only spherical coordinates are really implemented.
** At the moment, only numerical metrics in the 3+1 are really implementedIn the list of values, you do not need to fill in values for coordinates 2 and 3 if your metric is only 1D, neither values for metric components that are default, since the code already knows what to expect from the information at the header.
As an example, these are the first lines of the metric file of a spherically symmetric metric:
Spherical Boson Star 1D
1 1 T F F
T
F T T T F T T T T T
2050 2050 1 1
-0.57606074380309047 0.29084512426076159 1.0876471378629302
-0.56430440209282329 0.28990767238420562 1.0844094965852946
-0.55254806038255611 0.28898705485181941 1.0812190279524589
In this example, the only coordinate present is r, and the only non-deafault metric functions are alpha and g_{rr}. For this particular application, the first values of the radius are negative. The reason is that the simulation may include the origin, and these first values are meant to be used in the ghost cells, where the metric is also needed. For a black hole simulation, a few points inside the outer event horizon (for horizon-penetrating coordinates) or inside the expected inner boundary of the domain (for non horizon-penetrating coordinates) are sufficient.
As mentioned before, there is no need to introduce information on the metric derivatives, despite they are needed to compute the source terms, since they are calculated numerically by BHAC using second order finite differences. The accuracy of these derivatives depends, however, on the resolution at which the metric is given.
Finally, the resolution of the numerical metric does not need to coincide with that of the GRMHD simulation, and the same is true for the domain size. If the domain of the GRMHD simulation is larger, the default metric (Schwarzschild, M=1) will be applied beyond the end of the numerical one, which is usually not so bad because of Jebsen-Birkhoff’s theorem.