This is a general code cleanup, and also adds some new features to

this thorn.

[Unforunately, for the reasons explained in today's E-mail discussion
on the developers@cactuscode.org mailing list, it's not practical for
me to break this up into multiple smaller single-purposed commits.]


param.ccl
* add comments
* add additional interpolator parameters (described in detail below)
* add "output psi on 2D grid" parameters (described in detail below)
* add debugging parameters

doc/documentation.tex
* add "IDAxiBrillBH/" prefix to make LaTeX labels unique across thorns
* correct a few typos
* clarify explanation of how we solve on a 2-D grid and then interpolate
  to the Cactus 3-D grid
* explain the new interpolator and output-psi-on-2D-grid parameters
* add some comments on choosing the 2-D grid resolution parameters
* mention debug parameters

doc/TODO
* new file noting some problems with this thorn which I found, but
  didn't fix

src/IDAxiBrillBH.F
* use new CCTK_WARN(CCTK_WARN_ABORT, ...) instead of old CCTK_WARN(0, ...)
* more flexible interpolator parameters:
 	This thorn first solves the Brill-wave equation on an internal
 	(axisymmetric) 2D grid, then uses uses the standard Cactus
 	CCTK_InterpLocalUniform() local-interpolator API to interpolate
 	this to the 3D grid points.  Before this patch, this thorn
 	hard-wired the interpolation operator to "uniform cartesian",
 	only allowed orders 1, 2, and 3, and didn't allow any other
 	interpolator parameters to be set.  This patch adds two new
 	parameters to allow any interpolation operator and parameters
 	to be specified.  Either the old or the new parameters can be
 	used (see doc/documentation.tex for details of how this works);
 	the defaults leave the behavior of this thorn unchanged.
* add option to output psi on 2D grid
 	When using this thorn, I found it hard to choose the parameters
 	defining the resolution and extent of the internal 2D grid.
 	To help with this, I added an option to output the solution
 	on that grid to an ASCII data file, so it can be examined
 	and plotted.
* add debugging code
 	I've added debug options which, if enabled, print the values
 	of various internal quantities at specified 2D and/or 3D grid
 	points.  I've left these in the final "production" code as
 	possibly being useful for future debugging.
* add many comments
* correct or remove some old comments which were out of date
* adjust whitespace to make the code a bit more readable:
   I've changed code like
        allocate(cc(neb,nqb),ce(neb,nqb),cw(neb,nqb),cn(neb,nqb),cs(neb,nqb),
       $     rhs(neb,nqb),psi2d(neb,nqb),detapsi2d(neb,nqb),dqpsi2d(neb,nqb),
       $     detaetapsi2d(neb,nqb),detaqpsi2d(neb,nqb),dqqpsi2d(neb,nqb),
       $     etagrd(neb),qgrd(nqb))
   to make it clearer which arrays have the same size:
        allocate(          cc(neb,nqb))
        allocate(          ce(neb,nqb))
        allocate(          cw(neb,nqb))
        allocate(          cn(neb,nqb))
        allocate(          cs(neb,nqb))
        allocate(         rhs(neb,nqb))
        allocate(       psi2d(neb,nqb))
        allocate(   detapsi2d(neb,nqb))
        allocate(     dqpsi2d(neb,nqb))
        allocate(detaetapsi2d(neb,nqb))
        allocate(  detaqpsi2d(neb,nqb))
        allocate(    dqqpsi2d(neb,nqb))

        allocate(etagrd(neb))
        allocate(  qgrd(nqb))
* change some Fortran write statements to CCTK_INFO
 	[so they're properly synchronized with output from C routines,
 	even when stdout+stderr are redirected to a log file]
* include more information about what went wrong, in various error messages
 	[eg if the interpolator returns an error code, the code
 	now includes that error code in the error message]
* rename a few variables to make their purpose clearer:
 	ep1 --> error_at_this_grid_point
 	ep2 --> max_error_in_grid


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@69 0a4070d5-58f5-498f-b6c0-2693e757fa0f

1 files changed, 101 insertions, 29 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index c5ce742..1f6c948 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -9,7 +9,7 @@
 \begin{document}
 
 \title{IDAxiBrillBH}
-\author{Paul Walker, Steve Brandt}
+\author{Paul Walker, Steve Brandt,\\some cleanups by Jonathan Thornburg}
 \date{$ $Date$ $}
 
 \maketitle
@@ -17,6 +17,8 @@
 % Do not delete next line
 % START CACTUS THORNGUIDE
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 \begin{abstract}
   Thorn IDAxiBrillBH provides analytic initial data for a vacuum
   black hole spacetime:  a single Schwarzschild black hole in
@@ -25,6 +27,8 @@
   extrinsic curvature.
 \end{abstract}
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 \section{Purpose}
 
 The pioneer, Bernstein, studied a single black hole which is
@@ -40,10 +44,10 @@ $\hat{\gamma}_{ab}$ is some chosen conformal three metric.
 The Hamiltonian constraint reduces to
 \begin{equation}
 \hat \Delta \psi = \frac{1}{8}\psi \hat R,
-\label{eqn:conformal_hamiltonian}
+\label{IDAxiBrillBH/eqn:conformal-hamiltonian}
 \end{equation}
-where $\hat \Delta$ is the covariant Laplacian and $\hat R$ is Ricci
-tensor for the conformal three metric. This form allows
+where $\hat \Delta$ is the covariant Laplacian and $\hat R$ is the
+Ricci tensor for the conformal three metric. This form allows
 us to choose an arbitrary conformal three metric, and then solve an
 elliptic equation for the conformal factor, therefore satisfying the
 constraint equations ($K_{ij} = 0$ trivially satisfies the momentum
@@ -52,13 +56,13 @@ constraints in vacuum).  This approach was used to create
 Bernstein extended this to the black hole spacetime.  Using
 spherical-polar coordinates, one can write the 3-metric,
 \begin{equation}
-\label{eqn:sph-cood}
+\label{IDAxiBrillBH/eqn:sph-coord}
 ds^2 = \psi^4 (e^{2q} (dr^2 + r^2 d \theta^2) + r^2 \sin \theta d
 \phi^2),
 \end{equation}
 where $q$ is the Brill ``packet'' which takes some functional form.
-Using this ansatz with (\ref{eqn:conformal_hamiltonian}) leads to
-an elliptic equation for $\psi$ which must be solved
+Using this ansatz with (\ref{IDAxiBrillBH/eqn:conformal-hamiltonian})
+leads to an elliptic equation for $\psi$ which must be solved
 numerically. Applying the isometry condition on $\psi$ at a finite
 radius, and applying $M/2r$ falloff conditions on $\psi$ at the
 outer boundary (the ``Robin'' condition), along with a packet which
@@ -68,15 +72,16 @@ an incident gravitational wave.  The choice of $q=0$ produces the
 Schwarzschild solution.   The typical $q$ function used in
 axisymmetry, and considered here in the non-rotating case, is 
 \begin{equation}
-q = Q_0 \sin \theta^n \left [ \exp\left(\frac{\eta -
+\label{IDAxiBrillBH/eqn:Q}
+q = Q_0 \sin^n \theta \left [ \exp\left(\frac{\eta -
       \eta_0^2}{\sigma^2}\right ) + \exp\left(\frac{\eta +
       \eta_0^2}{\sigma^2}\right ) \right ].
 \end{equation}
 Note regularity along the axis requires that the exponent $n$ must be
 even.  Choose a logarithmic radial coordinate $\eta$, which is
 related to the
-asymptoticlly flat coordinate $r$ by $\eta = ln (2r/m)$, where $m$ is a
-scale parameter.  One can rewrite (\ref{eqn:sph-cood}) as
+asymptoticlly flat coordinate $r$ by $\eta = \ln (2r/m)$, where $m$ is a
+scale parameter.  One can rewrite (\ref{IDAxiBrillBH/eqn:sph-coord}) as
 \begin{equation}
 ds^2 = \psi(\eta)^4 [ e^{2 q} (d \eta^2 +  d\theta^2) +  \sin^2
 \theta d\phi^2].
@@ -86,20 +91,20 @@ In the previous Bernstein work, the above $r$ is transformed to a
 logarithmic radial coordinate
 
 \begin{equation}
-\label{eta_coord}
+\label{IDAxiBrillBH/eta-coord}
 \eta = \ln{\frac{2r}{m}}.
 \end{equation}
 
 The scale parameter $m$ is equal to the mass of the Schwarzschild
 black hole, if $q=0$.  In this coordinate, the 3-metric is
 \begin{equation}
-\label{eqn:metric_brill_eta}
+\label{IDAxiBrillBH/eqn:metric-brill-eta}
 ds^2 = \tilde{\psi}^4 (e^{2q} (d\eta^2+d\theta^2)+\sin^2 \theta
 d\phi^2),
 \end{equation}
 and the Schwarzschild solution is 
 \begin{equation}
-\label{eqn:psi}
+\label{IDAxiBrillBH/eqn:psi}
 \tilde{\psi} = \sqrt{2M} \cosh (\frac{\eta}{2}).
 \end{equation}
 We also change the notation of $\psi$ for the conformal factor is same 
@@ -108,7 +113,7 @@ factor $r^{1/2}$ in the conformal factor.  Clearly $\psi(\eta)$ and
 $\psi$ differ by a factor of $\sqrt{r}$.  The Hamiltonian 
 constraint is
 \begin{equation}
-\label{eqn:ham}
+\label{IDAxiBrillBH/eqn:ham}
 \frac{\partial^2 \tilde{\psi}}{\partial \eta^2} + \frac{\partial^2
   \tilde{\psi}}{\partial \theta^2} + \cot \theta \frac{\partial
   \tilde{\psi}}{\partial \theta} = - \frac{1}{4} \tilde{\psi}
@@ -122,14 +127,14 @@ we substitute
 \delta \tilde{\psi} & = & \tilde{\psi}+\tilde{\psi}_0 \\
                     & = & \tilde{\psi}-\sqrt{2m} \cosh(\frac{\eta}{2}).
 \end{eqnarray}
-to the equation~(\ref{eqn:ham}), then we can linearize it as
+to the equation~(\ref{IDAxiBrillBH/eqn:ham}), then we can linearize it as
 \begin{equation}
 \frac{\partial^2 \delta\tilde{\psi}}{\partial \eta^2} + \frac{\partial^2
   \delta\tilde{\psi}}{\partial \theta^2} + \cot \theta \frac{\partial
   \delta\tilde{\psi}}{\partial \theta} = - \frac{1}{4}
 (\delta\tilde{\psi} + \tilde{\psi}_0) (\frac{\partial^2 q}{\partial
   \eta^2} + \frac{\partial^2 q}{\partial \theta^2} -1).
-\label{eqn:ham_linear}
+\label{IDAxiBrillBH/eqn:ham-linear}
 \end{equation}
 For the boundary conditions, we use for the inner boundary condition
 an isometry condition:
@@ -142,19 +147,84 @@ and outer boundary condition, a Robin condition:
 \tilde{\psi})|_{\eta=\eta_{max}} = 0.
 \end{equation}
 
-%  [[ DPR: What is this: ?? ]]
-%This thorn provides
-% \begin{enumerate}
-%  \item CactusEinstein
-% \end{enumerate}
-
-\section{Comments}
-
-We calculate equation~(\ref{eqn:ham_linear}) with spherical
-coordinates.  However, Cactus needs Cartesian coordinates.  Therefore,
-we interpolate $\psi$ to the Cartesian grid by using an interpolator.
-Note that the interpolator has linear, quadratic, and cubic
-interpolation.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{2-D Grid and Interpolation Parameters}
+
+This thorn solves equation~(\ref{IDAxiBrillBH/eqn:ham-linear}) on a 2-D
+$(\eta,\theta)$ grid.  However, Cactus needs a 3-D grid, typically with
+Cartesian coordinates.  Therefore, this thorn interpolate $\psi$ and its
+$(\eta,\theta)$ derivatives to the Cartesian grid.
+
+The parameters \verb|neta| and \verb|nq| specify the resolution of
+this thorn's 2-D grid in $\eta$ and $\theta$ respectively.%%%
+\footnote{%%%
+	 Internally, this thorn uses ``$q$'' to refer
+	 to $\theta$ in Fortran code, with the $q$ function
+	 of~$(\protect\ref{IDAxiBrillBH/eqn:Q})$ being hidden
+	 in the Mathematica files (and not present in the Fortran
+	 code).  Noone seems to know \emph{why} the code does
+	 things this way\dots{}  Unfortunately, this renaming
+	 has leaked out into the parameter names\dots
+	 }%%%
+{}  The default values are a reasonable starting point, but you may
+need to increase them substantially if you need very high accuracy
+(very small constraint violations).
+
+To help judge what resolution may be needed, this thorn has an option
+to write out $\psi(\eta)$ and $\psi$ on the 2-D grid to an ASCII data file
+where it can be examined and/or plotted.  To do this, set the Boolean
+parameter \verb|output_psi2D|, and possibly also the string parameter
+\verb|output_psi2D_file_name|.
+
+This thorn uses the standard Cactus \verb|CCTK_InterpLocalUniform()|
+local interpolation system for this interpolation.  The interpolation
+operator is specified with the \verb|interpolator_name| parameter
+(this defaults to \verb|"uniform cartesian"|, the interpolation
+operator provided by thorn \textbf{CactusBase/LocalInterp}).
+
+The interpolation order and/or other parameters can be specified
+in either of two ways:%%%
+\footnote{%%%
+	 Notice that, for historical reasons, the
+	 interpolation parameter names are a bit
+	 inconsistent: \texttt{interpolat\underline{ion}\_order}
+	 versus \texttt{interpolat\underline{or}\_name}
+	 and \texttt{interpolat\underline{or}\_pars}.
+	 }%%%
+\begin{itemize}
+\item	The integer parameter \verb|interpolation_order| may be
+	used directly to specify the interpolation order.
+\item	More generally, the string parameter \verb|interpolator_pars|
+	may be set to any nonempty string (it defaults to the empty string).
+	If this is done, this parameter overrides \verb|interpolation_order|,
+	and explicitly specifies a parameter string for the interpolator.
+\end{itemize}
+Note that the default interpolator parameters specify \emph{linear}
+interpolation.  This is rather inaccurate, and (due to aliasing effects
+between the 2-D and 3-D grids) will give a fair bit of noise in the
+metric components.  You may want to specify a higher-order interpolator
+to reduce this noise.
+
+For example, for one test series where I (JT) needed very accurate
+initial data (I wanted the initial-data errors to be much less than
+the errors from 4th~order finite differencing on the 3-D Cactus grid),
+I had to use a resolution of $1000$ in $\eta$ and $2000$ in $\theta$,
+together with either 4th~order Lagrange or 3rd~order Hermite interpolation
+(provided by thorn \textbf{AEIThorns/AEILocalInterp}) to get sufficient
+accuracy.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Debugging Parameters}
+
+This thorn has options to print very detailed debugging information
+about internal quantities at selected grid points.  This is enabled
+by setting the integer parameter \verb|debug| to a positive value
+(the default is $0$, which means no debugging output).  See
+\verb|param.ccl| and the source code \verb|src/IDAxiBrillBH.F| for details.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \bibliographystyle{prsty}
 \begin{thebibliography}{10}
@@ -167,6 +237,8 @@ interpolation.
   K. Camarda, Ph.D thesis University of Illinois Urbana-Champaign, (1998)
 \end{thebibliography}
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 % Do not delete next line
 % END CACTUS THORNGUIDE