Roland Haas <>

August 3 2018


This thorn implements the radiative boundary condition originally introduced by the BSSN_MoL thorn.

1 Introduction

This thorn provides routines to implement radiative boundary conditions as described in section VI of [1]:

f = f0 + u(r vt)r, (1)

in its differential form

tf + vrf v(f f0) = 0, (2)

where v is the asymptotic propagation speed of the field f and f0 is its asymptotic value.

2 Physical System

Quoting [1]:

At the outer boundary we use a radiation (Sommerfeld) boundary condition. We start from the assumption that near the boundary all fields behave as spherical waves, namely we impose the condition

f = f0 + u(r vt)r. (3)

Where f0 is the asymptotic value of a given dynamical variable (typically 1 for the lapse and diagonal metric components, and zero for everything else), and v is some wave speed. If our boundary is sufficiently far away one can safely assume that the speed of light is 1, so v = 1 for most fields. However, the gauge variables can easily propagate with a different speed implying a different value of v.

In practice, we do not use the boundary condition (3) as it stands, but rather we use it in differential form:

tf + vrf v(f f0)r = 0. (4)

Since our code is written in Cartesian coordinates, we transform the last condition to

xi r tf + vif + vxi r2 f f0 = 0. (5)

We finite difference this condition consistently to second order in both space and time and apply it to all dynamic variables (with possible different values of f0 and v) at all boundaries.

There is a final subtlety in our boundary treatment. Wave propagation is not the only reason why fields evolve near a boundary. Simple infall of the coordinate observers will cause some small evolution as well, and such evolution is poorly modeled by a propagating wave. This is particularly important at early times, when the radiative boundary condition introduces a bad transient effect. In order to minimize the error at our boundaries introduced by such non-wavelike evolution, we allow for boundary behavior of the form:

f = f0 + u(r vt)r + h(t)rn, (6)

with h a function of t alone and n some unknown power. This leads to the differential equation

tf + vrf v r(f f0) = vh(t) rn+1 1 nv + h(t) rn h(t) rn forlarger, (7)

or in Cartesian coordinates

xi r tf + vif + vxi r2 f f0 xih(t) rn+1 . (8)

This expression still contains the unknown function h(t). Having chosen a value of n, one can evaluate the above expression one point away from the boundary to solve for h(t), and then use this value at the boundary itself. Empirically, we have found that taking n = 3 almost completely eliminates the bad transient caused by the radiative boundary condition on its own.

3 Numerical Implementation

The thorn implements to ingredients for the evolution.

  1. the boundary condition (8),
  2. a routine to extrapolate the contracted Christoffel symbols gjkΓjki (Xt1Xt3 in the code) into the outer boundaries at t = 0.

3.1 Extrapolation

At outer boundary points, as determined by GenericFD’s routine GenericFD_GetBoundaryInfo we use cubic polynomials along the normal of the boundary to extrapolate data from the interior into the boundary region. Due to the use of cubic polynomial we require at least four (4) points of valid data to be available for the extrapolation. This includes ghost zones, which are assumed to contain valid data in the interior of the domain.

3.2 Radiative boundary condition

The derivatives if appearing in (8) are approximated as 2nd order accurate finite difference expressions using centred expressions

f(x) f(x Δx) f(x + Δx) 2Δx (9)

and one sided expressions

f(x) ±3f(x) 4f(x Δx) + f(x 2Δx) 2Δx , (10)

for boundaries around the ± coordinate direction.

Finally we try to extrapolate for the part of the boundary that does not behave as a pure wave (i.e. Coulomb type terms caused by infall of the coordinate lines).

This we do by comparing the source term one grid point away from the boundary (which we already have), to what we would have obtained if we had used the boundary condition there. The difference gives us an idea of the missing part and we extrapolate that to the boundary assuming a power-law decay.

4 Using This Thorn

The thorn exports two aliased functions ExtrapolateGammas and NewRad_Apply that should be used in INITIAL and MoL_CalcRHS respectively to extrapolate the contracted Christoffel symbols and apply radiative boundary conditions to evolved variables respectively.

There call signatures are:

CCTK_INT FUNCTION                         \  
    ExtrapolateGammas                     \  
        (CCTK_POINTER_TO_CONST IN cctkGH, \  
         CCTK_REAL ARRAY INOUT var)  
CCTK_INT FUNCTION                         \  
    NewRad_Apply                          \  
        (CCTK_POINTER_TO_CONST IN cctkGH, \  
         CCTK_REAL ARRAY IN var,          \  
         CCTK_REAL ARRAY INOUT rhs,       \  
         CCTK_REAL IN var0,               \  
         CCTK_REAL IN v0,                 \  
         CCTK_INT IN radpower)

where var0 corresponds to f0, v0 is v and radpower is the assumed decay exponent n.

4.1 Obtaining This Thorn

The thorn is part of the EinsteinEvolve arrangement available on bitbucket.

git clone

4.2 Basic Usage

4.3 Special Behaviour

If the parameter z_is_radial is set it assumes that the z direction is a radial direction and uses this to evaluate the radial derivative in (8). This is the case for the multi-patch coordinate systems provided by the Llama [3] infrastructure.

4.4 Interaction With Other Thorns

This thorn requires thorn GenericFD which is part of Kranc [2].

4.5 Examples

The best example is likely to inspect the source code ML_BSSN and ML_BSSN_Helper, namely:

CCTK_INT FUNCTION                         \  
    NewRad_Apply                          \  
        (CCTK_POINTER_TO_CONST IN cctkGH, \  
         CCTK_REAL ARRAY IN var,          \  
         CCTK_REAL ARRAY INOUT rhs,       \  
         CCTK_REAL IN var0,               \  
         CCTK_REAL IN v0,                 \  
         CCTK_INT IN radpower)  

  LANG: C  
  # No need to sync here -- the RHS variables are never synced  
} "Apply NewRad boundary conditions to RHS"

#include "cctk_Arguments.h"  
#include "cctk_Functions.h"  
  const int radpower = 2; // for example  
  NewRad_Apply(cctkGH, gt12, gt12rhs, 0.0, 1.0, radpower);  

4.6 Support and Feedback

Please use the Einstein Toolkit mailing list to report issues and ask for help.

5 History

5.1 Thorn Source Code

5.2 Thorn Documentation

5.3 Acknowledgements

This documentation copies text from comments in the source code as well as the paper [1].


[1]   M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel and R. Takahashi, “Gauge conditions for long term numerical black hole evolutions without excision,” Phys. Rev. D 67, 084023 (2003) doi:10.1103/PhysRevD.67.084023 [gr-qc/0206072].

[2]   “Kranc: Kranc assembles numerical code“,

[3]   “The Llama code“,