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]: \begin {equation} f = f_0 + u(r - v t)/r\mbox {,} \label {eqn:EinsteinEvolve_NewRad_bc} \end {equation} in its differential form \begin {equation} \partial _t f + v \partial _r f - v\,(f - f_0) = 0\mbox {,} \label {eqn:EinsteinEvolve_NewRad_diff_bc} \end {equation} where \(v\) is the asymptotic propagation speed of the field \(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 fields behave as spherical waves, namely we impose the condition \begin {equation} f = f_0 + u(r-vt)/r . \label {eq:radiation} \end {equation} 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 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 (??) as it stands, but rather we use it in differential form: \begin {equation} \partial _t f + v \partial _r f - v \, (f - f_0)/r = 0 . \end {equation} Since our code is written in Cartesian coordinates, we transform the last condition to \begin {equation} \frac {x_i}{r} \, \partial _t f + v \partial _i f + \frac {v x_i}{r^2} \, \left ( f - f_0 \right ) = 0 . \end {equation} 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 \(f_0\) 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: \begin {equation} f = f_0 + u(r-vt)/r + h(t)/r^n , \label {eq:radpower} \end {equation} with \(h\) a function of \(t\) alone and \(n\) some unknown power. This leads to the differential equation
or in Cartesian coordinates \begin {equation} \frac {x_i}{r} \, \partial _t f + v \partial _i f + \frac {v x_i}{r^2} \, \left ( f - f_0 \right ) \simeq \frac {x_i h'(t)}{r^{n+1}} . \label {eqn:EinsteinEvolve_NewRad_num_bc} \end {equation}
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.
The thorn implements to ingredients for the evolution.
the boundary condition (??),
a routine to extrapolate the contracted Christoffel symbols \(g^{jk} \Gamma ^i_{jk}\) (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 (??) are approximated as 2nd order accurate finite difference expressions using centred expressions \begin {equation} f'(x) \approx \frac {f(x-\Delta x) - f(x+\Delta x)}{2\Delta x} \end {equation} and one sided expressions \begin {equation} f'(x) \approx \pm \frac {3 f(x) - 4 f(x\mp \Delta x) + f(x\mp 2\Delta x)}{2\Delta x}\mbox {,} \end {equation} 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 difference 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 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 \(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 (??). 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
<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"
<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);
}
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/.
| z_is_radial | Scope: private | BOOLEAN |
| Description: Consider the local z direction to be the radial one; use with Llama
| ||
| Default: no | ||
Implements:
newrad
Uses header:
GenericFD.h
Provides:
ExtrapolateGammas to
NewRad_Apply to
This section lists all the variables which are assigned storage by thorn EinsteinEvolve/NewRad. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.
NONE