August 3 2018

Abstract

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

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

$$f={f}_{0}+u\left(r-vt\right)\u2215r\text{,}$$ | (1) |

in its diﬀerential form

$${\partial}_{t}f+v{\partial}_{r}f-v\phantom{\rule{0.3em}{0ex}}\left(f-{f}_{0}\right)=0\text{,}$$ | (2) |

where $v$ is the asymptotic propagation speed of the ﬁeld $f$ and ${f}_{0}$ is its asymptotic value.

Quoting [1]:

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

$$f={f}_{0}+u\left(r-vt\right)\u2215r.$$ | (3) |

Where ${f}_{0}$ 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 suﬃciently far away one can safely assume that the speed of light is 1, so $v=1$ for most ﬁelds. However, the gauge variables can easily propagate with a diﬀerent speed implying a diﬀerent value of $v$.

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

$${\partial}_{t}f+v{\partial}_{r}f-v\phantom{\rule{0.3em}{0ex}}\left(f-{f}_{0}\right)\u2215r=0.$$ | (4) |

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

$$\frac{{x}_{i}}{r}\phantom{\rule{0.3em}{0ex}}{\partial}_{t}f+v{\partial}_{i}f+\frac{v{x}_{i}}{{r}^{2}}\phantom{\rule{0.3em}{0ex}}\left(f-{f}_{0}\right)=0.$$ | (5) |

We ﬁnite diﬀerence this condition consistently to second order in both space and time and apply it to all dynamic variables (with possible diﬀerent values of ${f}_{0}$ and $v$) at all boundaries.

There is a ﬁnal subtlety in our boundary treatment. Wave propagation is not the only reason why ﬁelds 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 eﬀect. In order to minimize the error at our boundaries introduced by such non-wavelike evolution, we allow for boundary behavior of the form:

$$f={f}_{0}+u\left(r-vt\right)\u2215r+h\left(t\right)\u2215{r}^{n},$$ | (6) |

with $h$ a function of $t$ alone and $n$ some unknown power. This leads to the diﬀerential equation

$$\begin{array}{rcll}{\partial}_{t}f+v{\partial}_{r}f-\frac{v}{r}\phantom{\rule{0.3em}{0ex}}\left(f-{f}_{0}\right)& =& \frac{vh\left(t\right)}{{r}^{n+1}}\left(1-nv\right)+\frac{{h}^{\prime}\left(t\right)}{{r}^{n}}& \text{}\\ & \simeq & \frac{{h}^{\prime}\left(t\right)}{{r}^{n}}\phantom{\rule{2em}{0ex}}for\phantom{\rule{1em}{0ex}}large\phantom{\rule{1em}{0ex}}r,& \text{(7)}\text{}\text{}\end{array}$$

or in Cartesian coordinates

This expression still contains the unknown function ${h}^{\prime}\left(t\right)$. Having chosen a value of $n$, one can evaluate the above expression one point away from the boundary to solve for ${h}^{\prime}\left(t\right)$, 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.

The thorn implements to ingredients for the evolution.

- the boundary condition (8),
- a routine to extrapolate the contracted Christoﬀel symbols ${g}^{jk}{\Gamma}_{jk}^{i}$ (Xt1…Xt3 in the code) into the outer boundaries at $t=0$.

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.

The derivatives ${\partial}_{i}f$
appearing in (8) are approximated as 2^{nd} order accurate ﬁnite diﬀerence expressions using centred expressions

$${f}^{\prime}\left(x\right)\approx \frac{f\left(x-\Delta x\right)-f\left(x+\Delta x\right)}{2\Delta x}$$ | (9) |

and one sided expressions

$${f}^{\prime}\left(x\right)\approx \pm \frac{3f\left(x\right)-4f\left(x\mp \Delta x\right)+f\left(x\mp 2\Delta x\right)}{2\Delta x}\text{,}$$ | (10) |

for boundaries around the $\pm $ 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 diﬀerence gives us an idea of the missing part and we extrapolate that to the boundary assuming a power-law decay.

The thorn exports two aliased functions ExtrapolateGammas and NewRad_Apply that should be used in INITIAL and MoL_CalcRHS respectively to extrapolate the contracted Christoﬀel 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)

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 ${f}_{0}$, v0 is $v$ and radpower is the assumed decay exponent $n$.

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

git clone git@bitbucket.org:einsteintoolkit/einsteinevolve.git

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.

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

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

<interface.ccl>

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)

USES FUNCTION NewRad_Apply

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)

USES FUNCTION NewRad_Apply

<schedule.ccl>

SCHEDULE ML_BSSN_NewRad IN MoL_CalcRHS AFTER ML_BSSN_EvolutionInterior

{

LANG: C

# No need to sync here -- the RHS variables are never synced

} "Apply NewRad boundary conditions to RHS"

SCHEDULE ML_BSSN_NewRad IN MoL_CalcRHS AFTER ML_BSSN_EvolutionInterior

{

LANG: C

# No need to sync here -- the RHS variables are never synced

} "Apply NewRad boundary conditions to RHS"

<NewRad.cc>

#include "cctk_Arguments.h"

#include "cctk_Functions.h"

void ML_BSSN_NewRad(CCTK_ARGUMENTS)

{

DECLARE_CCTK_ARGUMENTS;

const int radpower = 2; // for example

NewRad_Apply(cctkGH, gt12, gt12rhs, 0.0, 1.0, radpower);

}

#include "cctk_Arguments.h"

#include "cctk_Functions.h"

void ML_BSSN_NewRad(CCTK_ARGUMENTS)

{

DECLARE_CCTK_ARGUMENTS;

const int radpower = 2; // for example

NewRad_Apply(cctkGH, gt12, gt12rhs, 0.0, 1.0, radpower);

}

Please use the Einstein Toolkit mailing list users@einsteintoolkit.org to report issues and ask for help.

- Initial version of method likely by Miguel Alcubierre
- Current version of code by Erik Schnetter

- Initial version by Roland Haas using text by Miguel, Erik.

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“, http://kranccode.org/.

[3] “The Llama code“, https://llamacode.bitbucket.io/.