path: root/doc/documentation.tex

\documentclass{article}
\usepackage{graphicx}

\begin{document}

\title{Method of Lines}
\author{Ian Hawke}
\date{$ $Id$ $}
\maketitle

\abstract{The Method of Lines is a way of separating the time
  integration from the rest of an evolution scheme. This thorn is
  intended to take care of all the bookwork and provide some basic
  time integration methods, allowing for easy coupling of different
  thorns.} 

\section{Purpose}
\label{sec:purpose}

The Method of Lines (MoL) converts a (system of) partial differential
equation(s) into an ordinary differential equation containing some
spatial differential operator. As an example, consider writing the
hyperbolic system of PDE's
\begin{equation}
  \label{eq:mol1}
  \partial_t {\bf q} + {\bf A}^i({\bf q}) \partial_i {\bf B}({\bf q})
  = {\bf s}({\bf q})
\end{equation}
in the alternative form
\begin{equation}
  \label{eq:mol2}
  \partial_t {\bf q} = {\bf L}({\bf q}),
\end{equation}
which (assuming a given discretization of space) is an ODE.

Given this separation of the time and space discretizations, well
known stable ODE integrators such as Runge-Kutta can be used to do the
time integration. This is more modular (allowing for simple extensions
to higher order methods), more stable (as instabilities can now only
arise from the spatial discretization or the equations themselves) and
also avoids the problems of retaining high orders of convergence when
coupling different matter models.

MoL can be used for hyperbolic, parabolic and even elliptic problems
(although I definitely don't recommend the latter). As it currently
stands it is set up for systems of equations in the first order type
form of equation~(\ref{eq:mol2}). If you want to implement a
multilevel scheme such as leapfrog it is not obvious to me that MoL is
the thing to use. However if you have lots of thorns that you want to
interact, for example ADM\_BSSN and a hydro code plus maybe EM or a
scalar field, and they can easily be written in this sort of form,
then you probably want to use MoL.

This thorn is meant to provide a simple interface that will implement
the MoL inside Cactus as transparently as possible. It will initially
implement only the optimal Runge-Kutta time integrators (which are TVD
up to RK3, so suitable for hydro) up to fourth order and iterated
Crank Nicholson. All of the interaction with the MoL thorn should
occur directly through the scheduler. For example, all synchronization
steps should now be possible at the scheduler level.

For more information on the Method of Lines the most comprehensive
references are the works of Jonathan Thornburg~\cite{jt1,jt2} - see
especially section 7.3 of the thesis. From the CFD viewpoint the
review of ENO methods by Shu,~\cite{shu} has some information. For
relativistic fluids the paper of Neilsen and Choptuik~\cite{nc} is
also quite good. For Ian Hawke's viewpoint either see below or contact
by email.

\section{How to use}
\label{sec:use}


\subsection{Thorn users}
\label{sec:useruse}

For those who used the old version of MoL, this version is
unfortunately slightly more effort to use. That is, for most methods
you'll now have to set 4 parameters instead of just one. This
documentation is still in the process of conversion...

If you already have a thorn that uses the method of lines, then there
are four main parameters that are relevant to change the integration
method. The keyword {\tt MoL\_ODE\_Method} chooses between the
different methods. Currently supported are {\tt RK2}, {\tt ICN} and
{\tt Generic}. These are second order Runge-Kutta, Iterative Crank
Nicholson, and the generic Shu-Osher type Runge-Kutta methods. To
switch between the different types of generic methods there is also
the keyword {\tt Generic\_Type} which is currently restricted to {\tt
  RK} for the standard TVD Runge-Kutta methods (first to fourth order)
and {\tt ICN} for the implementation of the Iterative Crank Nicholson
method in generic form.

The parameter {\tt MoL\_Intermediate\_Steps} controls the number of
intermediate steps for the ODE solver. For the generic Runge-Kutta
solvers it controls the order of accuracy of the method.  For the {\tt
  ICN} methods this parameter controls the number of iterations taken,
which {\bf does not check for stability}. This parameter defaults to
3.

The parameter {\tt MoL\_Num\_Scratch\_Levels} controls the amount of
scratch space used. If this is insufficient for the method selected
there will be an error at parameter checking time. This parameter
defaults to 0, as no scratch space is required for the efficient ICN
and Runge-Kutta 2 solvers. For the generic solvers this must be at
least {\tt MoL\_Intermediate\_Steps - 1}.

A final parameter is {\tt MoL\_Memory\_Always\_On} which switches on
memory for the scratch space always if true and only during evolution
if false. This defaults to true for speed reasons; the memory gains
are likely to be limited unless you're doing something very memory
intensive at initialization or analysis.

\subsection{Thorn writers}
\label{sec:writeruse}

To port an existing thorn using the method of lines, or to write a new
thorn using it, should hopefully be relatively simple. As an example,
within the MoL arrangement is WaveMoL which duplicates the WaveToy
thorn given by CactusWave in a form suitable for use by MoL. In this
section, ``the thorn'' will mean the user thorn doing the physics.

We start with some terminology. Grid functions are split into four
cateogories.
\begin{enumerate}
\item The first are those that are evolved using a MoL form. That is,
  a right hand side is calculated and the variable updated using
  it. These we call {\it evolved} variables.
\item The second category are those variables that are set by a thorn
  at every intermediate step of the evolution, usually to respect the
  constraints. Examples of these include the primitive variables in a
  hydrodynamics code. Another example would be the gauge variables if
  these were set by constraints at every intermediate step (which is
  slightly artificial; the usual example would be the use of maximal
  slicing, which is only applied once every $N$ complete steps). These
  are known as {\it primitive} variables in analogy with the hydro
  code.
\item The third category are those variables that a thorn depends on
  but does not set or evolve. An example would include the metric
  terms considered from a thorn evolving matter. These are known as
  {\it dependent} variables.
\item The final category are those variables that do not interact with
  MoL. These would include temporary variables for analysis or setting
  up the initial data. These can safely be ignored.
\end{enumerate}

MoL needs to know every GF that falls in one of the first three
groups. If a GF is evolved by one thorn but is a dependent variable in
another (for example, the metric in full GR Hydro) then each thorn
should register the function as they see it. For example, the hydro
thorn will register the metric as a dependent variable and the
spacetime thorn will register the metric as an evolved variable. The
different variable categories are given the priority evolved,
primitive, dependent. So if a variable is registered as belonging in
two different categories, it is always considered by MoL to belong to
the category with the highest priority. 

Every thorn should register every grid function that it uses even if
you expect it to be registered again by a different thorn. For
example, a hydro thorn would register the metric variables as
dependent, whilst the spacetime evolution thorn would register them as
evolved (in ADM) or primitive (in ADM\_BSSN), both of which have
precedence. To register your GFs with MoL schedule a routine in the
bin {\tt MoL\_Register} which just contains the relevant function
calls.  For an evolved variable the GF corresponding to the update
term (${\bf L}({\bf q})$ in equation~(\ref{eq:mol2})) should be
registered at the same time, using the call {\tt MoL\_RegisterVar($q,
  L(q)$)}. For the primitive and dependent variables the function
calls are {\tt MoL\_RegisterPrimitive($q$)} and {\tt
  MoL\_RegisterDepends($q$)} respectively. Here $q$ is the {\it
  variable index} of the GF you wish to register, and $L(q)$ the
variable index of the associated right hand side GF. These are found
using the {\tt CCTK\_VarIndex*} functions.

There are both C and Fortran versions of these functions; see the
reference in section~\ref{sec:molfns}. For the C calls the prototypes
are given in the file {\tt MoL.h} which can be included without hard
paths by adding the line

{\tt USES INCLUDE: MoL.h}

\noindent to the {\tt interface.ccl} file of your thorn.

{\bf UPDATE} This has now been switched to using function
aliasing. The include file still exists, but is commented out in the
{\tt interface.ccl}. For an example of how to use the aliased
functions, see thorn WaveMoL.

Having done that, one routine (or group of routines) which we'll here
call {\tt Thorn\_CalcRHS} must be defined. This does all the finite
differencing that you'd usually do, applied to ${\bf q}$, and finds
the right hand sides which are stored in ${\bf L}$. This routine
should be scheduled in {\tt MoL\_CalcRHS}. The precise order that
these are scheduled should not matter, because no updating of any of
the user thorns ${\bf q}$ will be done until after all the RHSs are
calculated. {\bf Important note:} all the finite differencing must be
applied to the most recent time level ${\bf q}$ and not to the
previous time level ${\bf q}_p$ as you would normally do. Don't worry
about setting up the data before the calculation, as MoL will do that
automatically.

Finally, if you have some things that have to be done after each
update to an intermediate level, these should be scheduled in {\tt
  MoL\_PostStep}. Examples of things that need to go here include the
recalculation of primitive variables for hydro codes, the application
of boundary conditions\footnote{It is possible to alter the
  calculation of {\bf L} so that boundary conditions are automatically
  updated and do not need setting. This is slightly tricksy. For an
  example of how this would work see the new radiative boundary
  condition in ADM\_BSSN. For more on this see section 7.3.4
  of~\cite{jt1}.}, the solution of elliptic equations (although this
would be a very expensive place to solve them, some sets of equations
might require the updating of some variables by constraints in this
fashion). When applying boundary conditions the cleanest thing to do
is to write a routine applying the symmetries to the appropriate GFs
and, when calling it from the scheduler, adding the {\tt SYNC}
statement to the appropriate groups. An example is given by the
routine {\tt WaveToyMoL\_Boundaries} in thorn WaveMoL.

Points to note. The thorn routine {\tt Thorn\_CalcRHS} does not need
to know and in fact should definitely not know where precisely in the
MoL step it is. It just needs to know that it is receiving {\it some}
intermediate data stored in the GFs ${\bf q}$ and that it should
return the RHS ${\bf L}({\bf q})$. All the book-keeping to ensure that
it is passed the correct intermediate state at that the GFs contain
the correct data at the end of the MoL step will be dealt with by the
MoL thorn and the flesh. Also the synchronization of grids across
separate processors will be dealt with by the MoL thorn and the flesh.

Final point. Currently this MoL thorn is not guaranteed to work with
mesh refinement codes. Aside from a number of little details that need
to be worked out, currently MoL does not ensure that {\tt
  cctk\_delta\_time} is correctly set in the intermediate stages of
the calculation. This will hopefully change, but for the immediate
future don't expect this to work with any mesh refinement code.

Update. It is now possible to change the category of a GF in the
middle of a run. This was added solely to allow for mixed slicing for
Einstein. In this case the lapse is an evolved function if using
$1+\log$ slicing, a primitive variable if using maximal or static
slicing, and a dependent variable if evolved in some unknown
manner. Currently there exist functions for changing to any of the
four categories listed above.

{\bf I strongly suggest that these functions are not used}. They are
almost completely untested, and even those that are tested are far
from stable. Also, when moving to the version of this thorn that
supports Mesh Refinement it is quite possible that these functions may
be dropped.


\section{Example}
\label{sec:example}

As a fairly extended example of how to use MoL I'll outline how
ADM\_BSSN works in this context. The actual implementation of this is
given in the thorn {\tt AEIDevelopment/BSSN\_MoL}. 

As normal the required variables are defined in the {\tt
  interface.ccl} file, together with the associated source terms. For
example, the conformal factor and source are defined by

\begin{verbatim}
real ADM_BSSN_phi type=GF timelevels=2
{
  ADM_BS_phi
} "ADM_BSSN_phi"

real ADM_BSSN_sources type=GF
{
...,
  adm_bs_sphi,
...
}
\end{verbatim}
Also in this file we write {\tt USES INCLUDE: MoL.h}.

Once the sources are defined the registration with MoL is required,
for which the essential file is {\tt MoLRegister.c}. As we're
registering from C we ensure that we {\tt \#include "MoL.h"}. In the
ADM\_BSSN system the standard metric coefficients $g_{ij}$ are not
evolved, and neither are the standard extrinsic curvature components
$K_{ij}$.  However these are used by ADM\_BSSN in a number of places,
and are calculated from evolved quantities at the appropriate points.
In the MoL terminology these variables are {\it primitive}. As the
appropriate storage is defined in thorn Einstein, the actual calls
have the form

\begin{verbatim}
 ierr += MoL_RegisterPrimitive(CCTK_VarIndex("einstein::kxx"));
\end{verbatim}

\noindent The actual evolved variables include things such as the
conformal factor. This, and the appropriate source term, is defined in
thorn ADM\_BSSN, and so the call has the form

\begin{verbatim} 
 ierr += MoL_RegisterVar(CCTK_VarIndex("adm_bssn::ADM_BS_phi"),
                         CCTK_VarIndex("adm_bssn::adm_bs_sphi")); 
\end{verbatim}


As well as the evolved variables, and those constrained variables such
as the metric, there are the gauge variables. Precisely what status
these have depends on how they are set. If harmonic or 1+log slicing
is used then the lapse is evolved:

\begin{verbatim}
 ierr += MoL_RegisterVar(CCTK_VarIndex("einstein::alp"),
                         CCTK_VarIndex("adm_bssn::adm_bs_salp")); 
\end{verbatim}

\noindent If maximal or static slicing is used then the lapse is a
primitive variable\footnote{Note that this is actually a bit of a
  hack. The rational for {\it dependent} variables was to deal with
  maximal slicing. However it turned out that I hadn't thought it
  through correctly and that the treatment for primitive variables was
  required.}:

\begin{verbatim}
 ierr += MoL_RegisterPrimitive(CCTK_VarIndex("einstein::alp"));
\end{verbatim}

\noindent Finally, if none of the above apply we assume that the lapse
is evolved in some unknown fashion, and so it must be registered as a
dependent variable:

\begin{verbatim}
 ierr += MoL_RegisterDepends(CCTK_VarIndex("einstein::alp"));
\end{verbatim}

However, it is perfectly possible that we may wish to change how we
deal with the gauge during the evolution. This is dealt with in the
file {\tt PreLoop.F}. If the slicing changes then the appropriate
routine is called. For example, if we want to use 1+log evolution then
we call 

\begin{verbatim}
 call CCTK_VarIndex(lapseindex,"einstein::alp")
 call CCTK_VarIndex(lapserhsindex,"adm_bssn::adm_bs_salp")
 call MoL_ChangeVarToEvolved(ierr, lapseindex, lapserhsindex)
\end{verbatim}

\noindent It is not required to tell MoL what the lapse is changing
{\it from}, or indeed if it is changing at all; MoL will work this out
for itself.

Finally there are the routines that we wish to apply after every
intermediate step. These are {\tt ADM\_BSSN\_removetrA} which enforces
various constraints (such as the tracefree conformal extrinsic
curvature remaining trace free), {\tt ADM\_BSSN\_Boundaries} which
applies symmetry boundary conditions as well as various others (such
as some of the radiative boundary conditions). Note all the calls to
{\tt SYNC} at this point. We also convert from the updated BSSN
variables back to the standard ADM variables in {\tt
  ADM\_BSSN\_StandardVariables}, and also update the time derivative
of the lapse in {\tt ADM\_BSSN\_LapseChange}.

\section{The gory details}
\label{sec:gory}

This section is more for people wishing to maintain and extend the
code. {\bf NEEDS TOTAL REWRITE FOR MOL2}

There are three essential parts to the MoL thorn. The first is the
registration process where the MoL thorn is told just how many GFs it
should be looking to integrate, and where they and their RHSs
are. This is split into three separate routines.
\begin{itemize}
\item {\tt MoL\_Init\_Register} Called by the MoL thorn at {\tt
    postinitial}.  Initializes the number of variables to the maximum
  possible number (given by {\tt CCTK\_NumVars()}). Sets up two arrays
  of pointers to GFs (CCTK\_Real**), one for the data ({\tt qindex})
  and one for the RHSs ({\tt qrhsindex}). It also sets up index arrays
  for the primitive and dependent variables.
\item {\tt MoL\_Register} This group lets MoL knows where all the data
  is and should contain all the user thorn calls to {\tt
    MoL\_RegisterVar()}. This just associates the index of the
  GFs to be evolved with their RHSs, and increments the number of
  variables to be evolved. Also contains the routines to register the
  primitive and dependent variables.
\item {\tt MoL\_End\_Register} This could be used to reallocate the
  index arrays to the correct size. For the moment I've decided that
  the tiny memory gain is outweighed by the possiblity of problems
  using realloc.
\end{itemize}

The second part of the MoL thorn sets up all the scratch space and the
parameters for the generalized RK routines. This is also scheduled at
postinitial. This contains the routines

\begin{itemize}
\item {\tt MoL\_SetupScratch} Sets up the scratch space. Sets the
  scalar {\tt mol\_num\_scratch} which is the size of the required
  scratch space from the defaults, then allocates the scratch arrays. 
\item {\tt MoL\_RKSetup} Sets up the $\alpha$ and $\beta$ arrays for
  the generalized Runge Kutta step. Note that as we're currently only
  storing the most recent RHS arrays the $\beta$ arrays are single
  index only.
\end{itemize}

The final part of the MoL thorn is the evolution step. This is again
split into a number of separate routines and groups. Firstly we
describe a single evolution step, and what MoL expects.

By a single evolution step we here mean the execution of everything
inside the schedule bin {\tt EVOL}. MoL expects that every GF that it
knows about is allocated, initialized and lives at the same instant of
computational time. MoL also expects the ghost zones and boundaries to
be set correctly. Wherever possible the first two are checked. MoL
expects that the driver has rotated the timelevels so that the last
set of complete, consistent data is stored in ${\bf q}\_p$. MoL also
knows that the dependent variables may or may not have been evolved by
the time the MoL evolution group is executed.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=3cm]{MoLdia1}
\includegraphics[width=3cm]{MoLdia2}
\includegraphics[width=3cm]{MoLdia3}
\end{center}
\caption{How MoL treats the three different types of variables. The
  MoL step is performed at evolution after the driver has rotated the
  timelevels, which occurs right at the start of every {\tt
    CCTK\_EVOL} step. The physics thorns expect the most recent data
  to be at the current time level (the top solid line). So the first
  step for most types of variable is to copy or pointer switch the
  data from the previous time level (bottom solid line) where the most
  recent data exists after rotation. Figure (a) shows the {\it
    evolved} variables.  At each intermediate step the data is updated
  into the current time level and, if necessary, stored in scratch
  space.  Figure (b) shows the {\it primitive} variables. As these are
  set by the physics thorns into the current time level all that is
  required is the initial copy. Figure (c) shows the {\it dependent}
  variables for which we cannot know whether these variables are
  evolved before or after MoL at evolution, and hence whether the
  current data is already filled before MoL starts, we must store the
  current data in scratch space first, then do the copy from the
  previous level, and then at the end of the MoL step return the data
  to the initial state.}
\label{MoLvariables}
\end{figure}

MoL plays around with the timelevels during the intermediate
steps. This is required so that the latest level in the MoL evolution
is always stored in the same place so that the user thorns can easily
access them. This place is the current time level ${\bf q}$. This
ensures that MoL works correctly with all the standard boundary
condition routines and minimizes the effort of porting non timelevel
aware thorns (I hope!).

An outline of the schedule is as follows:
\begin{equation}
  \label{eq:schedoutline}
  \begin{array}[l]{l}
    \texttt{MoL\_StartStep} \\
    \texttt{MoL\_Step WHILE counter \{} \\
    \begin{array}[l]{l}
      \texttt{MoL\_CalcRHS \{ \}} \\
      \texttt{MoL\_Add} \\
      \texttt{MoL\_PostStep \{ \}} \\
    \end{array} \\
    \texttt{\}} \\
    \texttt{MoL\_EndStep \{ \}}
  \end{array}
\end{equation}

Each different type of variable is treated slightly differently by
MoL. Each is assumed to have at least 2 timelevels (although this is
checked, a fatal error occurs otherwise). Before entering the loop
over the intermediate steps MoL will first copy (pointer switch?) the
data into the current time levels. Before doing this the current
timelevel of any {\it dependent} variables is copied to scratch space,
as it may have already been updated. During the loop over the
intermediate steps only the {\it evolved} variables are directly
altered by MoL. When all user thorns have given their right hand side
GFs the evolved variables are updated into the current time level.
This data may also be copied to scratch space if required for later.
The primitive variables are assumed to be set by the user thorn,
either in the calculation of the RHS or during {\tt MoL\_PostStep} at
the end of each intermediate step. The dependent variables are assumed
to remain completely unchanged. After all the intermediate steps the
data in the current and previous timelevels is ``correct'' and
consistent for the evolved and the temporary variables. The data for
the current timelevel for the dependent variables is recovered from
scratch space. Precisely how the different variables are
treated is shown in figure~\ref{MoLvariables}.

\begin{itemize}
\item {\tt MoL\_StartStep} This ensures that the integer keeping
    track of where we are in the MoL step is set to the correct
    value. It also copies the previous data ${\bf q}_p$ into the
    current position ${\bf q}$ ready for the first step.
\item {\tt MoL\_Step} This is the main part of the thorn. The scheduler
  allows us to loop over this group the correct number of times to
  complete a single Cactus evolution step. Contained within this group
  are:
  \begin{itemize}
  \item {\tt MoL\_CalcRHS} The schedule group within which the user
    thorns will schedule their routines. These routines should
    calculate the GFs ${\bf L}({\bf q})$. The order these routines are
    scheduled within this group should not matter.
  \item {\tt MoL\_Add} This step performs the time integration
    depending on where in the MoL step we are. Updates directly into
    the current GF and copies to the scratch space if required.
  \item {\tt MoL\_PostStep} Another schedule group within which the
    user thorns can schedule such things as primitive variable
    recovery, boundary conditions, etc.
  \item {\tt MoL\_End} Just alters the scalar tracking the position in
    the MoL loop.
  \end{itemize}
\end{itemize}

Finally, there are the routines {\tt MoL\_FreeScratch} and {\tt
  MoL\_RKFree} which free up the memory that was explicitly taken. For
the moment these routines are scheduled at postevol.


\section{Adding new numerical integrators}
\label{sec:newmeth}

There are two obvious ways of adding new ODE integrators into MoL. The
first is to follow the route used in the efficient RK2 and ICN
methods. That is, you let the underlying infrastructure define the
scratch space but you do all the addition yourself. It's probably best
if you model your integrator on one of the efficient routines to start
with. 

The alternative is to use the generic integrator. To use this you just
need to add the correct set of $\alpha$ and $\beta$ coefficients so
that the generic routine can perform the additions. You'll also have
to set up the keyword parameters and so on.

\section{To do list}
\label{sec:todo}

\begin{itemize}
\item The documentation must be improved, especially inside the code
  itself. 
\item Errors are currently not handled well (if at all). This must be
  fixed.
\item A test suite is required. I'm not sure how to do this without
  using WaveMoL, but we could try.
\item In order to make the code work with a Mesh Refinement driver,
  the scratch spaces must be changed to be grid functions. This is
  currently impossible, but Tom Goodale will add the appropriate bits
  to the flesh to make it possible to do. At that point the entire
  code will probably be rewritten.
\end{itemize}

\section{Functions provided by MoL}
\label{sec:molfns}

{\bf Note for Cactus people: I'd really like there to be a generic
  style file so that I could use the same function environments as in
  the User's Guide. It would also be nice if the parsing of the
  interface files recognized function aliasing when that appears, but
  that's a long distant wish}.

All the functions listed below return error codes in theory. However
at this current point in time they always return 0 (success). Any
failure to register or change a GF is assumed fatal and MoL will
issue a level 0 warning stopping the code. This may change in future,
in which case negative return values will indicate errors.

\subsection*{MoL\_RegisterVar}

Tells MoL that the given GF is in the evolved category with the
associated update GF. 

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_RegisterVar(CCTK\_INT varindex,
  CCTK\_INT rhsindex) \\
  & {\bf Fortran} & MoL\_RegisterVar(CCTK\_INT ierr, CCTK\_INT varindex,
  CCTK\_INT rhsindex) 
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the GF to be evolved \\
  & {\bf rhsindex} & - & Index of the associated update GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_Register}.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_RegisterVar(CCTK\_VarIndex("wavetoymol::phi"), \\
  && \qquad CCTK\_VarIndex("wavetoymol::phirhs")); \\ 
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``wavetoymol::phi'')
  \\
  && call CCTK\_VarIndex(rhsindex, ``wavetoymol::phirhs'') \\
  && call MoL\_RegisterVar(ierr, varindex, rhsindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterDepends, MoL\_RegisterPrimitive, \\
  & MoL\_ChangeVarToEvolved.
\end{tabular}


\subsection*{MoL\_RegisterDepends}

Tells MoL that the given GF is in the dependent category.

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_RegisterDepends(CCTK\_INT varindex) \\
  & {\bf Fortran} & MoL\_RegisterDepends(CCTK\_INT ierr, CCTK\_INT varindex) 
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the dependent GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_Register}.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_RegisterDepends(CCTK\_VarIndex("einstein::alp")); \\ 
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'')
  \\
  && call MoL\_RegisterDepends(ierr, varindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterVar, MoL\_RegisterPrimitive, \\
  & MoL\_ChangeVarToDependent.
\end{tabular}


\subsection*{MoL\_RegisterPrimitive}

Tells MoL that the given GF is in the primitive category.

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_RegisterPrimitive(CCTK\_INT varindex) \\
  & {\bf Fortran} & MoL\_RegisterPrimitive(CCTK\_INT ierr, CCTK\_INT varindex) 
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the primitive GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_Register}.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_RegisterPrimitive(CCTK\_VarIndex("einstein::alp")); \\ 
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'')
  \\
  && call MoL\_RegisterPrimitive(ierr, varindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterVar, MoL\_RegisterDepends, \\
  & MoL\_ChangeVarToPrimitive.
\end{tabular}


\subsection*{MoL\_ChangeVarToEvolved}

Sets a GF to belong to the evolved category, with the associated
update GF. Not used for the initial setting.

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_ChangeVarToEvolved(CCTK\_INT varindex, CCTK\_INT
  rhsindex) \\ 
  & {\bf Fortran} & MoL\_ChangeVarToEvolved(CCTK\_INT ierr, CCTK\_INT varindex,
  CCTK\_INT rhsindex)  
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the GF to be evolved. \\
  & {\bf rhsindex} & - & Index of the associated update GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_PreStep}. \\
  & Note that this function was designed to allow mixed slicings for
  thorn Einstein. \\
  & This set of functions is largely untested and should be used with
  great care.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_ChangeVarToEvolved(CCTK\_VarIndex("einstein::alp"),\\
  && \qquad CCTK\_VarIndex(``adm\_bssn::adm\_bs\_salp'')); \\
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'') \\
  && call CCTK\_VarIndex(rhsindex, ``adm\_bssn::adm\_bs\_salp'') \\
  && call MoL\_ChangeVarToEvolved(ierr, varindex, rhsindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterVar, MoL\_ChangeVarToDependent, \\
  & MoL\_ChangeVarToPrimitive, MoL\_ChangeVarToNone.
\end{tabular}


\subsection*{MoL\_ChangeVarToDependent}

Sets a GF to belong to the dependent category. Not used for the
initial setting. 

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_ChangeVarToDependent(CCTK\_INT varindex) \\ 
  & {\bf Fortran} & MoL\_ChangeVarToDependent(CCTK\_INT ierr, CCTK\_INT varindex)  
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the dependent GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_PreStep}. \\
  & Note that this function was designed to allow mixed slicings for
  thorn Einstein. \\
  & This set of functions is largely untested and should be used with
  great care.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_ChangeVarToDependent(CCTK\_VarIndex("einstein::alp")); \\
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'') \\
  && call MoL\_ChangeVarToDependent(ierr, varindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterDepends, MoL\_ChangeVarToEvolved, \\
  & MoL\_ChangeVarToPrimitive, MoL\_ChangeVarToNone.
\end{tabular}

\subsection*{MoL\_ChangeVarToPrimitive}

Sets a GF to belong to the primitive category. Not used for the
initial setting. 

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_ChangeVarToPrimitive(CCTK\_INT varindex) \\ 
  & {\bf Fortran} & MoL\_ChangeVarToPrimitive(CCTK\_INT ierr, CCTK\_INT varindex)  
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the primitive GF.
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_PreStep}. \\
  & Note that this function was designed to allow mixed slicings for
  thorn Einstein. \\
  & This set of functions is largely untested and should be used with
  great care.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_ChangeVarToPrimitive(CCTK\_VarIndex("einstein::alp")); \\
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'') \\
  && call MoL\_ChangeVarToPrimitive(ierr, varindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_RegisterPrimitive, MoL\_ChangeVarToEvolved, \\
  & MoL\_ChangeVarToDependent, MoL\_ChangeVarToNone.
\end{tabular}

\subsection*{MoL\_ChangeVarToNone}

Sets a GF to belong to the ``unknown'' category. Not used for the
initial setting. 

\vskip 3mm

\begin{tabular}[l]{l l l}
  {\bf Synopsis} && \\
  & {\bf C} & CCTK\_INT ierr = MoL\_ChangeVarToNone(CCTK\_INT varindex) \\ 
  & {\bf Fortran} & MoL\_ChangeVarToNone(CCTK\_INT ierr, CCTK\_INT varindex)  
\end{tabular}

\begin{tabular}[l]{l l c l}
  {\bf Arguments} &&& \\
  & {\bf varindex} & - & Index of the GF to be unset. 
\end{tabular}

\begin{tabular}[l]{l l}
  {\bf Discussion} & \\
  & Should be called in a function scheduled in {\tt MoL\_PreStep}. \\
  & Note that this function was designed to allow mixed slicings for
  thorn Einstein. \\
  & This set of functions is largely untested and should be used with
  great care.
\end{tabular}

\begin{tabular}[l]{l l l}
  {\bf Examples} && \\
  & {\bf C} & ierr =
  MoL\_ChangeVarToNone(CCTK\_VarIndex("einstein::alp")); \\
  & {\bf Fortran} & call CCTK\_VarIndex(varindex, ``einstein::alp'') \\
  && call MoL\_RegisterNone(ierr, varindex)
\end{tabular}

\begin{tabular}[l]{l l }
  {\bf See Also} & \\
  & CCTK\_VarIndex, MoL\_ChangeVarToEvolved, MoL\_ChangeVarToDependent, \\
  & MoL\_ChangeVarToPrimitive.
\end{tabular}




\begin{thebibliography}{1}

\bibitem{jt1}
J. Thornburg.
\newblock {N}umerical {R}elativity in {B}lack {H}ole {S}pacetimes. 
\newblock Unpublished thesis, University of British Columbia.
\newblock 1993.
\newblock Available from \mbox{\tt
  http://www.aei.mpg.de/\~{}jthorn/phd/html/phd.html}. 

\bibitem{jt2}
J. Thornburg.
\newblock A {3+1} {C}omputational {S}cheme for {D}ynamic {S}pherically
{S}ymmetric {B}lack {H}ole {S}pacetimes -- {II}: {T}ime {E}volution.
\newblock Preprint {\tt gr-qc/9906022}, submitted to {\em Phys. Rev.}
{\bf D}. 

\bibitem{shu}
C. Shu.
\newblock {H}igh {O}rder {ENO} and {WENO} {S}chemes for
{C}omputational {F}luid {D}ynamics.
\newblock In T.~J. Barth and H. Deconinck, editors {\em High-Order
  Methods for Computational Physics}. Springer, 1999.
\newblock A related online version can be found under {\em Essentially
  {N}on-{O}scillatory and {W}eighted {E}ssentially {N}on-{O}scillatory
  {S}chemes for {H}yperbolic {C}onservation {L}aws} at {\tt
  http://www.icase.edu/library/reports/rdp/97/97-65RDP.tex.refer.html}. 

\bibitem{nc}
D.~W. Neilsen and M.~W. Choptuik.
\newblock Ultrarelativistic fluid dynamics.
\newblock {\em Class. Quantum Grav.}, {\bf 17}:\penalty0 733--759, 2000.

\end{thebibliography}

\include{interface}
\include{param}
\include{schedule}

\end{document}