describe recent code mods:

* works properly on multiprocessors
* sets excision mask (not properly described yet :( )


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@990 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 62 insertions, 50 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index bc7e415..69adc4f 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -152,16 +152,16 @@ outgoing null geodesics has zero expansion.  In terms of the usual
 $3+1$ variables, an apparent horizon satisfies the equation
 \begin{equation}
 H \equiv \del_i n^i + K_\ij n^i n^j - K = 0
-			\label{AEIDevelopment_AHFinderDirect_eqn-horizon}
+					\label{AHFinderDirect/eqn-horizon}
 \end{equation}
 where $n^i$ is the outward-pointing unit normal to the apparent horizon,
 and $\del_i$ is the covariant derivative operator associated with the
 3-metric $g_\ij$ in the slice.
 (See \cite{AEIDevelopment_AHFinderDirect_York-1989-in-Frontiers} for a
-derivation of~\eqref{AEIDevelopment_AHFinderDirect_eqn-horizon}.)
+derivation of~\eqref{AHFinderDirect/eqn-horizon}.)
 
 Thorn~\thorn{AHFinderDirect} finds an apparent horizon by numerically
-solving~\eqref{AEIDevelopment_AHFinderDirect_eqn-horizon}.  It requires
+solving~\eqref{AHFinderDirect/eqn-horizon}.  It requires
 as input the usual Cactus 3-metric $g_\ij$ and extrinsic curvature
 $K_\ij$, (and optionally the conformal factor $\psi$ if the
 \thorn{StaticConformal} metric semantics are used),
@@ -173,7 +173,7 @@ the irreducable mass associated with the area.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{What \thorn{AHFinderDirect} Needs}
-\label{AEIDeveopment_AHFinderDirect_sect-what-AHFinderDirect-needs}
+\label{AHFinderDirect/sect-what-AHFinderDirect-needs}
 
 There are some restrictions on the spacetime, or more precisely
 on each slice where you want to find apparent horizons, which are
@@ -215,14 +215,12 @@ necessary in order for \thorn{AHFinderDirect} to work:
 	(where the Cactus 3-D grid resolution is too low) is that
 	the interpolated $g_{ij}$ fails to be positive definite.
 \item	At the moment \thorn{AHFinderDirect} and the Cactus interpolators
-	don't know anything about excision, so if the apparent horizon
+	don't know how to avoid an excised region, so if the apparent horizon
 	(or any trial horizon surface as the algorithm is iterating
 	towards the apparent horizon) gets too close to an excised
 	region, you'll get garbage results as the interpolator tries
-	to interpolate data from the excised region.  I hope to
-	properly support excision (at least in the sense of checking
-	for the trial horizon surface being too close to an excised
-	region) soon.
+	to interpolate data from the excised region.  I plan to fix
+	this sometime soon.
 \item	\thorn{AHFinderDirect} requires that any apparent horizon
 	it's going to (try to) find must be a \defn{Strahlk\"{o}rper}
 	(literally ``ray body'', or more commonly ``star-shaped region'')
@@ -252,8 +250,9 @@ here's what works and what doesn't:
 \item	Obviously, your Cactus configuration must include
 	\thorn{AHFinderDirect}, and your \code{ActiveThorns} parameter(s)
 	must activate it.  \thorn{AHFinderDirect} inherits from \thorn{Grid},
-	\thorn{ADMBase}, and \thorn{StaticConformal}, so you'll need
-	those (or more precisely some thorns providing them), too.
+	\thorn{ADMBase}, \thorn{SpaceMask}, and \thorn{StaticConformal},
+	so you'll need those (or more precisely some thorns providing them),
+	too.
 \item	\verb|Grid::domain = "full"|, \verb|"bitant"|,
 	\verb|"quadrant"|, and \verb|"octant"| are supported.
 	Alas, at present rotating (or more precisely nonlocal)
@@ -270,22 +269,20 @@ here's what works and what doesn't:
 	are only provided by the thorns \thorn{PUGHInterp} and
 	\thorn{LocalInterp} respectively, so you need to have
 	these thorns compiled into your configuration and activated.
-\item	\thorn{AHFinderDirect} works in multi-processor Cactus runs,
-	but at present it's quite inefficient for large numbers of
-	processors.  In fact, it exhibits ``anti-scaling'' with the
-	number of processors: the run time is {\em larger\/} as the
-	number of processors increases.  I hope to fix this inefficiency
-	soon, but for now, the run time relative to the single-processor
-	case is roughly
-	\begin{flushleft}
-	\begin{tabular}{cc}
-	$N$		& $t_N / t_1$		\\
-	\hline
-	2		& 1.4			\\
-	4		& 2.0			\\
-	$N>4$		& $1.55 + 0.12N$	%%%\\
-	\end{tabular}
-	\end{flushleft}
+\item	\thorn{AHFinderDirect} works fine in single- or multi-processor
+	Cactus runs.
+\item	\thorn{AHFinderDirect} can set an excision mask based on
+	each horizon's shape, using either the old-style (\verb|CCTK_REAL|)
+	mask (compatible with \thorn{AHFinder}) or the new-style
+	(\verb|CCTK_INT|) mask bit-fields defined by \thorn{SpaceMask}.
+	However, at the moment the mask is only set based on the
+	actual patch system
+	(see section~\ref{AHFinderDirect/sect-multipatch-system})
+	\eg{} if \verb|grid::domain| and the local coordinate origin
+(see section~\ref{AHFinderDirect/sect-parameters/local-coordinate-origin})
+	are such that the patch system has octant symmetry, then the
+	mask is only set to excise that octant, not the full horizon
+	shape.  I consider this to be a bug, and I plan to fix it soon.
 \end{itemize}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -334,7 +331,7 @@ possibly set explicitly.
 Note that all of the ``\code{[}$n$\code{]}'' parameters are Cactus
 array parameters, which you need to specify separately in your parameter
 file for each apparent horizon.  The example in
-section~\ref{AEIDevelopment_AHFinderDirect_sect-examples} should make
+section~\ref{AHFinderDirect/sect-examples} should make
 this clear.
 
 %%%%%%%%%%%%%%%%%%%%
@@ -403,23 +400,28 @@ this clear.
 	\mbox{}\\
 		Also print a single line for each Newton iteration giving
 		the 2-norm and $\infty$-norm of the $H(h)$ function defined
-		by~\eqref{AEIDevelopment_AHFinderDirect_eqn-horizon}.
+		by~\eqref{AHFinderDirect/eqn-horizon}.
 		This is the default.
 	\item[\code{"algorithm details"}]
 	\mbox{}\\
 		Print lots of detailed messages tracing what the code
 		is doing.
+	\item[\code{"algorithm debug"}]
+	\mbox{}\\
+		Print even more detailed messages tracing what the code
+		is doing, mainly useful for debugging purposes.
 	\end{description}
 \end{description}
 
 %%%%%%%%%%%%%%%%%%%%
 
 \subsection{Choosing the Local Coordinate Origin for each Apparent Horizon}
+\label{AHFinderDirect/sect-parameters/local-coordinate-origin}
 
 For each apparent horizon you want to find, you need to specify
 the Cactus $(x,y,z)$ coordinates of a local coordinate system origin.
 As described in
-section~\ref{AEIDeveopment_AHFinderDirect_sect-what-AHFinderDirect-needs},
+section~\ref{AHFinderDirect/sect-what-AHFinderDirect-needs},
 each apparent horizon must be a Strahlk\"{o}rper with respect to its
 local coordinate system origin.
 
@@ -462,7 +464,7 @@ I hope to add such a provision soon.%%%
 \thorn{AHFinderDirect} requires an initial guess
 for the apparent horizon's coordinate position and shape
 (that is, for the $h(\text{angle})$ function defined in
-section~\ref{AEIDeveopment_AHFinderDirect_sect-what-AHFinderDirect-needs}),
+section~\ref{AHFinderDirect/sect-what-AHFinderDirect-needs}),
 for each apparent horizon you want to find.  Unlike some other apparent
 horizon finders (\eg{}~the curvature flow method in \thorn{AHFinder}),
 for \thorn{AHFinderDirect} there's no restriction on whether the initial
@@ -798,7 +800,7 @@ at some point this thorn may be enhanced to do this.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Examples}
-\label{AEIDevelopment_AHFinderDirect_sect-examples}
+\label{AHFinderDirect/sect-examples}
 
 There are a few example parameter files in the \code{par/} directory,
 including Kerr initial data, Misner initial data, and Misner time-evolution
@@ -868,7 +870,7 @@ AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{How \thorn{AHFinderDirect} Works}
-\label{AEIDevelopment_AHFinderDirect_sect-how-ahfinderdirect-works}
+\label{AHFinderDirect/sect-how-ahfinderdirect-works}
 
 \thorn{AHFinderDirect} uses the apparent horizon (henceforth \defn{horizon})
 finding algorithm of
@@ -882,7 +884,7 @@ grid:
 
 As described above, I parameterizes the horizon shape by
 $r = h(\text{angle})$ for some single-value function $h: S^2 \to \Re^+$.
-The apparent horizon equation~\eqref{AEIDevelopment_AHFinderDirect_eqn-horizon}
+The apparent horizon equation~\eqref{AHFinderDirect/eqn-horizon}
 then becomes a 2-D elliptic PDE on $S^2$ for the function $h$.  I finite
 difference this in angle to obtain a system of simultaneous nonlinear
 algebraic equations for $h$ at the angular grid points, and solve this
@@ -918,27 +920,38 @@ Computationally, this algorithm has 3 main parts:
 	described in
 \cite{AEIDevelopment_AHFinderDirect_Thornburg-1996-apparent-horizon-finding},
 	so this computation is quite fast.  The Jacobian is a highly
-	sparse matrix, but at the moment we neglect this and store it
-	as a full (dense) matrix.  (This should change soon.)
+	sparse matrix; \thorn{AHFinderDirect} has code to store it
+	as either a dense matrix (for debugging purposes), or a sparse
+	matrix (the default).  Which option is used is determined by
+	a compile-time configuration in \verb|src/include/config.h|.
 \item	Solving the nonlinear equations $H(h) = 0$ by a global Newton's
-	method or a variant.  At the moment this thorn uses a simple
-	Newton's-method solver, and solves the linear system using
-	$\sf LU$ decomposition routines from LAPACK
-	(this should change to a sparse $\sf LU$ decomposition/solution
-	soon).
+	method or a variant.  How this is done depends on how the Jacobian
+	is stored.  At present,
+	\begin{itemize}
+	\item	If \thorn{AHFinderDirect} is configured to store the
+		Jacobian as a dense matrix, then LAPACK is used to solve
+		the linear equations.
+	\item	If \thorn{AHFinderDirect} is configured to store the
+		Jacobian as a sparse matrix, then an
+		incomplete-$\sf LU$-decomposition--conjugate-gradient
+		solver is used.
+	\end{itemize}
+	By default only the sparse-matrix code is configured, so LAPACK
+	isn't used and there's no need to link with the LAPACK library.
 \end{itemize}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{The Multipatch System}
+\label{AHFinderDirect/sect-multipatch-system}
 
 Perhaps the most unusual feature of \thorn{AHFinderDirect} is the
 ``multipatch'' system used to cover $S^2$ without coordinate singularities.
 In general there are 6~patches, one each covering a neighborhood of
 the $\pm z$, $\pm x$, and $\pm y$ axes, but this may be reduced in
 the presence of suitable symmetries.  For example,
-figure~\ref{AEIDevelopment_AHFinderDirect_fig-3patch}
-on page~\pageref{AEIDevelopment_AHFinderDirect_fig-3patch} shows a system
+figure~\ref{AHFinderDirect/fig-3patch}
+on page~\pageref{AHFinderDirect/fig-3patch} shows a system
 of 3~patches covering the $+xyz$~octant of $S^2$.  This would be
 suitable for finding an apparent horizon with mirror symmetry about
 the (local) $z=0$~plane, and either 90~degree periodic rotation symmetry
@@ -956,7 +969,7 @@ $x$~and $y$~axes.
 	The angular resolution is 5~degrees.  Notice that the
 	patches overlap by several ``ghost zone'' grid points.
 	}
-\label{AEIDevelopment_AHFinderDirect_fig-3patch}
+\label{AHFinderDirect/fig-3patch}
 \end{figure}
 
 To allow easy angular finite differencing within the patch system,
@@ -967,10 +980,9 @@ each patch is extended beyond its nominal extent by a ``ghost zone''%%%
 	 would call these ``patch zones'' or ``symmetry
 	 zones''.
 	 }%%%
-{} (2~grid points wide in
-figure~\ref{AEIDevelopment_AHFinderDirect_fig-3patch}).  Angular grid
-function values in the ghost zone can be obtained by interpatch
-interpolation%%%
+{} (2~grid points wide in figure~\ref{AHFinderDirect/fig-3patch}).
+Angular grid function values in the ghost zone can be obtained by
+interpatch interpolation%%%
 \footnote{%%%
 	 Due to the way the patch coordinates are defined,
 	 adjacent patches always share a common ``perpendicular''
@@ -982,8 +994,8 @@ finite differencing within the nominal extent of each patch can proceed
 normally, ignoring the patch boundaries.  \thorn{AHFinderDirect} can
 be configured at compile time to use either 2nd~order or 4th~order
 angular finite differencing (3~point or 5~point angular molecules);
-the default is 4th~order (5~point).  See \verb|src/include/config.h|
-if you want to change this.
+the default is 4th~order (5~point).  This is configured at compile time
+in \verb|src/include/config.h|.
 
 By default \thorn{AHFinderDirect} will automagically choose a patch
 system type for each apparent horizon searched for, based on the local