Adding custom analytic spacetimes

Test case

We will show how to add an analytic spactime by implementing the Ellis wormhole metric. In order to make it slightly less trivial, and since it is a good practice in spherical geometry, we will us a logarithmic radial coordinate (or actually, a similar transformation that allows us to cover the whole spacetime with a single chart).

The metric is:

where $\rho \in (-\infty,\infty)$, $n$ is a parameter that sets the scale of the throat, and $r=+\sqrt{\rho^2 + n^2}$ is the areal radius.

The coordinate transformation will be similar to that in Combi et al. (2024),

where $c_1$ and $c_2$ are constants. We will set $c_1=2, c_2=1$.

In terms of the new coordinate $\eta$, the metric is:

The 3+1 decomposition of this metric is straightforward:

$\alpha=1$, $\beta^i=0$, and

$\gamma_{ij}={\rm diag}\{4 \cosh^2 \eta, (4\sinh^2\eta + n^2),(4\sinh^2\eta + n^2)\sin^2\theta)\}$.

Implementation

We start by creating a new coordinates file. The easiest way is by copying an existent one, ideally, the one closest to the one we wish to implement.
In this case, it is probably that of logarithmic spherical coordinates, mod_coord_splog.t, written by Pushpita Das.

Let’s call our new file mod_coord_ellis.t.

cd $BHAC_DIR/src/geometry
cp mod_coord_splog.t mod_coord_ellis.t

Now we open the file and start filling in the general information on the new coordinates. We can call them ‘ellis‘. Since the 3+1 decomposition is trivial, we don’t initialize from the 4-metric.

character*20, parameter              :: coord="ellis"
logical, save                        :: init_from_g4 = .false.

The metric has a single parameter (ncoordpar=1), which is $n$.
We allocate this parameter at index n_=1 of the array coordpar.
The actual value of the parameter can be set elsewhere, ideally, in initglobaldata_usr() in amrvacusr.t.

integer, parameter                   :: ncoordpar = 1, n_ = 1
double precision, save               :: coordpar(ncoordpar)

No we go to the ‘contains‘ of the module from top to bottom.

The subroutine init_coord can be used for any process that needs to be executed before everything else in the module, for example, to set pointers to pre-existent coordinate transformations, or to read a file with a tabulated numerical metric. In this case, we will leave it empty.

For this particular file, the two following functions (rsp(s) and s(rsp)) contain the transformation from areal coordinates to logarithmic radial coordinates and its inverse. We will replace them by the transformation between our $\rho$ (rsp, for ‘radius spherical’) and our $\eta$.

  !===============================
  double precision function eta(rsp)
    use mod_amrvacdef

    double precision :: rsp

    s = asinh(0.5d0*rsp)

  end function s
  !================================
  double precision function rsp(eta)

    double precision :: s

    rsp = 2.d0*sinh(s)

  end function rsp
  !===============================

The subroutine get_sqrtgamma_analytic() can be used to provide an analytic expression for the determinant of the spatial metric. In this case, such expression is easy to find, as the metric is diagonal.
There are two pre-processor options: when the $\theta$ coordinate is present (^IFZIN) and when it is not (^IFZOUT). So we code:

    {^IFZIN
    sqrtgamma = (2.0d0*cosh(x1)) * ((rsp(x1))**2 + coordpar(n_)**2) * abs(sin(x^Z))
    }{^IFZOUT
    sqrtgamma = (2.0d0*cosh(x1)) * ((rsp(x1))**2 + coordpar(n_)**2)
    }

If the expression for the metric determinant is hard to find or too complicated, it is possible to set instead

if(present(is_analytic)) is_analytic = .false.

and set sqrtgamma to any value. In such case, the determinant will be calculated numerically from the spatial metric $\gamma_{ij}$.

In the next three subroutines, get_alpha(), get_beta and get_g_component(), we provide the lapse and the components of the contravariant shift and the covariant spatial metric.

In the case where the metric is initialized from the 4-metric (init_from_g4 = .true.), it will be necessary to fill only get_g_component() providing the $g_{0\mu}$ components, and the lapse and shift will be calculated numerically.

These routines also require to provide the metric derivatives in each direction, which are assigned to the optional output argument dgdk.

There are other two optional arguments, iszero and dgdk_iszero. These are to let the code know when some components (or their derivatives) are zero, so that they point to a single zero pointer in order to avoid wasting memory by storing arrays of zeros.

For our present case, the routines get_alpha() and get_beta remain unchanged, and we only need to change the spatial components of the metric and their derivatives.