Adding initial data

From Einstein Toolkit Documentation
Jump to: navigation, search

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.


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.


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 Kruskal coordinates are beautifully described in Sean M. Carroll's "A No-Nonsense Introduction to General Relativity", 2001, at .

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: Media: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
  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: Media:InitialData.F90

CCL Files: Cactus Declarations

In the last step (now that all the hard work is done), we describe to Cactus how this torn is called, what its variables are, what parameters it takes, and when its routines need to be called.


This thorn is called "Kruskal", and it implements an interface that we also call "Kruskal" (for the lack of a better idea).


This thorn requires access to a few other variables, in particular the grid points' coordinates, and the ADMBase variables which it calculates:

INHERITS: ADMBase StaticConformal grid

It also calculates and stores the Kruskal and Schwarzschild coordinates of every grid point for convenience, so we define grid functions for these:

CCTK_REAL Schwarzschild_radius TYPE=gf
} "Schwarzschild radius"

The complete file is here: Media:interface.ccl.txt


There are two kinds of parameters that we have to deal with:

  • Parameters that choose which kind of initial data should be set up (which we extend to include "Kruskal")
  • Parameters for ourselves, where e.g. the ADM mass of the Schwarzschild spacetime is set, and the slicing is chosen

Let's begin with the initial data choice parameters, which are defined in ADMBase. We obtain access to them, and extend their range (the set of allowed values) to include "Kruskal"). We can do this because ADMBase allows us to do so -- private parameters could not be extended.


EXTENDS KEYWORD initial_lapse
  "Kruskal" :: "Kruskal lapse"

EXTENDS KEYWORD initial_shift
  "Kruskal" :: "Kruskal shift"

EXTENDS KEYWORD initial_data
  "Kruskal" :: "Kruskal metric"

USES KEYWORD metric_type

We also set a few parameters of our own:


BOOLEAN calculate_Schwarzschild_radius "Calculate the Schwarzschild radius?"
} no

CCTK_REAL mass "Black hole mass"
  0:* :: ""
} 1.0

KEYWORD slicing "Choice of slicing"
  "Kruskal"       :: "v = t"
  "Schwarzschild" :: "v = u tanh (t / 4M)"
  "Wald-Iyer"     :: "v = t cos theta"
} "Kruskal"

CCTK_REAL delta "Finite differencing distance"
  (0:* :: ""
} 1.0e-6

The complete file is here: Media:param.ccl.txt


Finally we declare which variables need to have storage, and which routines need to be called when.

We allow users to store the Schwarzschild and Kruskal coordinates of each grid point if the parameter calculate_Schwarzschild_radius is set, so we need to allocate storage for this group:

if (calculate_Schwarzschild_radius)
  STORAGE: Schwarzschild_radius

We also describe under which circumstances our initial data routine should be called, and when it should be called:

if (CCTK_EQUALS (initial_lapse, "Kruskal")
    || CCTK_EQUALS (initial_shift, "Kruskal")
    || CCTK_EQUALS (initial_data, "Kruskal"))
  SCHEDULE Kruskal_InitialData IN ADMBase_InitialData
    LANG: Fortran
  } "Set up Kruskal initial data"

ADMBase provides the schedule group "ADMBase_InitialData" for, well, setting up initial data. This ensures that e.g. coordinates are set up beforehand, and boundary conditions or AMR prolongation or restriction occur afterwards. There is a similar HydroBase schedule group for hydrodynamics variables.

Scheduling the initial data routine takes no options, leaving all schedule options at their default values. This means e.g. that this routine is called once for each refinement level (corresponding to "OPTIONS=local"). Also, no grid functions need to be synchronised since all grid points are set, including boundary and ghost points.

The complete file is here: Media:schedule.ccl.txt

Setting up and Running an Example

You can obtain a tarball of the complete thorn here: Media:Kruskal.tar.gz.

This section is missing. There are example parameter files in the par subdirectory. It is particularly interesting to run "kruskal-wald-iyer.par", which shows that "numerical apparent horizons" are found after all, although true apparent horizons do not exist. Speak about .