Adding initial data

Writing an Initial Data thorn for the Einstein Toolkit

This page explains how to write a new initial data thorn for the Einstein Toolkit. These initial data will set up ADMBase or HydroBase variables, so that they are usable for all other thorns that are based Einstein Toolkit.

Introduction

You need:

• a prescription or routine that actually generates the initial data,
• to decide whether this are initial data for the spacetime, for hydrodynamics, or both,
• these instructions.

Example

Let us assume, for the sake of simplicity, that we are going to set up initial data for the spacetime only. Let us pick a concrete example: the Kruskal coordinates. The Kruskal coordinates describe a single, static black hole (the Schwarzschild spacetime), and their advantage over many other coordinate systems is that they cover the whole, extended Schwarzschild spacetime (i.e. including both asymptotically flat ends, the worm hole, the black and white hole horizons), and that they do not introduce coordinate singularities. Of course there are curvature singularities inside the white and black holes, but there are no "artificial" singularities, e.g. at the horizons. The Krukal coordinates are beautifully described in Sean M. Carroll's "A No-Nonsense Introduction to General Relativity", 2001, at http://pancake.uchicago.edu/~carroll/notes/grtinypdf.pdf .

To make things a bit more interesting, we are also going to apply a coordinate transformation to the Kruskal coordinates, so that we can evaluate them on arbitrary slices. Of course that means that the actual coordinates are not the Kruskal coordinates. This means that the slicing can in principle be chosen arbitrarily within the extended Schwarzschild spacetime. Other coordinates (e.g. ingoing Eddington-Finkelstein coordinates) are regular only a in part of the spacetime. This will allow us to examine arbitrary slicings of Schwarzschild -- even slicings that are not spherically symmetric -- which is very interesting when one studies trapped surfaces or apparent horizons.

In particular, we want to implement the following slicings:

• Kruskal
• Schwarzschild
• the Wald-Iyer slicing introduced in Phys. Rev. D 44, R3719 (1991), "Trapped surfaces in the Schwarzschild geometry and cosmic censorship"

As described in this paper, this slicing has the property that it comes arbitrarily close to the (future) curvature singularity in the Schwarzschild spacetime, but does not contain any apparent horizons. We want to examine this numerically.

Kruskal Solution

As mentioned above, we begin by implementing a Fortran module that evaluates the Kruskal solution at a particular event (t,x,y,z). How to do this is outside the scope of this tutorial; one example implementation is here: File:kruskal.F90.

This module contains a routine that has as its input arguments 4 real numbers t, x, y, and z, describing the spacetime event where the four-metric needs to be evaluated. Its output arguments are alpha, beta^i, and g_ij, which this routine has to fill in. For convenience, this routine also sets the two arguments u and v to the Kruskal coordinates of this event, and two other arguments st and sr to its Schwarzschild coordinates. (Since the metric is spherically symmetric, angular coordinates can be omitted.)

Internally, the routine transforms the tuple (t,x,y,z) to (u,v) and (st,sr) according to the chosen slicing. In particular, the different slicings lead to:

• Kruskal: v = t
• Schwarzschild: v = u tanh (t / 4M)
• Wald-Iyer: v = t z/r

where M is the ADM mass and r is the Pythagorean of (x,y,z). Details do not matter for the tutorial here can be found in the source code. It should be mentioned that this coordinate transformation requires a non-linear root finding which is implemented with the help of the GSL (GNU Scientific Library).

This module is the heart of this initial data thorn, and most of the physical and numerical work has been invested here. What remains it to wrap this routine, so that the Einstein Toolkit understands it.

In many other cases, the initial data will be calculated not pointwise, but by an external solver on an auxiliary grid, and there will be an interpolation routine that evaluates the solution at arbitrary coordinates (t,x,y,z). It is sufficient if the solution is only known at t=0.

Loop over all grid points

The next step is to populate the ADMBase (and/or HydroBase) grid functions. This requires a routine which is scheduled in Cactus's INITIAL bin, and which iterates over all grid points, evaluating the solution there. The solution is evaluated by calling the Kruskal module described above.

Basic Structure

The boilerplate form of such a routine is as follows:

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"

subroutine Kruskal_InitialData (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  integer   :: i, j, k
  do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
           ! DO THE ACTUAL WORK
        end do
     end do
  end do
end subroutine Kruskal_InitialData

This routine has to fill in all ADMBase variables, which are:

• lapse alpha
• shift beta^i
• metric g_ij
• extrinsic curvature K_ij
• potentially: time derivative of lapse
• potentially: time derivative of shift

Whether the latter two are necessary depends on whether the evolution system uses them; for example, the BSSN does use the time derivative of the shift. If they are not set, many thorns will assume that they are zero, which may or may not be what is intended.

Excursion: Extrinsic Curvature

Reviewing what we have so far, we notice that the Kruskal module does not calculate the extrinsic curvature K_ij. The reason is that calculating K_ij requires calculating the time derivative of g_ij, and the coordinate transformation implemented in the Kruskal module is already quite complicated. We are therefore going to calculate K_ij via finite differencing. Note that many initial data prescriptions have explicit expressions for K_ij, so that this step may not be necessary.

This finite differencing is implemented, at each grid point, by evaluating the initial data several times with small offsets in the time and space directions. Straight-forward second order finite differencing then yields all the derivatives of lapse, shift, and three-metric, and the extrinsic curvature is then calculated in a straightforward manner via the ADM relation

d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k

The time derivatives of lapse and shift could be calculated in a similar manner. This is not implemented in this thorn, because its solution is not intended for a time evolution; instead, horizons can be located directly on the initial data.

Setting only the desired quantities

The Einstein Toolkit allows mix-and-matching between different three-metrics, lapse, and shift conditions. (This is only rarely used these days.) Therefore, this routine must set metric, lapse, and shift only when these are actually selected by the user in the parameter file; otherwise, it would overwrite data set by other thorns.

The user chooses initial conditions via the ADMBase parameters

• ADMBase::initial_data (metric and extrinsic curvature)
• ADMBase::initial_lapse
• ADMBase::initial_shift

Each thorn providing initial data needs to extend the range of these parameters (see below). It may only set up its initial data if the user chose so. In our case, we examine these parameters in the beginning

  logical   :: want_lapse, want_shift, want_data
  want_lapse = CCTK_EQUALS (initial_lapse, "Kruskal")
  want_shift = CCTK_EQUALS (initial_shift, "Kruskal")
  want_data  = CCTK_EQUALS (initial_data,  "Kruskal")

and set the grid functions in the loop only if these conditions hold:

           ! Store initial data
           
           if (want_data) then
              gxx(i,j,k) = gama(1,1)
              gxy(i,j,k) = gama(1,2)
              gxz(i,j,k) = gama(1,3)
              gyy(i,j,k) = gama(2,2)
              gyz(i,j,k) = gama(2,3)
              gzz(i,j,k) = gama(3,3)
              
              kxx(i,j,k) = kk(1,1)
              kxy(i,j,k) = kk(1,2)
              kxz(i,j,k) = kk(1,3)
              kyy(i,j,k) = kk(2,2)
              kyz(i,j,k) = kk(2,3)
              kzz(i,j,k) = kk(3,3)
           end if
           
           if (want_lapse) then
              alp(i,j,k) = alfa
           end if
           
           if (want_shift) then
              betax(i,j,k) = beta(1)
              betay(i,j,k) = beta(2)
              betaz(i,j,k) = beta(3)
           end if

For reference, the complete initial data routine is here: File:InitialData.F90