 \documentstyle[12pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P A G E   L A Y O U T . . . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textwidth6.25in
\textheight8.5in
\evensidemargin0in
\oddsidemargin0in
\topmargin-.5in
\parindent0em
\parskip=12pt
\def\baselinestretch{1.2}
\pagestyle{empty}
\input{psfig}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\landscape
%\special{landscape}
\begin{document}
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% F O N T S . . .   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M A C R O S . . . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\BEGINSLIDE}[1]{\vspace*{#1}}
\newcommand{\ENDSLIDE}[1]{\vspace*{#1}\pagebreak}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
{\Large
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B E G I N    S L I D E S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.0in}
\section{Binary Black Hole Grand Challenge} 
\begin{center}
	 {\bf Black hole collisions have 2 major interest to numerical relativist }
 \end{center}
\begin{enumerate}
\item Black hole collisions should be strong generators of gravitational
      radiation.  LIGO(Laser Interferometric Gravitational Wave Observatory)
      should be able to detect this radiation by the end of the century.
\item This is an extremely difficult problem for numerical relativist.  
		We will make use of the most powerful computers and the sophistication
      of AMR.
      It should provide us with the first ``general'' solution of Einstein
      Equations.
\end{enumerate}
\begin{center} Groups involved \end{center}
\begin{itemize}
	\item Physics
		\begin{itemize}
			\item {\bf Austin},
               {\bf Cornell},
               {\bf Illinois},
			      {\bf North Carolina},
			      {\bf Northwestern},
			      {\bf Penn State}, 
			      {\bf Pittsburgh},
		\end{itemize}
	\item Computer Science
		\begin{itemize}
			\item {\bf Austin},
			{\bf Illinois}
			{\bf Syracuse}
		\end{itemize}
\end{itemize}
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.0in}
\section{4 dimensional spacetime into 3+1}
\begin{itemize}
\item The four dimensional spacetime is split into 3 dimensional hypersurfaces
which evolve in time.
\item Breaking spacetime up this way gives us 10 non-linear PDE's to solve
\item These equations are divided into two classes.
	\begin{itemize}
		\item 4 Initial Value Equations (Elliptic)
		\item 6 Coupled Hyperbolic equations (Hyperbolic)
	\end{itemize}
\end{itemize}
One of the most important components in the solution to the
               evolution equation is the choice of gauges.
\begin{itemize}
	\item The choice of gauge is the choice of coordinate systems.
	\item Some gauges can allow all the computational grids points to be
         ``sucked'' inside a black hole, where as other gauge choices can
         allow the hole to move on the computational domain.
\end{itemize}

{\bf NOTE:}
In Classical Physics evolutions are dictated by momentum variables and
position variables.  Numerical Relativity uses the extrinsic curvature, $K_{ij}$
to act like the momentum, and the metric $g_{ij}$ to act as the position.
\ENDSLIDE{0.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.0in}
\section{Initial Value Problem for 2 black holes}
\begin{itemize}
	\item The 4 equations break down into two sets of equations
		\begin{itemize}
			\item The Momentum Constraint Equation ( vector equation)
			\begin{equation}
  			\bar{D}_j\bar{K}^{ij} = 0.
			\end{equation}
			\item The Hamiltonian Constraint Equation (scalar equation)
			\begin{equation}
  			\bar{\nabla}^2\psi = -\frac{1}{8}\psi^{-7}\bar{K}_{ij}\bar{K}^{ij}
			\end{equation}
		\end{itemize}
\end{itemize}
The momentum equation can be solved with semi-analytical techniques for 
2 black holes.
The Hamiltonian Constraint equation solver solves this constraint equation
along with the boundary conditions:

The boundary condition at infinity follows directly from asymptotic
 and conformal flatness:
\begin{equation}
  {\partial\psi\over\partial r} = {1 - \psi\over r},
\end{equation}

The boundary condition around the holes is, ($\alpha=1,2$)
\begin{equation}
   n^i_\alpha\bar{D}_i\psi\biggl|_{B_\alpha} =
         -\frac{\psi}{2r_\alpha}\biggr|_{B_\alpha},
\end{equation}
where $r_\alpha = |\vec{x} - \vec{C}_\alpha|$ and
$\vec{n}_\alpha = (\vec{x} - \vec{C}_\alpha)/r_\alpha$ (these are
Euclidean operators in the flat background space).

The first 3D solutions were found by three separate groups
\begin{enumerate}
	\item G. Cook   (Cornell University)
	\item S. Klasky, M. Choptuik, R. Matzner (University Texas at Austin) 
	\item M. Dubal, S. Oliveira, R. Matzner (University Texas at Austin) 
\end{enumerate}
\begin{center} {\bf  This problem has been solved} \end{center}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.0in}
\section{Methods to try to solve the Evolution Equations}
\begin{itemize}
\item There are many methods to try to solve the evolution equations.
	\begin{itemize}
		\item ADM formalisms
		\item Hyperbolic formalisms
	\end{itemize}
\end{itemize}

In order to solve these equations we must use AMR.  This is because
the problem varies over many length scales.  

Most difficult part of the solution will be the boundary
conditions.  

The outer boundary problem will be solved by using a matching condition
between a Cauchy(evolution) code, and a characteristic, or perturbation
code. 

The Cauchy code can ``handle'' the strongly nonlinear regime, and
the Characteristic code can ``handle'' only the weakly nonlinear regime.
The boundary condition is only known at infinity, thus we need two codes.

\ENDSLIDE{-0.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section{Methods to try to solve the Evolution Equations}
The inner boundary is a problem for several reasons:
	\begin{enumerate}
		\item We will excise the points inside the apparent horizon.
		\item We would like to excise the points inside the event horizon.
            The problem is that the event horizon is only known at $t=\infty$.
            Thus, we need a way to exclude the singularity of the hole.
		\item The apparent horizon is found by looking for {\em Trapped surfaces}.
				This breaks down into solving a very difficult elliptic equation.
		\item When the hole moves in the computational domain, coordinate points
            that were previously excised will now appear
            inside the computational
            domain.
	\end{enumerate}
Many groups have explored special ways to difference around the holes
(causal differencing).
A new formalism, J. York, the ``hyperbolic'' formalism was developed such that
relativist would not have to {\em worry} about the gauge freedom. 
\ENDSLIDE{ 0.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section{The Texas, Syracuse  {\tt T2} approach }

\begin{center}
	Here we describe one such ADM method {\tt T2}
   that Texas and Syracuse are working on.
\end{center}

The evolution equations that the {\tt T2} solves are as follows:
\begin{eqnarray}
\partial_t g_{ij} & = & - 2 \alpha K_{ij} + D_i \beta_j + D_j \beta_i \\
\label{evolk.eqn}
\partial_t K^i_j  & = & {\cal L}_{\beta} K^i_j - D^i D_j \alpha +
\alpha \left( R^i_j + K K^i_j \right)
\label{evol.eqn}
\end{eqnarray}
where $g_{ij}$ is the 3-metric,

$K^i_j$ is the extrinsic curvature,
$D_i$ is the covariant derivative computational with $g_{ij}$,

 $\alpha$ is the lapse function, $\beta^i$ is the shift vector 

and ${\cal L}_{\beta}$ is the Lie
derivative along the vector field $\beta$. 

${\cal L}_{\beta} K^i_j$ is given
by the following equation:
\begin{equation}
{\cal L}_{\beta} K^i_j = \beta^k D_k K^i_j - K^k_j D_k \beta^i
+ K^i_k D_j \beta_k
\label{lie.term}
\end{equation}
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section{Variable Storage}
{\large
\begin{tabular}{ll}
\hline
Grid Function & Number of independent components \\ \hline
$g_{ij}$  & $6\times 2$ independent components(two = n+1 timestep) \\
$g^{ij}$  & 6 independent components(required for various stages) \\
$K^i_j$   & $9\times 2$ components (multiply by two for n+1 timestep) \\
$\Gamma^i_{jk}$ & 18 components (stored instead of $g_{ij,k}$) \\
$R_{ij}$  & 6 components (Raised when needed in eqn.~(\ref{evolk.eqn})) \\
$\alpha$  & scalar function \\
$\beta^i$ & 3 components ($\beta_i$ can be obtained using $g_{ij}$) \\
\hline
\end{tabular}
}
Hence the lower limit on required variables for the solution to the evolution
is $64$ grid variables.
Plans are in the work for writing the equations using
$K_{ij}$ thus eliminating 6 variables.
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section{Modules for the connection coefficients}
The calculation of the connection coefficients can be broken up into several
stages. This obviates the necessity for storing the first derivatives of the
3-metric which saves us $18$ grid functions.

The connection coefficients are given by
\begin{equation}
\Gamma^i_{jk} = \frac{1}{2} g^{il} \left( g_{lj,k} + g_{lk,j} - g_{jk,l}
\right)
\label{connect.1}
\end{equation}

The calculation can be broken up into three separate sets using $9$ temporary
variables.

There are 3 steps:
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\vskip -0.5in
\section{3 Steps for Connection coefficients}
{\large
For the first term in eqn.~(\ref{connect.1})
\begin{eqnarray}
\Gamma^i_{jk} & = & \frac{1}{2} g^{il} g_{lj,k} \nonumber \\
& = & \frac{1}{2} \left(g^{i1} g_{1j,k} + g^{i2} g_{2j,k} + g^{i3} g_{3j,k}
\right) \label{connect.2}\\
\end{eqnarray}
We call these parts C1-C3.
For the second term in eqn.~(\ref{connect.1})
\begin{eqnarray}
\Gamma^i_{jk} & = & \Gamma^i_{jk} + \frac{1}{2} g^{il} g_{lk,j} \nonumber \\
& = & \Gamma^i_{jk} + \frac{1}{2} \left(g^{i1} g_{1k,j} + g^{i2} g_{2k,j}
+ g^{i3} g_{3k,j} \right) \label{connect.3}\\
\end{eqnarray}
We call these parts C4-C6.
For the third term in eqn.~(\ref{connect.1})
\begin{eqnarray}
\Gamma^i_{jk} & = & \Gamma^i_{jk} - \frac{1}{2} g^{il} g_{jk,l} \nonumber \\
& = & \Gamma^i_{jk} - \frac{1}{2} \left(g^{i1} g_{jk,1} + g^{i2} g_{jk,2}
+ g^{i3} g_{jk,3} \right) \label{connect.4}
\end{eqnarray}
We call these parts C7-C9.
}
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section {Ricci tensor}
The calculation of the Ricci tensor can be similarly broken up into several
segments with the use of $9$ temporary grid functions once more. The Ricci
tensor is given by
\begin{equation}
R_{ij} = \Gamma^a_{ij,a} - \Gamma^a_{ia,j} +
\Gamma^a_{ab} \Gamma^b_{ij} - \Gamma^a_{bj} \Gamma^b_{ia}
\label{ricci.1}
\end{equation}
The evaluation of the Ricci tensor
is broken as follows:
\begin{eqnarray}
R_{ij} & = & \Gamma^a_{ij,a} \nonumber \\
& = & \Gamma^1_{ij,1} + \Gamma^2_{ij,2} + \Gamma^3_{ij,3}
\label{ricci.2}
\end{eqnarray}
\begin{eqnarray}
R_{ij} & = & R_{ij} - \Gamma^a_{ia,j} \nonumber \\
& = & R_{ij} - \left(\Gamma^1_{i1,j} + \Gamma^2_{i2,j} + \Gamma^3_{i3,j}
\right)
\label{ricci.3}
\end{eqnarray}
\begin{eqnarray}
R_{ij} & = & R_{ij} + \Gamma^a_{ab} \Gamma^b_{ij} - \Gamma^a_{bj} \Gamma^b_{ia}
\label{ricci.4}
\end{eqnarray}
where in eqns.~(\ref{ricci.1},\ref{ricci.2}) terms like $\Gamma^1_{ij,1}$ and
$ \Gamma^1_{i1,j}$ are calculated using temporary arrays.
\ENDSLIDE{0.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.0in}
\section{Evolution equation for the 3-metric}
\begin{eqnarray}
\partial_t g_{ij} & = & - 2 \alpha K_{ij} + D_i \beta_j + D_j \beta_i. 
\end{eqnarray}

The evaluation of the second and third terms is carried out by once again
using $9$ temporary grid functions. Recall that $\beta^i$, the contravariant
form of the shift, is stored. Since $g_{ij}$ is compatible with $D_i$ we can
rewrite $D_i \beta_j = g_{jk} D_i \beta^k$. Hence we are left with calculating
\begin{eqnarray}
\label{gshift.eqn}
D_i \beta_j = g_{jk} \left( \partial_i \beta^k - \Gamma^k_{il}\beta^l \right).
\end{eqnarray}

Here once again we use a $9$ grid functions to store the components of
$\partial_i \beta^k$ first and then add $\Gamma^k_{il}\beta^l$ to form
eqn.~(\ref{gshift.eqn}). Then evaluating eqn.~(\ref{evolg.eqn}) becomes a
matter of using the temporary variables for $D_i \beta_j$.

Use second order finite differencing for time derivative.
5 Modules for these routines
{\large
\begin{verbatim}
C in evolve_g1
      call diff(nx,ny,nz,h,beta1,t11,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,t12,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,t13,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t21,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t22,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t23,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t31,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t32,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t33,'z',ml,dl,nl)
C in evolve_g2
      FORALL (i=1:3,j=1:3)
         t(i,j) = t(i,j)+c(1,i,j)*beta(1)+c(2,i,j)*beta(2)
     .            +c(3,i,j)*beta(3)
      END FORALL
      FORALL (i=1:3,j=1:3)
         t(i,j) = t(i,1) * g(1,j) + t(i,2) * g(2,j) + t(i,3) * g(3,j)
      END FORALL
C in evolve_g3
      FORALL (j=1:3,i=1:3,i<=j)
         t(i,j) = t(i,j) + t(j,i) 
         t(j,i) = t(i,j)
      END FORALL
C in evolve_g4
      FORALL (i=1:3,j=1:3)
         t(i,j) = t(i,j) - 2.0d0*alpha(1) * ( g(1,i) * k(1,j) +
     .            g(2,i) * k(2,j) + g(3,i) * k(3,j)  )
      END FORALL
C in evolve_g5
      FORALL (j=1:3,i=1:3,i<=j)
         g(i,j,1) = 2.0d0*dt*t(i,j) + g(i,j,-1)
      END FORALL
\end{verbatim}
}
\ENDSLIDE{-1.00in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.5in}
\section{Evolution equation for $K^i_j$}
\vskip -1in
\begin{eqnarray}
\nonumber
\partial_t K^i_j  & = & {\cal L}_{\beta} K^i_j - D^i D_j \alpha +
\alpha \left( R^i_j + K K^i_j \right)
\end{eqnarray}
We break this equation up into 5 parts(5 modules)
{\large
\begin{verbatim}
C evolve_k1
      call diff(nx,ny,nz,h,alpha,temp1,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,alpha,temp2,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,alpha,temp3,'z',ml,dl,nl)
C evolve_k2
      call diff(nx,ny,nz,h,temp1,t11,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,temp1,t12,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,temp1,t13,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,temp2,t22,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,temp2,t23,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,temp3,t33,'z',ml,dl,nl)
      t21 = t12
      t31 = t13
      t32 = t23
      FORALL (j=1:3,i=1:3,i<=j)
         t(i,j) = t(i,j) - c(1,i,j) *dalpha(1) -
     .              c(2,i,j)*dalpha(2) - c(3,i,j) * dalpha(3)
         t(j,i) = t(i,j)
      END FORALL
      call rais02i1() 
      FORALL (i=1:3,j=1:3)
         t(i,j) = -t(i,j)
      END FORALL
\end{verbatim}
}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.25in}
\section{Evolution equation for $K^i_j$ continued}
\large{
\begin{verbatim}
! Evolve K3
*     Lie_beta(K^1_1)
      call diff(nx,ny,nz,h,k11,temp1,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,k11,temp2,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,k11,temp3,'z',ml,dl,nl)
      t11 = t11 + beta1 * temp1 + beta2*temp2+beta3*temp3

      call diff(nx,ny,nz,h,beta1,temp1,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,temp2,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,temp3,'z',ml,dl,nl)
      t11 = t11 - k11 * temp1 - k21 * temp2 - k31 * temp3

      call diff(nx,ny,nz,h,beta1,temp1,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,temp2,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,temp3,'y',ml,dl,nl)
      t11 = t11 + k11 * temp1 + k12 * temp2 + k13 * temp3
C Evolve K4
      trace_k(1) = k11 + k22 + k33
      FORALL (i=1:3,j=1:3)
         t(i,j) = t(i,j) + alpha(1) *(
     .              gup(i,1) * r(1,j) + gup(i,2) * r(2,j) +
     .              gup(i,3) * r(3,j) + trace_k(1) * k(i,j) )
      END FORALL
C Evolve K5
      FORALL (i=1:3,j=1:3)
         k(i,j,1) = 2.0d0*dt*t(i,j) + k(i,j,-1)
      END FORALL
\end{verbatim}
}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.25in}
\section{Derivative routines, Vertex Centered}
Spatial derivative routines contain routines for Inner boundary, Outer boundary
and interior. Time derivatives, $O(h^2)$ use a 3 level leap frog scheme.
All derivative routines are second order for the first derivative in 
$x,y,z$ directions.

To handle the hole, 2 non-distributed lists are created, {\tt m(nl),d(nl)}
where {\tt m(l)} contains the {i,j,k} location of the ``hole'' points
and {\tt d(l)} contains either (1,-1,0) indicating either forward,
backwards or centered. 

The derivatives around the hole
are determined to second order.  Second derivatives, {\it for example $f_{xx}$,
are only first order accurate}. The interior second derivatives are 
second order accurate.
Routine looks as follows:
\large{
\begin{verbatim}
      subroutine diff(h,f,df,dir,m,d,nl)
      use params
      implicit none
      real*8 h
      real*8, dimension(nx,ny,nz) :: f,df
      real*8, dimension(3*nl)     :: d
      integer, dimension(3*nl)    :: m
      character dir
!HPF$ PROCESSORS PROC(nproc)
!HPF$ DISTRIBUTE f *(BLOCK,*,*) ONTO *PROC
!HPF$ ALIGN(:,:,:) WITH *f(:,:,:)         :: df
      call diff_int(h,f,df,dir)
      call diff_ob(h,f,df,dir)
      if (nl.gt.1) call diff_ib(h,f,df,dir,m,d,nl)
      return
      end
\end{verbatim}
}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.25in}
{\large
\begin{verbatim}
      subroutine diffx_ib(h,f,df,m,d,nl)
      use params
      implicit none
      integer nl
      real*8, dimension(nx,ny,nz) :: f,df
      real*8, dimension(3*nl)     :: d
      integer, dimension(3*nl)    :: m
      real*8 h,h2
      integer i,j,k,l

      h2 = 0.5d0 / h

      do l = 1 , nl
         i  = m(l)
         j  = m(nl+l)
         k  = m(2*nl+l)
         d1 = d(l)
         d2 = d(nl+l)
         d3 = d(2*nl+l)
         da = abs(d1)
         df(i,j,k) = h2*(d1*(-3.0d0*f(i,j,k)+4.0d0*f(i+int(d1),j,k)
     $                       -      f(i+2*int(d1),j,k)) +
     $               (1.0d0-da)*(f(i+1,j,k)-f(i-1,j,k) ) )
      end do
      return
      end
\end{verbatim}
}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{-1.25in}
{\large
\begin{verbatim}
      subroutine diffx_ob(h,f,df)
      use params
      implicit none
      real*8, dimension(nx,ny,nz) :: f,df
      real*8  h,h2,hi
      integer i,j,k
!HPF$ PROCESSORS PROC(nproc)
!HPF$ DISTRIBUTE f *(BLOCK,*,*) ONTO *PROC
!HPF$ ALIGN(:,:,:) WITH *f(:,:,:)         :: df
      hi = 1.0d0/h
      h2 = 0.5d0 / h
      FORALL (k=1:nz,j=1:ny)
         df(1,j,k) = h2 * ( -3.0d0*f(1,j,k)
     .             + 4.0d0*f(1+1,j,k) - f(1+2,j,k) )
         df(nx,j,k) = h2 * (  3.0d0*f(nx,j,k)
     .              - 4.0d0*f(nx-1,j,k) + f(nx-2,j,k) )
      END FORALL
      FORALL (i=2:nx-1,k=1:nz)
            df(i,1,k) = h2 *(f(i+1,1,k)-f(i-1,1,k) )
            df(i,ny,k) = h2 *(f(i+1,ny,k)-f(i-1,ny,k) )
      END FORALL
      FORALL (i=2:nx-1,j=1:ny)
            df(i,j,1) = h2 *(f(i+1,j,1)-f(i-1,j,1) )
            df(i,j,nz) = h2 *(f(i+1,j,nz)-f(i-1,j,nz) )
      END FORALL

      return
      end
\end{verbatim}
}
\ENDSLIDE{-1.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\BEGINSLIDE{0.0in}
\section{Future Plans} 
\begin{itemize}
\item Moving the hole on the domain, use some sort of Alcubierre and Schutz
method(implicit scheme using ADI). ( Huq \& Matzner)
\item Apparent Horizon solver written in a Data Parallel language.
(Huq, Klasky \& Matzner)
\item Causal differencing around the inner boundary. 
\item Stability of methods around the inner boundary is being looked at by
\item AMR in {\tt DAGH} and {\tt HPF} : Haupt,Huq,Klasky,Parashar
Correll
\end{itemize}
\ENDSLIDE{0.0in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E N D    S L I D E S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
}
\end  {document}
