% *======================================================================*
%  Cactus Thorn template for ThornGuide documentation
%  Author: Ian Kelley
%  Date: Sun Jun 02, 2002
%  $Header$
%
%  Thorn documentation in the latex file doc/documentation.tex
%  will be included in ThornGuides built with the Cactus make system.
%  The scripts employed by the make system automatically include
%  pages about variables, parameters and scheduling parsed from the
%  relevant thorn CCL files.
%
%  This template contains guidelines which help to assure that your
%  documentation will be correctly added to ThornGuides. More
%  information is available in the Cactus UsersGuide.
%
%  Guidelines:
%   - Do not change anything before the line
%       % START CACTUS THORNGUIDE",
%     except for filling in the title, author, date, etc. fields.
%        - Each of these fields should only be on ONE line.
%        - Author names should be separated with a \\ or a comma.
%   - You can define your own macros, but they must appear after
%     the START CACTUS THORNGUIDE line, and must not redefine standard
%     latex commands.
%   - To avoid name clashes with other thorns, 'labels', 'citations',
%     'references', and 'image' names should conform to the following
%     convention:
%       ARRANGEMENT_THORN_LABEL
%     For example, an image wave.eps in the arrangement CactusWave and
%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
%   - Graphics should only be included using the graphicx package.
%     More specifically, with the "\includegraphics" command.  Do
%     not specify any graphic file extensions in your .tex file. This
%     will allow us to create a PDF version of the ThornGuide
%     via pdflatex.
%   - References should be included with the latex "\bibitem" command.
%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
%   - Do not use \appendix, instead include any appendices you need as
%     standard sections.
%   - For the benefit of our Perl scripts, and for future extensions,
%     please use simple latex.
%
% *======================================================================*
%
% Example of including a graphic image:
%    \begin{figure}[ht]
% 	\begin{center}
%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
% 	\end{center}
% 	\caption{Illustration of this and that}
% 	\label{MyArrangement_MyThorn_MyLabel}
%    \end{figure}
%
% Example of using a label:
%   \label{MyArrangement_MyThorn_MyLabel}
%
% Example of a citation:
%    \cite{MyArrangement_MyThorn_Author99}
%
% Example of including a reference
%   \bibitem{MyArrangement_MyThorn_Author99}
%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
%   1--16. {\tt http://www.nowhere.com/}}
%
% *======================================================================*

% If you are using CVS use this line to give version information
% $Header$

\documentclass{article}

% Use the Cactus ThornGuide style file
% (Automatically used from Cactus distribution, if you have a
%  thorn without the Cactus Flesh download this from the Cactus
%  homepage at www.cactuscode.org)
\usepackage{../../../../doc/latex/cactus}

\begin{document}

% The author of the documentation
\author{Roland Haas {\textless}rhaas@ncsa.illinois.edu{\textgreater}}

% The title of the document (not necessarily the name of the Thorn)
\title{NewRad}

% the date your document was last changed, if your document is in CVS,
% please use:
%    \date{$ $Date$ $}
\date{August 3 2018}

\maketitle

% Do not delete next line
% START CACTUS THORNGUIDE

% Add all definitions used in this documentation here
%   \def\mydef etc

% Add an abstract for this thorn's documentation
\begin{abstract}
This thorn implements the radiative boundary condition originally introduced
by the BSSN\_MoL thorn.
\end{abstract}

% The following sections are suggestive only.
% Remove them or add your own.

\section{Introduction}
This thorn provides routines to implement radiative boundary conditions as
described in section VI of~\cite{EinsteinEvolve_NewRad_Alcubierre:2002kk}:
\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.

\section{Physical System}
% RH: this is copied verbatim from Alcubierre:2002kk
Quoting~\cite{EinsteinEvolve_NewRad_Alcubierre:2002kk}: 

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~(\ref{eq:radiation})
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
\begin{eqnarray}
\partial_t f + v \partial_r f
- \frac{v}{r} \, (f - f_0) &=&
\frac{v h(t)}{r^{n+1}} \left( 1 - n v \right)
+ \frac{h'(t)}{r^n} \nonumber \\
&\simeq& \frac{h'(t)}{r^n} \qquad \mathrm{for\ large\ } r ,
\end{eqnarray}
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.

\section{Numerical Implementation}
The thorn implements to ingredients for the evolution. 
\begin{enumerate}
\item the boundary condition~\eqref{eqn:EinsteinEvolve_NewRad_num_bc},
\item a routine to extrapolate the contracted Christoffel symbols
$g^{jk} \Gamma^i_{jk}$ (\code{Xt1}\ldots\code{Xt3} in the code) into the
outer boundaries at $t=0$.
\end{enumerate}

\subsection{Extrapolation}
At outer boundary points, as determined by \code{GenericFD}'s routine
\code{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 \emph{includes} ghost zones, which are assumed to contain valid
data in the interior of the domain. 

\subsection{Radiative boundary condition}
The derivatives $\partial_i f$ appearing in
\eqref{eqn:EinsteinEvolve_NewRad_num_bc} are approximated as
2\textsuperscript{nd} 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\mp2\Delta
x)}{2\Delta x}\mbox{,}
\end{equation}
for boundaries around the $\pm$ coordinate direction. 

% copied from newrad.cc so likely originally written by Miguel
% Alcubierre.
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.

\section{Using This Thorn}
The thorn exports two aliased functions \code{ExtrapolateGammas} and
\code{NewRad\_Apply} that should be used in \code{INITIAL} and
\code{MoL\_CalcRHS} respectively to extrapolate the contracted
Christoffel symbols and apply radiative boundary conditions to evolved
variables respectively.

There call signatures are:
\begin{verbatim}
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)
\end{verbatim}
where \code{var0} corresponds to $f_0$, \code{v0} is $v$ and
\code{radpower} is the assumed decay exponent $n$.

\subsection{Obtaining This Thorn}
The thorn is part of the EinsteinEvolve arrangement available on
bitbucket.
\begin{verbatim}
git clone git@bitbucket.org:einsteintoolkit/einsteinevolve.git
\end{verbatim}

\subsection{Basic Usage}

\subsection{Special Behaviour}
If the parameter \code{z\_is\_radial} is set it assumes that the $z$
direction is a radial direction and uses this to evaluate the radial
derivative in \eqref{eqn:EinsteinEvolve_NewRad_num_bc}. This is the case
for the multi-patch coordinate systems provided by the
\code{Llama}~\cite{EinsteinEvolve_NewRad_llamacode} infrastructure.

\subsection{Interaction With Other Thorns}
This thorn requires thorn \code{GenericFD} which is part of
\code{Kranc}~\cite{EinsteinEvolve_NewRad_kranccode}.

\subsection{Examples}
The best example is likely to inspect the source code \code{ML\_BSSN}
and \code{ML\_BSSN\_Helper}, namely:

\begin{verbatim}
<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
\end{verbatim}

\begin{verbatim}
<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"
\end{verbatim}

\begin{verbatim}
<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);
}
\end{verbatim}

\subsection{Support and Feedback}
Please use the Einstein Toolkit mailing list
\url{users@einsteintoolkit.org} to report issues and ask for help.

\section{History}

\subsection{Thorn Source Code}
\begin{itemize}
\item Initial version of method likely by Miguel Alcubierre
\item Current version of code by Erik Schnetter
\end{itemize}

\subsection{Thorn Documentation}
\begin{itemize}
\item Initial version by Roland Haas using text by Miguel, Erik.
\end{itemize}

\subsection{Acknowledgements}
This documentation copies text from comments in the source code as well
as the paper~\cite{EinsteinEvolve_NewRad_Alcubierre:2002kk}.

\begin{thebibliography}{9}
\bibitem{EinsteinEvolve_NewRad_Alcubierre:2002kk} 
  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 {\bf 67}, 084023 (2003)
  doi:10.1103/PhysRevD.67.084023
  [gr-qc/0206072].
\bibitem{EinsteinEvolve_NewRad_kranccode} ``Kranc: Kranc assembles
numerical code``, \url{http://kranccode.org/}.
\bibitem{EinsteinEvolve_NewRad_llamacode} ``The Llama code``, \url{https://llamacode.bitbucket.io/}.

\end{thebibliography}

% Do not delete next line
% END CACTUS THORNGUIDE

\end{document}
