path: root/doc/slice_evolver.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentstyle[a4wide]{article}
\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Documentation for the slice evolver in Cactus thorn Exact}

\author{Carsten Gundlach\\ Faculty for Mathematical Studies\\
University of Southampton\\ Southampton SO17 1BJ, UK}

\date{6 June 02}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This text gives all the formulas that are used in the slice evolver of
Cactus thorn Exact. I wrote it originally for a joint paper with
Miguel Alcubierre and Paul Walker that we never finished. The two last
sections are the original code documentation I wrote.

Here are some ways of using exact solutions in numerical relativity:

1) Use preferred slices of exact solutions to obtain initial data,
lapse and shift. Evolve with this lapse and shift, and with exact
boundary data or some ad hoc BC.

2) Obtain initial data (solutions of the constraints) without
symmetries by using an arbitrary slice in an exact solution. 

3) Evolve these data with an arbitrary lapse and shift, and at the
same time evolve the slice with the same lapse and shift.

a) Test accuracy of the code.

b) Test stability and accuracy of inner boundary conditions.

c) Test effect of outer boundary conditions on introducing
``spurious'' waves.

4) Test gauge conditions on their own, without the need for stable BC
and a stable evolution scheme, by evolving the slice only, and faking
the evolution of the Cauchy data. See if these can smooth a highly
distorted slice of Kerr or Schwarzschild.

The main idea of thorn Exact is that each physically different
spacetime is stored only in one set of coordinates, namely those that
are the standard for it and which are adapted to its symmetries
(``preferred coordinates''). One can then obtain it in any other
coordinate system, in particular one that hides the symmetries, by a
numerical coordinate transformation. As we are interested in a 3+1
split of the metric, and in Cauchy data, the fundamental object is a
slice $X^A(x^i)$ embedded in the spacetime. The four coordinates $X^A$
are the preferred ones, and the three coordinates $x^i$ are the ones
intrinsic to the slice.  From $X^A(x^i)$ and the four-metric $g_{AB}$
in the preferred coordinates, we compute the three-metric and
extrinsic curvature on the slice. We can then either evolve the $g$
and $K$ with the usual methods, or, if we want to test only a gauge
condition, we can evolve the slice itself, and read off $g$ and $K$ in
its new position. 

Since Miguel and I wrote the original code, several physical
spacetimes have been implemented in a number of coordinate
systems. However, if you want an exact solution in completely generic
coordinates, or if you want to test a gauge condition, the slice
evolver is the way to go.

Alas, the \verb|"Minkowski/conf wave"| model doesn't work with the
arbitrary-slice evolver.  There's no warning, you'll just silently
get wrong results.  Sigh\dots

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Extracting Cauchy data on an arbitrary slice}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An arbitrary slice in the exact solution, together with a coordinate
system $x^i$ on this slice is given through four scalar functions
$X^A$ by $x^A=X^A(x^i)$, where $x^A$ are the preferred coordinates of
the exact solution, and $(t,x^i)$ are the coordinates of the numerical
code.  The induced 3-metric on the slice is simply
\begin{equation}
\label{3-metric}
g_{ij} = {X^A}_{,i} {X^B}_{,j} g_{AB},
\end{equation}
where we use a comma to denote partial derivatives. The index or the
context  will always make
clear if these are partial derivatives in the $x^A$ or $(t,x^i)$
coordinates. Next we construct the unit normal vector to the slice
as
\begin{equation}
n^A = N \bar n^A, \qquad N = (\bar n^A \bar n^B g_{AB})^{-1/2}
\end{equation}
where a non-unit normal vector $\bar n^A$ is 
\begin{equation}
\bar n^A = g^{AB} \epsilon_{BCDE} {X^C}_{,1} {X^D}_{,2} {X^E}_{,3}.
\end{equation}
Here $\epsilon_{BCDE}$ is any totally antisymmetric tensor, for
example the one given by $\epsilon_{0123}=1$ in coordinates $X^A$. We
only have to make sure that $\bar n^A$ is future-pointing.

The extrinsic curvature $K_{ij}$ of $\Sigma$ is calculated easily if
we extend the coordinates $x^i$ on $\Sigma$ to geodesic coordinates
$(x^i,t)$ on a neighbourhood of $\Sigma$. In these coordinates,
\begin{equation}
\label{gdot}
K_{ij} = -{1\over 2}{\partial\over\partial t} g_{ij}  \hbox{ (for } 
\alpha=1 \hbox{ and } \beta^i = 0 \hbox{ only).}
\end{equation}
In the same coordinates we also have 
\begin{equation}
{\partial X^A \over \partial t} = n^A \hbox{ (for } 
\alpha=1 \hbox{ and } \beta^i = 0 \hbox{ only).}
\end{equation}
Substituting (\ref{3-metric}) into (\ref{gdot}), and using the trivial
commutation relation
\begin{equation}
\left[{\partial\over\partial t}, {\partial\over\partial
x^i}\right] = 0, 
\end{equation}
we obtain
\begin{equation}
\label{Kij1}
K_{ij} = -{1\over 2}\left({n^A}_{,i} {X^B}_{,j} + {X^A}_{,i} {n^B}_{,j}\right)
g_{AB}  -{1\over 2}{X^A}_{,i} {X^B}_{,j} n^C g_{AB,C}.
\end{equation}
Note that we are now back within a single slice $\Sigma$, and this
result is valid without reference to a particular coordinate system on
spacetime beyond the slice.  This gives us a practical algorithm for
calculating $g_{ij}$ and $K_{ij}$, given the four functions $X^A$ on a
numerical $x^i$ grid. We obtain the ${X^A}_{,i}$ by finite differencing on
the grid, calculate $g_{ij}$ and $n^A$ algebraically for each point on
the grid, finite difference again to obtain ${n^A}_{,i}$, and
calculate $K_{ij}$ again algebraically. The metric coefficients
$g_{AB}(x^C)$ are given in closed form, and are accessed in practice
through a function call. We can obtain the $g_{AB,C}$ in closed form,
but in order to save labor and the scope for error, we actually obtain
them by finite differencing. The error incurred is negligible, as we
are not bound to a numerical grid, but can use a tiny finite
differencing interval $\Delta x^C$.

The term ${n^A}_{,i}$ contains both second derivatives of $X^A$ and
first derivatives of $g_{AB}$. While the former must necessarily be
calculated by finite differencing on the grid, the latter can be
obtained to higher accuracy using the analytic expressions for
$g_{AB,C}$. (This is useful at least when $X^A_{,i}$ are small.)
Therefore we prefer to eliminate ${n^A}_{,i}$ in favor of
${X^A}_{,ij}$ for numerical purposes. We do this by taking the
derivative of the identity
\begin{equation}
n^A {X^B}_{,j} g_{AB} = 0,
\end{equation}
and obtain
\begin{equation}
\label{Kij2}
K_{ij} =  {X^A}_{,ij} n^B g_{AB}  + {1\over 2} \left(
{X^A}_{,i} {X^C}_{,j} n^B + {X^C}_{,i} {X^B}_{,j} n^A  
- {X^A}_{,i} {X^B}_{,j} n^C 
\right)g_{AB,C}.
\end{equation}

In the hyperbolic Bona-Mass\'o formulation of the Einstein equations,
the spatial derivatives of the metric form part of the initial
data. Again we split these up into derivatives of the slicing and
the exact metric:
\begin{equation}
g_{ij,k}= \left({X^A}_{,ik} {X^B}_{,j} 
+ {X^A}_{,jk} {X^B}_{,i} \right) g_{AB}+ 
{X^A}_{,i} {X^B}_{,j} {X^C}_{,k}  g_{AB,C}.
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Evolving a slice in time}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To advance a given slice in time, for a given lapse $\alpha(x^i)$ and shift
$\beta^i(x^j)$, we use the definition of the vector field
$\partial/\partial t$ to obtain
\begin{equation}
{\partial  X^A \over\partial t}= \alpha n^A + \beta^i {X^A}_{,i},
\end{equation}
with $n^A$ defined above in terms of $X^A$ and their first and second
derivatives.  This is a partial differential equation for the four
scalar functions $X^A$. In a numerical relativity situation, we want
to evolve this equation for an interval $\Delta t$ in which $\alpha$
and $\beta^i$ are considered as fixed functions of the coordinates
$x^i$. To compare this scheme with a second order convergent numerical
relativity code, the numerical error must be of $O(\Delta t^3)$. One
way of doing this is using a Runge-Kutta like approach. In this we
evolve forward in time by $(\Delta t)/2$, recompute ${X^A}_{,i}$ and
$n^A$, and leapfrog over this half step. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Using the preferred slicing}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Things simplify if we use the preferred slicing and spatial
coordinates 
\begin{equation}
X^0 = t, \qquad X^I = x^i.
\end{equation}
Then we  have $g_{ij}=g_{IJ}$, and the expression for the extrinsic
curvature simplifies to
\begin{equation}
\label{Kij3}
K_{ij} = {1\over 2\alpha} \left(-g_{IJ,0} + g_{IJ,K} \beta^K
+ {\beta^K}_{,I} g_{KJ} + {\beta^K}_{,J} g_{KI} \right).
\end{equation}
Here we have used the expressions for the lapse and shift in the
preferred slicing
\begin{equation}
\alpha=(-g^{00})^{-1/2}, \qquad \beta^I = -{g^{0I}\over g^{00}}.
\end{equation}
These expressions are useful in their own right, as we may want to use
the analytic value of the lapse and shift in the preferred slicing for
the numerical evolution.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Code documentation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Rudimentary documentation for Cactus thorn exact.

CG 4.11.97, revised 28.1.98 (PW 28.1.98, 4.9.98 on blend BC PW 2.2.98 on Comp)

This thorn checks cactus against an exact solution.

1) An exact solution with a boost and a rotation Killing vector. The
solution comes out of Eqns. (2.1), (2.8) and (2.10) of J. Bicak,
C. Honselaers and B. G. Schmidt, Proc. Royal Soc. London A 390,
411-419 (1983).

2) Schwarzschild in Eddington-Finkelstein coordinates.

One can separately ask the thorn to provide initial data, the lapse
and shift, and boundary data, from the exact solution. This means that
cactus can be tested without the complications of, say, maximal
slicing, or implementing an outer boundary conditions. But one can
also test these, using the thorn only to obtain initial data.

Look at exact\_rfr.F. There are three decisions:

1) Set initial data (gij and Kij, and Vs and Ds of BM) and initial
lapse and shift from the exact solution? Yes, if model =
"exact". In this case, register exactinitialize with rfr as
CACTUS\_INITIAL.  Unless lapseinit = "exact" and shift = "exact"
are also set, there will be an error message, but other values will be
tolerated. IMPORTANT: The thorn then assumes that after it has been
called, something else will overwrite the lapse and shift.

2) Are lapse and shift to be set from the exact solution during
runtime? Yes, if slicing = "exact" and shift = "exact". If one
is set, the other must be, too. Register exactgauge with
rfr as	CACTUS\_PRESTEP. Other gauge choices, eg. minmax or geodesic,
are possible.

3) Are boundary data to be obtained from the exact solution. Yes, if
bound = "exact". Then register exactboundary with rfr as
CACTUS\_BOUND. Another choice could be "flat".

4) Blending boundary conditions are available. These will set the
values of g, K, v, and D to an exact solution outside some radius, the
evolved value inside some radius, and a linear blend in between. They
are controlled by various parameters starting with exblend. See
exactblendbound.F for more documentation and comment. To activate
blended boundary conditions use

bound = "exact blend"

This blending technique is the same as the GC; on a radius between r1 and   
r2 a linear combination of the numerical and exact solution is placed on  
the grid, scaled in r (so at r1 you get all numerical, at r2, all exact,  
and half way between 50/50). Above r2 the exact solution is simply 
injected onto the grid.        

You can also blend in a cartesian direction (eg, only boundary layers
within n * dx of the grid) using bound = "exact cartblend". Control n
with exblend\_width, with $< 0$ meaning width * dx, and $> 0$ meaning width
in coordinate space.

Generally, this thorn could be broken by other thorns overwriting what
it does when they should not, or by not overwriting the initial lapse
and shift if they should. Therefore CAREFUL with the order of calling.

Exact Cauchy data, lapse and shift, and BM variables are obtained from
a routine that provides the 4-metric and its inverse. Its output is
differentiated numerically where necessary. This poses no problem, as
we are numerically differentiating analytic expressions with a very
small epsilon (not on the grid!)

Additionally this routine will compute convergence coefficients and
differences between the exact solution in runtime, a functionality
previously in thorn\_Comparison.  To use this feature, you need to run
in convergence mode (convergence="yes") and then activate the
following parameters:

comparison = "yes"
comp\_fields = "gxx gyy" 

If you list the comp\_fields in the slice fields then differences
between the fields and the exact solutions will be output. The
comparison is done against the exactmodel.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Code structure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Code structure for faking a GR time evolution by evolving a surface
embedded in an exact solution. 

model = "slice" extracts Cauchy data from an arbitrary embedded slice 
in an exact solution. The solution in question is selected by, eg.
exactmodel = "fakebinary". model="exact" is now in effect a subclass of
model = "slice", although implemented independently.

system = "slice" then fakes a time evolution, by evolving the slice in
the exact solution, for arbitrary given lapse and shift, and
extracting Cauchy data at each time step. model = "slice" is then
implied and redundant.

The use of this is to test BM and ADM against exact solutions, but
mainly to test new elliptic gauge conditions in a 3D setting without
having to evolve BH spacetimes. Note that the "exact solution" need
only be a 4-metric given in a subroutine, but it need not obey the
Einstein equations. Examples so far are Miguel's bowl metric and Kip's
fake binary metric (just coded up by Bernd in thorn\_exact.)

Notation: $x^A = X^A(x^i,t)$. A=1,4 corresponding to X,Y,Z,T. 

Grid functions: $x^A$ are stored in GFs slicex, slicey, slicez,
slicet.  Temporary storage slicetmpx, slicetmpy, slicetmpz,
slicetmpt. These two groups are referred to in the following as SLICE
and TMP.

{\obeylines

CACTUS\_INITIAL slice\_initialize: 
\quad	initialize $x^A$ in SLICE
\quad	call slice\_data
\quad	call boundaries

CACTUS\_EVOL slice\_evolve: 
\quad	Loop over all points:
\quad\quad		evolve from SLICE to TMP1 using time derivative in TMP2
\quad\quad\quad			(Lax step, or forward in time to mid point)
\quad	end do	
\quad	call linextraponebound on SLICE
\quad	Loop over interior points:
\quad\quad		calculate $x^A_{,i}$ from TMP1	
\quad\quad		include slice\_normal
\quad	end do	
\quad	call linextraponebound on TMP2
\quad	Loop over all points:
\quad\quad		update SLICE using TMP2
\quad\quad\quad			(leapfrog, or centered in time, using $dx^A/dt$ in temp)
\quad	end do
\quad	call linextraponebound on SLICE
\quad	call slice\_data	

slice\_data:
\quad	Loop over interior points
\quad\quad		calculate $x^A_{,i}$ from SLICE
\quad\quad		include slice\_normal
\quad\quad		calculate $g_ij$
\quad\quad		calculate $x^A_{,ij}$ from SLICE
\quad\quad		calculate $g_{AB,C}$
\quad\quad		calculate $K_ij$
\quad\quad		(if BM) calculate $D_{ijk}$, $V_i$
\quad	end do	
\quad	call linextraponebound on TMP2
			
slice\_normal: to be used inside a loop over interior points.
\quad		calculate $n_A$
\quad		call exactmetric
\quad		calculate $n^A$
\quad		calculate $dx^A/dt$ in TMP2

}
						


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%