\chapter{Users Manual}\label{app:demo}
Here is a users guide to use WaVeMaKeR version $\epsilon$.
This version of WaVeMaKeR contains routines for computing the velocity
and velocity potential on a three dimensional periodic surface
wave. We also include copies of various driving routines which we have
used as well as some additional programs which may be of interest.
\section{Background}
The  routines are used to approximate integrals of the form 
\begin{equation}
I({\bf x})=\int Q({\bf x}^{\prime})\cdot \nabla_{\bf x}G({\bf
x-x^{\prime}})d\alpha^{\prime}
\label{defeq} 
\end{equation}
where
\begin{equation}
G({\bf x- x^{\prime}})=-{1 \over |{\bf x - x^{\prime}}|}.
\end{equation}
Here ${\bf x}=(x,y,z)$ and $Q({\bf x})=(q_1,q_2,q_3)$ are functions of
$\alpha=(\alpha_1,\alpha_2)$, the surface parameterization.
The equations defining the velocity potential and the velocity of
surface waves can be
defined in terms of integrals which take the form of $I$. 

The program was developed from the Parallel Multipole Treecode
Algorithm (PMTA). PMTA was developed at the Electrical Engineering
Department of Duke University. As a derivative of PMTA, we include the
documentation PMTAdoc.tex. 

There are a number of assumptions which we have used in developing the
code. We assume that the integrand of $I$ is periodic in the $x$ and $y$
dimensions. The code works by discretizing the integrand in a cube of
side CubeLength and then duplicating the box periodically in the $x$
and $y$ dimensions. Note that if, for example, the integrand is $2\pi$
periodic in 
$x$ and $4\pi$ periodic in $y$, the box must be of side
$4\pi$. Additionally, for one routine, we assume that the integrand
is described by $n \times n$ equally spaced particles, with $n$
particles in each dimension. 

\section{Using the PMTA Routines}

\subsection{Compiling}
Our current version allows compilation on either a workstation or on a
CRAY/SGI shared memory parallel machine. Because the code is in an
early stage of development, we have separate listings for each
architecture.  This is because
each platform has its own unique multiprocessing compiler directives. 
We should note that the driving routines are essentially the same. It
is the PMTA 
library routines that differ from platform to platform. 
We may later incorporate all of the changes into a single code. The
compilation of the code is 
similar for both versions. The PMTA routines should be compiled using
an ANSI C compiler. We have had success using {\tt gcc} on Sun/Sparc
platforms, and using {\tt cc} on SGI/CRAY platforms. 

It is important that the variables SERIAL and SURFACE are defined to
the preprocessor. We use the compile options -DSERIAL and
-DSURFACE. We include a Makefile which demonstrates how to compile
programs. 

\subsection{PMTA routines}
The use of the PMTA routines involves three steps - Initialization,
Calculation, and cleanup. The call to initialize the PMTA routines is 
InitializeCoulomb() with a variety of calling options (described
later). The next step is to call Coulomb(). The final step is
clean up after all is done. We now describe the variables to these
three routines. 

The information about the surface is contained in the variable {\tt
ParticleTable} which is allocated by:
\begin{verbatim} 
   ParticleTable = (PTablePtr) malloc(NumParticles * sizeof(PTable));} 

   for(i=0;i<NumParticles;i++){
       ParticleTable[i].p.x = 
       ParticleTable[i].p.y = 
       ParticleTable[i].p.z =      <----- Set to appropriate values
       ParticleTable[i].p.q1 = 
       ParticleTable[i].p.q2 = 
       ParticleTable[i].p.q3 =
    }

\end{verbatim} 
where {\tt NumParticles} is the total number of particles on the
interface, $x$, $y$, and $z$ are the particle positions, and q1,q2,and
q3 are the components of $Q$ in \ref{defeq}.

The routine which initializes is called as follows:
\begin{verbatim}

InitializeCoulomb(
      CubeLength,         /* Length of Cube                                      */
      NumLevels,          /* Number of tree levels: 3 -> 0,1,2 -> 8^2 leaves     */
      FFT,                /* FFT: 0=no 1=yes                                     */
      Mp,                 /* Number of terms, p                                  */
      PFMA+EPMTA,         /* Algorithm: 0:B&H 1:FMA 2: EPMTA (includes B&H & FMA)*/
      Theta,              /* B&H Theta  (0.5 - 4.0) required for Algorithm 0 & 2 */
      NumProcessors,      /* Number of Processors: 1 for serial                  */
      Direct,             /* O(n^2) direct computation 0=no 1=yes                */
      NumParticles,       /* NumParticles                                        */
      BALANCED,           /* Use load balancing: 0=no 1=yes                      */
      Verbose,            /* Verbose: 0=no 1=yes                                 */
      Pmon,               /* Pmon performance: 0=no 1=yes (ignored for serial    */
      ParticleTable,      /* ParticleTable                                       */
      Mk,                 /* Compute 3^k X 3^k X 3^k region of unit cells        */
      SphereRadius,       /* macroscopic sphere radius */
      EllipseX,           /* macroscopic ellipse dimensions */
      EllipseY,
      EllipseZ
    );

\end{verbatim}

Before calling this routine, the {\tt ParticleTable} should be
initialized to the starting values. We now describe the variables.
\begin{verbatim}

VARIABLE        DESCRIPTION AND RANGE
(type)

CubeLength      The length of the center cube, centered at
(double)        the origin, containing all of the particles. 

NumLevels       The number of levels in the oct tree. 
(int)           This value should be greater than 1.

FFT             To use FFT version - we always set this to zero.
(int)

Mp              The number of terms in the multipole expansion.
(int)           This value should be greater than zero.

PFMA            Specifies the Fast Multipole Algorithm. Should always
(int)           have the value 1.

EPMTA           We don't use this. Should be 0.
(int)

Theta           We don't use this. Should be 0.0.
(double)

NumProcessors   In the serial version, set to 1. On the parallel version
(int)           set to some power of 2 (e.g. 1,2,4,8,...). This is NOT
                the number of processors the machine will use. It
                breaks the work load into NumProcessors pieces. This
                number should generally be somewhat larger than the
                number of processors you are using because the pieces
                aren't load balanced. The parallel versions
                dynamically schedule the processor load. 

Direct          We don't use this. Should always be zero. 
(int) 		

NumParticles    The number of particles on the surface.
(int)

BALANCED        We don't use this. Set it to 1.
(int)

Verbose         Prints out a bunch of stuff. Set it to 0.
(int)

Pmon            We don't use this. Set it to 0.
(int)

ParticleTable 	The table of particles, previously defined.
(PTable)

Mk              Set up periodic calculations. Should always be 1.
(int)

SphereRadius, EllipseX, EllipseY, and EllipseZ are not used, and
should be less than 0 in value.
\end{verbatim}
A typical call to this routine, after defining {\tt NumLevels},{\tt
Mp}, {\tt NumParticles}, and {\tt ParticleTable} would look like

\begin{verbatim}

	InitializeCoulomb(CubeLength, NumLevels, 0, Mp, 1, 0.0,
                    1, 0, NumParticles, 1, 0, 0,
                    ParticleTable, 1, -1.0,-1.0,-1.0,-1.0);

\end{verbatim}

Admittedly it seems a little odd to keep features that are never used,
but this seemed to be the easiest way to write the code. Perhaps a
future version will be simplified. 

The actual calculations are computed by calling the routine
\begin{verbatim}
 Coulomb(ParticleTable);
\end{verbatim}
The value of the integral at the grid point {\tt i} is returned in
three parts which must be summed. The parts are 
\begin{verbatim}
	ParticleTable[i].f.x
	ParticleTable[i].f.y
	ParticleTable[i].f.z
\end{verbatim}

The sum of these three terms is an approximation to the following sum
(labeled $I_A$):
\begin{equation}
I_A({\bf x_i})=\sum_{\bf i\ne j} Q({\bf x_j}) \cdot \nabla_{\bf
x}G({\bf x_i-x_j}).
\end{equation}
The user should note that this quantity should then be multiplied by
the mesh spacing. Using this form, it is possible to choose different
mesh spacings for each spatial dimension.

When all calculations are finished, the routine 
\begin{verbatim}
FreeCoulomb();
\end{verbatim}
frees the memory allocated. NOTE: On the SGI Indy platform, this call
sometimes results in a Segmentation Error. We have not encountered
this difficulty on other platforms.

\subsection{Original Subroutines}

We have added a number of original subroutines to the code in order to
maximize efficiency. (Not all of these routines are essential to using
the code.) 
These routines are readcoef(), coulomb(), adjust\_mem(),
align(),dnumlev(),zero\_all(), align() and  inumlev().

We here describe all but coulomb():
\begin{verbatim}

ROUTINE                 DESCRIPTION

readcoef(factor)   This routine reads and scales the coefficients
                   used for  
                   computing the far field contribution. The
                   argument factor should be set to 
                   factor = 2.0*pi/CubeLength;

adjust_mem()       This routine takes care of the dynamic memory
                   allocation. This should be called after each
                   time step before doing any more calculations.

zero_all()         This routine resets certain values in the oct 
                   tree. It is a housecleaning routine.

align()            Permits changing the number of particles
                   without reinitializing the program.

dnumlev()          These routines permit changing the number of 
inumlev()          oct tree levels without reinitializing	
                   the program.

\end{verbatim}

The last four routines are used by the routine coulomb(), and are not
essential to the operation of the code.  However, the first two are
essential to our adaptation.

The subroutine coulomb() is used to compute the integral for either
first or third order accuracy. For a third order accurate calculation,
the integral must be computed with two different mesh spacings. This
requires 5 calls to the subroutine Coulomb(). For a fixed number of
terms in the multipole expansion, the best value of NumLevels differs
for different values of NumParticles. In writing coulomb(), we use the
subroutines align, dnumlev(), and inumlev() to take advantage of this.

To see how these are used, we include two segments from the routine
coulomb() 
\begin{verbatim} 
...
  n/=2;
  NumParticles/=4;    <---- decrease the number of particles
  NumLevels--;        <---- decrease the number of levels 
  align();            <---- housekeeping for particle number
  dnumlev();          <---- housekeeping for level number

  for(i=0; i< NumParticles; i++){
    xi=i/n;
    xj=i%n;
    xi*=2;
    xj*=2;
    ParticleTable[i].p.x = x[xi];
    ParticleTable[i].p.y = y[xj];
    ParticleTable[i].p.z = eta[xi][xj];
    ParticleTable[i].q1 = q1[xi][xj];
    ParticleTable[i].q2 = q2[xi][xj];
    ParticleTable[i].q3 = q3[xi][xj];    
  }
  zero_all();         <---- reset values in the oct tree
  ParticleCells();    <---- distribute particles in the oct tree
  CellPlists();       <---- housekeeping in the oct tree

  Coulomb(ParticleTable);
...
  NumParticles*=4;
  NumLevels++;
  n*=2;
  align();
  inumlev();

  for(i=0; i< NumParticles; i++){
    xi=i/n;
    xj=i%n;
    ParticleTable[i].p.x = x[xi];
    ParticleTable[i].p.y = y[xj];
    ParticleTable[i].p.z = eta[xi][xj];
    ParticleTable[i].q1 = q1[xi][xj];
    ParticleTable[i].q2 = q2[xi][xj];
    ParticleTable[i].q3 = q3[xi][xj];    
  }
  zero_all();
  ParticleCells();
  CellPlists();

  Coulomb(ParticleTable);
...	
\end{verbatim}
In the first segment, we want to perform calculations at the coarser
mesh scale. We divide the number of particles by 4, decrease the
number of levels, do some housekeeping and then do the computations. 
When we want to compute at the finer level, we reverse the process.

The call to the subroutine coulomb has the form
\begin{verbatim}
coulomb(x,y,phi,z,q1,q1,q1,ParticleTable,n,ord);
\end{verbatim}
Here x,y, and z are the coordinate positions of the particles. The
variables x and y are either arrays or vectors, depending on whether
the fixed grid or Lagrangian method is being used. The variables
qx1, qy1, and qz1 are arrays which contain the discrete values of the
components of $Q$ in Equation \ref{defeq}. ParticleTable has been defined
previously. The variable n is the number of particles in each
dimension. The variable ord takes the value 1 or 3 for first or third
order quadrature. The variable phi is an array containing the values
of the integral at the mesh points. The arrays x, y, eta, q1, q2, and
q3 must be arrays of dimension $n\times n$. 

The values returned by the routine Coulomb() are scaled in the routine
coulomb() by a factor of $h^2/4\pi$ where $h={\tt
CubeLength}/n$. Hence, the routine Coulomb() computes the following
sum (denoted by $I_C$):
\begin{equation}
I_C({\bf x_i})={h^2 \over 4\pi}\sum_{\bf i\ne j} Q({\bf x_j}) \cdot
\nabla_{\bf 
x}G({\bf x_i-x_j}).
\end{equation}

\subsection{Dynamic Memory Allocation}

At the end of each time step, the particle positions and charge should
be updated. Then the routines ParticleCells() and adjust\_mem() should
be called. This allocates the memory necessary for the next time step.
The code should look something like
\begin{verbatim}
...
for(i=0; i< NumParticles; i++){
      xi=i/n;
      xj=i%n;
      ParticleTable[i].p.x = x[xi][xj];
      ParticleTable[i].p.y = y[xi][xj];
      ParticleTable[i].p.z = z[xi][xj];
    }
    ParticleCells();
    adjust_mem();
...
\end{verbatim}

\subsection{Driving Routines}
In order to make the routines more accessible, we have included
several simple driving programs which demonstrate the use of the
routines. The implementation should be something like the following:
\begin{verbatim}

#include ``PMTA.h''
main()
{
       declarations and assignments;
     
       InitializeCoulomb(arguments);
 
       for(timestep=1;timestep<=NN;timestep++){
           
          Coulomb(arguments);
or
          coulomb(arguments);

          update variables;

	  for(i=0; i< NumParticles; i++){
               xi=i/n;	
               xj=i%n;
               ParticleTable[i].p.x = x[xi][xj];
               ParticleTable[i].p.y = y[xi][xj];
               ParticleTable[i].p.z = z[xi][xj];
           }
           ParticleCells();
           adjust_mem();

       } /* time steps */

       FreeCoulomb();
       
}
\end{verbatim}

It is important to note that the particles must remain in the cell
centered at the origin. It is the users responsibility to verify that
this is so. The code does not currently check to verify this. 

We have included several sample codes which were used for some of our
earlier experiments. It should be observed that these are only samples
and should not be considered the only approaches to using the
routines. In particular, we have made use of Numerical Recipes
routines for computing spatial derivatives. While these are acceptable
for our present studies, there are more efficient FFT routines
available for specific high performance architectures. We include one
routine which computes waves for a fixed grid, and one which computes
waves for a Lagrangian grid. The programs are, respectively, wavesf.c
and wavesb.c, with supporting routines in wavlibf.c and wavlibb.c.

\section{Additional programs}
In addition to these basic routines, we include descriptions of
several other programs which may be of interest when using the
WaVeMaKeR routines.

The first of these computes the coefficients required for the periodic
calculations. This is a Fortran program which should only need to be run
once. (The coefficients are contained in the file {\tt coef1} so it
should not be necessary to run this program.) 
We have labeled this program percoefcalc1.f. All of the
necessary subroutines are included in the source file. This program
should only take a few minutes to run on a workstation. An appropriate
value of the variable $p$ is 100.

The program eqorth31.c is written in C and computes an equal
orthogonal coordinate system for a 
given surface. Note that an explicit form of the surface and its
derivatives is necessary for this to work. This program is still in
fairly rough form. However, the basic idea should be evident. We have
used this for waves of amplitude $.5$. For higher amplitudes, the
results might not be satisfactory. The derivative routines and Poisson
solver use the numerical recipes routines nrutil.c rlft3.c
fourn.c. Because of copyright restrictions, we cannot distribute these
routines. 

The next collection of programs computes the smoothing factors
required for performing the Fourier smoothing. This set of programs is
still in development, and again, is provided primarily for
informational purposes. The program collection includes the following:
frho311.f coeffdri1.f big01.f big0.f big1.f big2.f big3.f
mkrscompinc.f mkcompinc.f frssumd.f fsumd.f. 

The program frho311.f computes the smoothing factors for the third
order quadrature rule, and coeffdri1.f
computes the coefficients of the series used by frho311.f. 
The factors computed are the following:
\begin{eqnarray}
\rho^1({\bf k}h,\alpha)&=&{1\over 2\pi}\sum_{\bf j\ne
0}{j_1\sin(j_1k_1h)\cos(j_2k_2h) \over (j_1^2+\alpha j_2^2)^{3/2}} \\
\rho^2({\bf k}h,\alpha)&=&{1\over 2\pi}\sum_{\bf j\ne
0}{j_2\sin(j_2k_2h)\cos(j_1k_1h) \over (j_1^2+\alpha j_2^2)^{3/2}} \\
\bar{\rho}({\bf k}h,\alpha)&=&{1\over 2\pi |{\bf k}|}\sum_{\bf j\ne
0}{1-\cos(j_1k_1h)\cos(j_2k_2h) \over (j_1^2+\alpha j_2^2)^{3/2}} \\
\tilde{\rho}&=&{k_1\rho^1+k_2\rho^2 \over |{\bf k}|}
\end{eqnarray}
The output includes these variables, as well as the ratio
$\tilde{\rho}/\bar{\rho}$. The number of modes used is set in the
declarations at the beginning of the program. At the present, we
assume an equal number of modes in each spatial dimension.

We have included a
Makefile to use for compiling the program coeffdri1.f . It is
important to note that 
the big*.f files are written in optimized Fortran, and do not need to
be compiled using an optimizer. 
The program coeffdri1.x should be run first, followed by frho311.x. 
The value of $\alpha$ should be stored as the second line in the file
coeff.in. The first line should be 1.0. 
The programs coeffdri1.x and frho311.x should take only a few seconds
to execute on a workstation, despite the enormous size of the programs.
In its present form the calculations should have at worst single
precision accuracy. 

The last programs which we detail computes coefficients of three
dimensional steady waves. The program newtcont.f computes steady waves
using continuation in the wave height. The program newtcont1.f
computes the coefficients using continuation in the amplitude of the
mode $(3,1)$.



