path: root/doc/ThornGuide.tex

\documentstyle{report}
\newcommand{\parameter}[1]{{\it #1}}

\begin{document}

\chapter{IDAnalyticBH}

\begin{tabular}{@{}ll}
Code Authors & Gerd Lanfermann \\
Maintained by & Cactus Developers \\
Documentation Authors & 
\end{tabular}

\section{Introduction}
  
\subsection{Purpose of Thorn}

Thorn EllBase provides the basic functionality to register a class of
ellitpic equations and register solver for that particular class.
The solvers are called by the user through a unique interface, which calls the 
elliptic solver requested by passing through the name under which the solver 
routine is  registered.
EllBase defines four basic elliptic classes and offers two solver 
(SOR and PETSc) for each class.

The elliptic classes supported by this thorn in three dimensions are
\begin{enumerate}
\item{\em flat:} solves a linear elliptic equation is flat space.
\item{\em metric:} solves a linear elliptic equation for a given metric
\item{\em conformal metric:} solves a linear elliptic equation for a
given metric and a conformal factor.
\item{\em generic:} solves a linear elliptic equation by passing the 
	stencil functions. There is support for a maximum of 27 stencil 
	functions ($3^3$).
\end{enumerate}

\subsection{Technical Specification}

\begin{itemize}

\item{Implements} ellbase
\item{Inherits from} grid
\item{Tested with thorns} EllTest

\end{itemize}

\subsubsection{Using EllBase}

EllBase by itself does not provide any ellitpic solving capabilities. 
It merely provides the registration structure and calling interface.

The idea of a unified calling interface can be motivated as follows: 
assume you a have elliptic problem which conforms to one of the elliptic 
classes defined in EllBase.

In your routine you provide the necessay information needed by a solver
for that particular elliptic class: you provide the function to solve for, 
the coefficient matrix, etc. You can now chose a solver, which is registered 
for that elliptic class. Later you  might want to compile with a different 
solver for this elliptic class: just change the name of the solver in your
elliptic interface call. If somebody improves a solver you have been using,
there is no need for you to change any code on your side: the interface 
will hide all of that. Another advantage is that your code will compile 
and run, even though certain solvers are not compiled in. In these case, you 
will have to do some return value checking to offer alterntives.

\subsubsection{Registration Mechanism}

Before a user can successfully apply a elliptic solver to one of his problems,
two things need to be done by the author who programms the solver.
\begin{itemize}
\item{\bf Register a class of elliptic equations} Depending on the elliptic 
problem This provides the unique calling, the solving routines needs to have 
specific input data. The interface, which is called by the user, has to 
reflect these arguments. EllBase already offers 
serveral of these interfaces, but if you need to have a new one, you can 
provide your own.

\item{\bf Register a solver for a particular elliptic equation classe}
Once a class of ellitpic equations has been made available as described 
above, the author can now register solvers for that particular class. 
Later a user will access the solver calling the interface with the 
arguments needed for the elliptic class and a name, under which a solver 
for this elliptic problem has been registered.
\end{itemize}

The registration process is part of the authors thorn, not part of EllBase. 
There is no need to change code in EllBase. Usually, a author of solver 
routines will register the routines that register an elliptic equation class 
and/or an elliptic solver at STARTUP (see \ref{se:RFR} for information on
 scheduling). If a author registers both, class and solver, you must make 
sure, that the elliptic class is registered {\em before} the solver 
registration takes place. 

\subsubsection{EllBase Programming Guide}

Here we give a step by step guide on how to implement an new elliptiv 
solver class, its interface and provide a solver for this class.  Since 
some of the functionality needed in the registration code can only be 
achieved in C, a basic knowledge of C is helpful.

\begin{itemize}
\item{\bf Assumption}: The elliptic equation we want to solve is of type XXX. 
The solver for this elliptic equation will be called ``FastSOR_solver''
and will be written in Fortran. Since Fortran cannot be called directly, we
need to have C wrapper function ``FastSOR_wrapper''. The elliptic equation 
class will be called ``SimpleEllClass'': it will be flat space solver, 
that only 
takes the coefficient matrix $M$. (Note that this solver class is already 
provided by
EllBase.)
The name of the demonstration thorn will be 
``ThornFastSOR''. Since I will only demonstrate the registration priciple 
and calling structure, I leave it to the interested reader to write a really 
fast SOR solver.

\item{\bf Elliptic class declaration}: {\tt SimpleEllThorn/src/SimpleEll_Class.c} 
\begin{verbatim}
\end{verbatim}

\item{\bf Elliptic solver interface}: {\tt src/SimpleEll_Interface.c}
\begin{verbatim}

#include ``cctk.h''
#include ``cctk_parameters.h''

#include ``cctk_FortranString.h'' 
#include ``StoreNamedData.h''

static pNamedData *SimpleEllSolverDB;

void Ell_SimpleEllSolverRegistry(void (*solver_func), const char *solver_name)
{
  StoreNamedData(&SimpleEllSolverDB,solver_name,(void*)solver_func);
}
\end{verbatim}
The routine above registers the solver (or better the function pointer of the solver routine ``*solve_func'') for the equation class 
{\em SimpleEllClass} by the name {\tt solver_name} in the database {\tt
SimpleEllSolverDB}. This database is declared in statment {\tt static pNamedData...}.


Next, we write our interface in the same file {\tt ./SimpleEll_Interface.c}:
\begin{verbatim}
void Ell_SimpleEllSolver(cGH *GH, int *FieldIndex, int *MIndex, 
		         CCTK_REAL *AbsTol, CCTK_REAL *RelTol,
                         const char *solver_name) {

/* prototype for the equation class wrapper:
   grid hierarchy(*GH), ID-number of field to solve for (*FieldIndex),
   two arrays of size three holding convergence information (*AbsTol, *RelTol)
*/
  void (*fn)(cGH *GH, int *FieldIndex, int *AbsTol, int *RelTol);

  /* derive the function name from the requested name and hope it is there */
  fn = (void(*)) GetNamedData(LinConfMetricSolverDB,solver_name);
  if (!fn) CCTK_WARN(0,''Ell_SimpleEllSolver: Cannot find solver! ``);

  /* Now that we have the function pointer to our solver, call the 
     solver and pass through all the necessary arguments */
  fn( GH, FieldIndex, MIndex, AbsTol, RelTol);
}

\end{verbatim}
The interface {\tt Ell_SimpleEllSolver} is called from the user side. It receives a pointer to the grid hierarchy, the ID-number of the field to solver for, two arrays which the used upload with convergence test info, and finally, the name of the solver the user want to employ {\tt *solver_name}. {\bf Note:} all these quantities are referenced by pointers, hence the ``*''.

Within the interface, the solver_name is used to get the pointer to function which was registered under this name.
Once the function is known, it called with all the arguments passed to 
the interface.

To allow calls from Fortran, the interface in C needs to be ``wrapped''. 
(This wrapping  is different from the one necessary to make to actual solver 
accessible by the elliptic registry).

\begin{verbatm}
/* Fortran wrappr for the routine Ell_SimpleEllSolver */
void FMODIFIER FORTRAN_NAME(Ell_SimpelEllSolver)
     (cGH *GH, int *FieldIndex, int *MIndex, 
     int *AbsTol, int *RelTol, ONE_FORTSTRING_ARG) {
  ONE_FORTSTRING_CREATE(solver_name);

  /* Call the interface */
  Ell_SimpleEllSolver(GH, FieldIndex, MIndex, AbsTol, RelTol, solver_name);
  free(solver_name);
}
\end{verbatim}

\item{\bf Elliptic solver}:{\tt ./src/FastSOR_solver.F}\\
Here we show the first lines of the Fortran code for the solver:

\begin{verbatim}
      subroutine FastSOR_solver(_CCTK_FARGUMENTS,
     .   Mlinear_lsh,Mlinear, 
     .   var,
     .   abstol,reltol)

      implicit none

      _DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      INTEGER CCTK_Equals

      INTEGER Mlinear_lsh(3)
      CCTK_REAL Mlinear(Mlinear_lsh(1),Mlinear_lsh(2),Mlinear_lsh(3))
      CCTK_REAL var(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3))
	
      INTEGER Mlinear_storage

c     We have no storage for M if they are of size one in each direction
      if ((Mlinear_lsh(1).eq.1) .and. 
     .   (Mlinear_lsh(2).eq.1)  .and.
     .   (Mlinear_lsh(3).eq.1)) then
         Mlinear_storage=0
      else
         Mlinear_storage=1
      endif


\end{verbatim}
This Fortran solver receives the following arguments: the ``typical''
 CCTK_ARGUMENTS({\tt  (see \ref{sec:cctk_arguments}): {\tt \_CCTK\_FARGUMENTS},
the {\em size} of the coefficient matrix: {\tt Mlinear_lsh}, the coefficient
matrix {\tt Mlinear}, the variable to solve for: {\tt var}, and the two arrays with convergence information. 

In the declaration section, we delcare: the cctk arguments, the Mlinear size array, the coefficient matrix, by the 3-dim. size array, the variable to solve for. Why do we pass the size of Mlinear explicitly and do not use the 
{\tt cctk_lsh} (processor local shape of a grid function) as we did for {\tt var} ? The reason is the following: while we can expect the storage of {\tt var} to be {\em on} for the solve, there is no reason (in a more general elliptic case) to assume, that the coefficient 
matrix has storage allocated, perhaps it is not needed at all! In this case, we have to protect ourself against referencing empty arrays. For this reason, we also employ the flag {\tt Mlinear_storage}.

\item{\bf Elliptic solver wrapper}:{\tt ./src/FastSOR_wrapper.c}\\
The Fortran solver can not be used within the elliptic registry directly. 
Instead the Fortran code is called through a wrapper:
\begin{verbatim}

void FastSOR_wrapper(cGH *GH, int *FieldIndex, int *MIndex, 
		     int *AbsTol,int *RelTol) {

  CCTK_REAL *Mlinear=NULL, *var=NULL;
  int Mlinear_lsh[3];
  int i;

  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex);

  if (*MIndex>0) Mlinear   = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*MIndex);


  if (GH->cctk_dim>3)
    CCTK_WARN(0,''This elliptic solver implementation does not do dimension>3!'');  
  
  for (i=0;i<GH->cctk_dim;i++) {
    if((*MIndex<0))  Mlinear_lsh[i]=1;
    else             Mlinear_lsh[i]=GH->cctk_lsh[i];
  }

  /* call the fortran routine */
  FORTRAN_NAME(SimpleEll_Solver)(_PASS_CCTK_C2F(GH),
         Mlinear_lsh, Mlinear, var, 	
	 AbsTol, RelTol);
}
\end{verbatim}

The wrapper {\tt FastSOR_wrapper} takes these arguments: the indices 
of the field to solve for ({\tt FieldIndex}) and the coefficient matrix 
({\tt MIndex}), the two arrays containing
convergence information ({\tt AbsTol, RelTol}). 
In the body of the programm we provide two CCTK_REAL pointers to the 
data section of the field to solver ({\tt var, Mlinear}) by means of {\tt Get_VarDataPtrI}. For Mlinear, we only do this, if the index is non-negative. 
A negative index is a signal by the user that the coefficient matrix has 
no storage allocated.(For more general elliptic equation cases, e.g. no 
source terms.)
 To make this information of a possibly empty matrix available to Fortran, we load a 3-dim. and pass this array thorugh to Fortran. See discussion above.


\item{\bf Elliptic solver startup}: {\tt ./src/Startup.c}\\
The routine below in {\tt Startup.c} performs the registration of our solver wrapper {\tt FastSOR_wrapper} under the name ``{\em fastsor}'' for the elliptic class ``{\em Ell_SimpleEll}''. We do not register with the solver interface 
{\tt Ell_SimpleEllSolver} directly, but with the class.
\begin{verbatim}
#include ``cctk.h''
#include ``cctk_parameters.h''

void FastSOR_register(cGH *GH) {

  /* protoype of the solver wrapper: */
  void FastSOR_wrapper(cGH *GH, int *FieldIndex, int *MIndex,
                      	int *AbsTol, int*RelTol);

  Ell_RegisterSolver(FastSOR_wrapper,''fastsor'',''Ell_SimpleEll'');
}
\end{verbatim}
Note that more solver registration code could be put here (registration 
for other classes, e.g.)


\item{\bf Elliptic solver scheduling}: {\tt schedule.ccl}
We schedule the registratio of the fast SOR solver at FIXME, by this time, 
{\em the elliptic class} {\tt Ell_SimpleEll} has already been registered. 
\begin{verbatim}
schedule FastSOR_register at CCTK_INITIAL
{
  LANG:C
} ``Register the fast sor solver''
\end{verbatim}



\end{document}