Homogenise parameters, changed documentation

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@82 a678b1cf-93e1-4b43-a69d-d43939e66649

9 files changed, 131 insertions, 133 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index a17be7a..2703ebb 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -13,9 +13,9 @@ cartesian grid), as well as data with an angular dependency.}
 
 \section{Purpose}
 
-The purpose of this thorn is to create initial data for a Brill wave
-spacetime.  It does so by starting from a three--metric of the form
-originally considered by Brill
+The purpose of this thorn is to create (time symmetric) initial data
+for a Brill wave spacetime.  It does so by starting from a
+three--metric of the form originally considered by Brill
 \begin{equation}
 ds^2 = \Psi^4 \left[ e^{2q} \left( d\rho^2 + dz^2 \right)
 + \rho^2 d\phi^2 \right] =\Psi^4 \hat{ds}^{2},
@@ -25,7 +25,29 @@ where $q$ is a free function subject to certain regularity and
 fall-off conditions, $\rho=\sqrt{x^2+y^2}$ and $\Psi$ is a conformal
 factor to be solved for.
 
-Substituting this metric into the Hamiltonian constraint gives an
+Thorn {\tt IDBrillData} provides three choices for the $q$ function: 
+an exponential form, ({\tt IDBrillData::q\_function = "exp"})
+\begin{equation}
+q = a \; \frac{\rho^{2+b}}{r^2} \left( \frac{z}{\sigma_z} \right)^2
+e^{-(\rho - \rho_0^2)} \left[ 1 + d \frac{\rho^m}{1 + e \rho^m}
+\cos^2 \left( n \phi + \phi_0 \right) \right]
+\end{equation}
+a generalized form of the $q$ function first written down by Eppley
+({\tt IDBrillData::q\_function = "eppley"}) 
+\begin{equation}
+q = a \left( \frac{\rho}{\sigma_\rho} \right)^b \frac{1}{1 + \left[
+\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}}\left[ 1 + d \frac{\rho^m}{1 + e \rho^m}
+\cos^2 \left( n \phi + \phi_0 \right) \right]
+\end{equation}
+and the (default) Gundlach $q$ function which includes the Holz form
+({\tt IDBrillData::q\_function = "gundlach"})
+\begin{equation} 
+q = a \left( \frac{\rho}{\sigma_\rho} \right)^b e^{-\left[
+\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}} \left[ 1 + d \frac{\rho^m}{1 + e \rho^m}
+\cos^2 \left( n \phi + \phi_0 \right) \right]
+\end{equation}
+
+Substituting the metric into the Hamiltonian constraint gives an
 elliptic equation for the conformal factor $\Psi$ which is then
 numerically solved for a given function $q$:
 \begin{equation}
@@ -39,14 +61,9 @@ where the conformal Ricci scalar is found to be
 Assuming the initial data to be time symmetric means that the momentum
 constraints are trivially satisfied.
 
-The thorn considers several different forms of the function $q$
-depending on certain parameters that will be described below.
-
-Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase}
-parameter {\tt initial\_data} to be {\tt brilldata}.
-
-In the case of axisymmetry, the Hamiltonian constraint can be written
-as an elliptic equation for $\Psi$ with just the flat space Laplacian,
+In the case of axisymmetry (that is $d=0$ in the above expressions for
+$q$), the Hamiltonian constraint can be written as an elliptic
+equation for $\Psi$ with just the flat space Laplacian,
 \begin{equation}
 \nabla_{flat} \Psi + \frac{\Psi}{4} (\partial_z^2 q + \partial_\rho^2 q) = 0
 \end{equation}
@@ -55,82 +72,60 @@ ADMBase::initial\_data = "brilldata2D"} then this elliptic equation
 is solved rather than the equation above.
 
 
-\section{Parameters for the thorn}
+\section{Generating Initial Data with IDBrillData}
 
-The thorn is controlled by the following parameters:
+Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase}
+parameter {\tt initial\_data} to be {\tt brilldata}, or for the case of
+axisymmetry {\tt brilldata2D} can also be used.
 
-\begin{itemize}
+The parameter {\tt IDBrillData::q\_function} chooses the form of the
+$q$ function to be used, defaulting to the Gundlach expression.
 
-\item brill\_q (INT):  Form of the function $q$ [0,1,2] (default 2):
+Additional {\tt IDBrillData} parameters for each form of $q$ fix the 
+remaining freedom:
 
 \begin{itemize}
 
-\item brill\_q = 0:
-\[
-q = a \; \frac{\rho^{2+b}}{r^2} \left( \frac{z}{\sigma_z} \right)^2
-e^{-(\rho - \rho_0^2)}
-\]
+\item Exponential $q$: {\tt IDBrillData::q\_function = "exp"}
 
-\item  brill\_q = 1:
-\[
-q = a \left( \frac{\rho}{\sigma_\rho} \right)^b \frac{1}{1 + \left[
-\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}}
-\]
+$(a, b,\sigma_z,\rho_0)=$ ({\tt exp\_a, exp\_b, exp\_sigmaz,exp\_rho0})
 
-\item  brill\_q = 2:
-\[
-q = a \left( \frac{\rho}{\sigma_\rho} \right)^b e^{-\left[
-\left( r^2 - r_0^2 \right) / \sigma_r^2 \right]^{c/2}}
-\]
+\item Eppley $q$: {\tt IDBrillData::q\_function = "eppley"}
 
-\item If one specifies 3D data (see above), the function $q$ is multiplied
-by an additional factor with an angular dependency:
-\[
-q \rightarrow q \left[ 1 + d \frac{\rho^m}{1 + e \rho^m}
-\cos^2 \left( n \phi + \phi_0 \right) \right]
-\]
+$(a, b,\sigma_\rho, r_0,\sigma_r,c)=$ ({\tt eppley\_a, eppley\_b, eppley\_sigmarho, eppley\_r0, eppley\_sigmar, eppley\_c})
 
-\end{itemize}
-
-\item brill\_a (REAL): Amplitude (default 0.0).
-
-\item brill\_b (REAL):  $b$ in above expressions (default 2.0).
-
-\item brill\_c (REAL):  $c$ in above expressions (default 2.0).
-
-\item brill\_d (REAL):  $d$ in above expressions (default 0.0).
+\item Gundlach $q$: {\tt IDBrillData::q\_function = "gundlach"}
 
-\item brill\_e (REAL):  $e$ in above expressions (default 1.0).
+$(a, b,\sigma_\rho, r_0,\sigma_r,c)=$ ({\tt gundlach\_a, gundlach\_b, gundlach\_sigmarho, gundlach\_r0, gundlach\_sigmar, gundlach\_c})
 
-\item brill\_m (REAL):  $m$ in above expressions (default 2.0).
+\item Non-axisymmetric part for each choice of $q$
 
-\item brill\_n (REAL):  $n$ in above expressions (default 2.0).
+$(d, m, e, n, \phi0)=$ ({\tt brill3d\_d, brill3d\_m, brill3d\_e, brill3d\_n, brill3d\_phi0})
 
-\item brill\_r0 (REAL):  $r_0$ in above expressions (default 0.0).
-  
-\item brill\_rho0 (REAL):  $\rho_0$ in above expressions
-  (default 0.0).
-  
-\item brill\_phi0 (REAL):  $\phi_0$ in above expressions
-  (default 0.0).
-  
-\item brill\_sr (REAL):  $\sigma_r$ in above expressions
-  (default 1.0).
-  
-\item brill\_srho (REAL):  $\sigma_\rho$ in above
-  expressions (default 1.0).
- 
 \end{itemize}
 
-The elliptic solver is controlled by the additional parameters:
+Note that the default $q$ expression is
+$$
+q = {\tt gundlach\_a} \quad \rho^2 e^{-r^2} 
+$$
+
+{\tt IDBrillData} can use the elliptic solvers (type LinMetric)
+provided by {\tt CactusEinstein/EllSOR}, {\tt AEIThorns/BAM\_Elliptic}
+or {\tt CactusElliptic/EllPETSc} to solve the equation resulting from 
+the Hamiltonian constraint. 
+In all cases the parameter {\tt thresh} sets the threshold for the elliptic
+solve. The choice of elliptic solver is made 
+through the parameter {\tt brill\_solver}:
 
 \begin{itemize}
   
-\item {\tt solver} (KEYWORD): Elliptic solver used to solve the
-  hamiltonian constraint [sor/petsc/bam] (default "sor").
+\item {\tt sor}: Understands the Robin boundary condition, additional 
+parameters control the maximum number of iterations ({\tt sor\_maxit}).
   
-\item {\tt thresh} (REAL): Threshold for elliptic solver (default
-  0.00001).
+\item {\tt bam}: {\tt BAM\_Elliptic} does not properly implement the
+elliptic infrastructure of {\tt EllBase}, and the {\tt BAM\_Elliptic}
+parameter to use the Robin boundary condition must be set independently
+of {tt IDBrillWave::brill\_bound}.
 
 \end{itemize}
 
@@ -143,6 +138,29 @@ the derivatives of the conformal factor are not calculated, so that
 only {\tt staticconformal::conformal\_storage = "factor"} is
 supported.
 
+\section{References}
+
+\subsection{Specification of Brill Waves}
+
+\begin{enumerate}
+
+\item Dieter Brill, {\bf Ann. Phys.}, 7, 466, 1959.
+
+\item Ken Eppley, {\bf Sources of Gravitational Radiation}, edited by L. Smarr (Cambridge University Press, 
+Cambridge, England, 1979), p. 275.
+
+\end{enumerate}
+
+\subsection{Numerical Evolutions of Brill Waves}
+
+\begin{enumerate}
+
+\item {\it Gravitational Collapse of Gravitational Waves in 3D Numerical Relativity}, 
+ Miguel Alcubierre, Gabrielle Allen, Bernd Bruegmann, Gerd Lanfermann, Edward Seidel, Wai-Mo Suen, Malcolm Tobias,
+{\bf Phys. Rev. D61}, 041501, 2000.
+
+\end{enumerate}
+
 % Automatically created from the ccl files 
 % Do not worry for now.
 
diff --git a/param.ccl b/param.ccl
index 1313ac5..4586bab 100644
--- a/param.ccl
+++ b/param.ccl
@@ -20,35 +20,35 @@ private:
 
 # Parameters for elliptic solve
 
-KEYWORD brill_solver "Which elliptic solver to use"
+KEYWORD solver "Which elliptic solver to use"
 {
   "sor"    ::  "Use SOR solver"
   "petsc"  ::  "Use PETSc solver"
   "bam"    ::  "Use bam solver"
 } "sor"
 
-KEYWORD brill_bound "Which boundary condition to use"
+KEYWORD bound "Which boundary condition to use"
 {
-  "const" :: "constant boundary: set brill_const_v0"
-  "robin" :: "Robin boundary: set brill_robin_falloff, brill_robin_inf"
-} "const"
+  "const" :: "constant boundary: set const_v0"
+  "robin" :: "Robin boundary: set robin_falloff, robin_inf"
+} "robin"
 
-INT brill_robin_falloff "Fall-off of Robin BC"
+INT robin_falloff "Fall-off of Robin BC"
 {
   0: :: "any positive integer value"
 } 1
 
-REAL brill_const_v0 "Value of constant BC"
+REAL const_v0 "Value of constant BC"
 {
   : :: "anything goes"
 } 1.0
 
-REAL brill_robin_inf "Value at infinity of Robin BC"
+REAL robin_inf "Value at infinity of Robin BC"
 {
   : :: "anything goes"
 } 1.0
 
-REAL brill_thresh "How far (absolute norm) to go"
+REAL thresh "How far (absolute norm) to go"
 {
   0.0: :: "Positive number please"
 } 0.00001
@@ -181,7 +181,7 @@ REAL brill3d_phi0 "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 
 # Additional parameters
 
-REAL brill_rhofudge "delta rho for axis fudge"
+REAL rhofudge "delta rho for axis fudge"
 {
   0: :: "Positive please"
 } 0.00001
@@ -189,7 +189,7 @@ REAL brill_rhofudge "delta rho for axis fudge"
 INT sor_maxit "Maximum number of iterations"
 {
 0:* :: "Positive"
-} 100
+} 10000
 
 BOOLEAN output_coeffs "output coefficients for elliptic solve"
 {
diff --git a/src/brilldata.F b/src/brilldata.F
index 3c91b97..e38afaa 100644
--- a/src/brilldata.F
+++ b/src/brilldata.F
@@ -75,7 +75,7 @@ c     Set up background metric and coefficients for linear solve.
 
 c     Tolerances for elliptic solve.
 
-      AbsTol(1)= brill_thresh
+      AbsTol(1)= thresh
       AbsTol(2)= -1.0D0
       AbsTol(3)= -1.0D0
  
@@ -94,14 +94,14 @@ c     Elliptic solver
 c     Just in case some solvers do not use the standard interface
       conformal_state = 0
 
-      if (CCTK_EQUALS(brill_solver,"sor")) then
+      if (CCTK_EQUALS(solver,"sor")) then
          call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit")
          call Ell_LinMetricSolver(ierr,cctkGH, 
      .   metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor")
-      else if (CCTK_EQUALS(brill_solver,"petsc")) then
+      else if (CCTK_EQUALS(solver,"petsc")) then
          call Ell_LinMetricSolver(ierr,cctkGH, 
      .   metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"petsc")
-      else if (CCTK_EQUALS(brill_solver,"bam")) then
+      else if (CCTK_EQUALS(solver,"bam")) then
          call Ell_LinMetricSolver(ierr,cctkGH, 
      .   metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"bam")
       end if
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index 467c28d..bd13c79 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -21,10 +21,9 @@
       DECLARE_CCTK_FUNCTIONS
 
       integer i,j,k
-      integer nx,ny,nz
 
       CCTK_REAL x1,y1,z1,rho1,rho2
-      CCTK_REAL phi,psi4,e2q,rhofudge
+      CCTK_REAL phi,psi4,e2q
       CCTK_REAL zero,one
 
       CCTK_REAL brillq,phif
@@ -34,22 +33,12 @@ c     Numbers.
       zero = 0.0D0
       one  = 1.0D0
 
-c     Set up grid size.
-
-      nx = cctk_lsh(1)
-      ny = cctk_lsh(2)
-      nz = cctk_lsh(3)
-
-c     Parameters.
-
-      rhofudge = brill_rhofudge
-      
 c     Replace flat metric left over from elliptic solve by
 c     Brill metric calculated from q and Psi.
 
-      do k=1,nz
-         do j=1,ny
-            do i=1,nx
+      do k=1,cctk_lsh(3)
+         do j=1,cctk_lsh(2)
+            do i=1,cctk_lsh(1)
 
                x1 = x(i,j,k)
                y1 = y(i,j,k)
@@ -84,13 +73,23 @@ c              The individual coefficients can be read off as
 c                 This fudge assumes that q = O(rho^2) near the axis. Which
 c                 it should be, or the data will be singular.
 
-                  gxx(i,j,k) = 0.0d0
-                  gyy(i,j,k) = 0.0d0
-                  gzz(i,j,k) = 0.0d0
-                  gxy(i,j,k) = 0.0d0
+                  gxx(i,j,k) = zero
+                  gyy(i,j,k) = zero
+                  gzz(i,j,k) = zero
+                  gxy(i,j,k) = zero
 
                end if
 
+               gxz(i,j,k) = zero
+               gyz(i,j,k) = zero
+
+               kxx(i,j,k) = zero
+               kyy(i,j,k) = zero
+               kzz(i,j,k) = zero
+               kxy(i,j,k) = zero
+               kxz(i,j,k) = zero
+               kyz(i,j,k) = zero
+
             end do
          end do
       end do
@@ -99,9 +98,9 @@ c                 it should be, or the data will be singular.
 
          conformal_state = 1
 
-         do k=1,nz
-            do j=1,ny
-               do i=1,nx
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
 
                   psi(i,j,k) = brillpsi(i,j,k)
                   
@@ -113,9 +112,9 @@ c                 it should be, or the data will be singular.
 
          conformal_state = 0
 
-         do k=1,nz
-            do j=1,ny
-               do i=1,nx
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
 
                   psi4  = brillpsi(i,j,k)**4
 
@@ -130,19 +129,5 @@ c                 it should be, or the data will be singular.
 
       end if
 
-c     In any case,
-
-      gxz = zero
-      gyz = zero
-
-c     Vanishing extrinsic curvature completes the Cauchy data.
-
-      kxx = zero
-      kyy = zero
-      kzz = zero
-      kxy = zero
-      kxz = zero
-      kyz = zero
-
       return
       end
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
index f2e5f13..185bef5 100644
--- a/src/setupbrilldata2D.F
+++ b/src/setupbrilldata2D.F
@@ -103,7 +103,7 @@ c              will not be regular on the axis.
 c     Synchronization and boundaries.
 
 c      call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
-      call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
+c      call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
 
       return
       end
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
index 97cf28b..b77d886 100644
--- a/src/setupbrilldata3D.F
+++ b/src/setupbrilldata3D.F
@@ -47,7 +47,7 @@ c      c        \         z      rho                    phi          phi    /
 
       CCTK_REAL x1,y1,z1,rho1,rho2
       CCTK_REAL phi,phif,e2q
-      CCTK_REAL brillq,rhofudge,eps
+      CCTK_REAL brillq,eps
       CCTK_REAL zero,one,two,three
 
 c     Define numbers
@@ -57,10 +57,6 @@ c     Define numbers
       two   = 2.0D0
       three = 3.0D0
 
-c     Parameters.
-
-      rhofudge = brill_rhofudge
-
 c     Epsilon for finite differencing.
 
       eps = cctk_delta_space(1)
diff --git a/test/test_brilldata_1.par b/test/test_brilldata_1.par
index 5d673ba..0d15a58 100644
--- a/test/test_brilldata_1.par
+++ b/test/test_brilldata_1.par
@@ -27,7 +27,7 @@ idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
-idbrilldata::brill_solver 			= "bam"
+idbrilldata::solver 				= "bam"
 bam_elliptic::bam_bound                 	= "bamrobin"
 
 # Output.
diff --git a/test/test_brilldata_2.par b/test/test_brilldata_2.par
index 630cd64..4e0b101 100644
--- a/test/test_brilldata_2.par
+++ b/test/test_brilldata_2.par
@@ -28,8 +28,7 @@ idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
-idbrilldata::brill_solver 			= "bam"
-
+idbrilldata::solver 			 	= "bam"
 bam_elliptic::bam_bound                 	= "bamrobin"
 
 # Output.
diff --git a/test/test_wsor.par b/test/test_wsor.par
index a9b0337..4c04382 100644
--- a/test/test_wsor.par
+++ b/test/test_wsor.par
@@ -26,8 +26,8 @@ idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
-idbrilldata::brill_solver 			= "sor"
-idbrilldata::brill_bound               		= "robin"
+idbrilldata::solver 				= "sor"
+idbrilldata::bound               		= "robin"
 idbrilldata::sor_maxit				= 100
 # Warning -- since this doesn't converge for so few iterations, the
 # result will depend on number of cpus, even beyond the sor tolerance