Lots of tidying changes, parameter names have changed but hopefully clearer now.

Few more documentation changes still to come

Constructing 3D data doesn't work, but this was also true before Einstein2 changes.


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

13 files changed, 269 insertions, 240 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 49731f9..a17be7a 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -3,13 +3,13 @@
 \begin{document}
 
 \title{IDBrillData}
-\author{Carsten Gundlach}
-\date{6 September 1999}
+\author{Carsten Gundlach, Gabrielle Allen}
+\date{$Date$}
 \maketitle
 
-\abstract{This thorn creates initial data for Brill wave spacetimes.
-It can create both axisymmetric data (in a 3D cartesian grid), as
-well as data with an angular dependency.}
+\abstract{This thorn creates time symmetric initial data for Brill
+wave spacetimes.  It can create both axisymmetric data (in a 3D
+cartesian grid), as well as data with an angular dependency.}
 
 \section{Purpose}
 
@@ -22,40 +22,38 @@ ds^2 = \Psi^4 \left[ e^{2q} \left( d\rho^2 + dz^2 \right)
 \label{eqn:brillmetric}
 \end{equation}
 where $q$ is a free function subject to certain regularity and
-fall-off conditions and $\Psi$ is a conformal factor to be solved for.
+fall-off conditions, $\rho=\sqrt{x^2+y^2}$ and $\Psi$ is a conformal
+factor to be solved for.
 
-The thorn considers several different forms of the function $q$
-depending on certain parameters that will be described below.
-Substituting the metric above into the Hamiltonian constraint results
-in an elliptic equation for the conformal factor $\Psi$ that can be
-solved numerically once the function $q$ has been specified. The
-initial data is also assumed to be time-symmetric, so the momentum
+Substituting this 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}
+\hat{\nabla} \Psi - \frac{\Psi}{8} \hat{R} = 0
+\end{equation}
+where the conformal Ricci scalar is found to be
+\begin{eqnarray}
+\hat{R} = -2 \left(e^{-2q} (\partial^2_z q + \partial^2_\rho q) + 
+\frac{1}{\rho^2} (3 (\partial_\phi q)^2 + 2 \partial_\phi q)\right)
+\end{eqnarray}
+Assuming the initial data to be time symmetric means that the momentum
 constraints are trivially satisfied.
 
-% [[ DPR: The code does not use these two different initial_data
-% keywords, but the axisym parameter instead.  It should probably use
-% the two different initial_data values.  I have changed the doc below
-% to correspond with the current code: ]]
+The thorn considers several different forms of the function $q$
+depending on certain parameters that will be described below.
 
-The thorn is activated by choosing the CactusEinstein/ADMBase parameter
-``initial\_data'' to be:
-%in one of the following two ways:
+Brill initial data is activated by choosing the {\tt CactusEinstein/ADMBase}
+parameter {\tt initial\_data} to be {\tt brilldata}.
 
-\begin{itemize}
-  
-\item initial\_data = ``brilldata'': %Axisymmetric 
-Brill wave initial data
-%  (but calculated in a cartesian grid!).
-  
-%\item initial\_data = ``brilldata3D'': Brill wave initial data with an
-%  angular dependency.
-
-\end{itemize}
+In the case of axisymmetry, 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}
+If the initial data is chosen to be {\tt
+ADMBase::initial\_data = "brilldata2D"} then this elliptic equation
+is solved rather than the equation above.
 
-To choose axisymmetric data, set the parameter {\tt axisym} to true
-(the default).  Note that the data is computed on a cartesian grid in
-any event.  {\tt axisym = "no"} adds an angular dependence factor,
-which is detailed below.
 
 \section{Parameters for the thorn}
 
@@ -128,14 +126,23 @@ The elliptic solver is controlled by the additional parameters:
 
 \begin{itemize}
   
-\item solver (KEYWORD): Elliptic solver used to solve the
-  hamiltonian constraint [sor/petsc/bam] (default ``sor'').
+\item {\tt solver} (KEYWORD): Elliptic solver used to solve the
+  hamiltonian constraint [sor/petsc/bam] (default "sor").
   
-\item thresh (REAL): Threshold for elliptic solver (default
+\item {\tt thresh} (REAL): Threshold for elliptic solver (default
   0.00001).
 
 \end{itemize}
 
+\section{Notes}
+
+Thorn {\tt IDBrillData} understands both the {\tt physical} and {\tt
+static conformal} {\tt metric\_type}. In the case of a conformal
+metric being chosen, the conformal factor is set to $\Psi$. Currently
+the derivatives of the conformal factor are not calculated, so that
+only {\tt staticconformal::conformal\_storage = "factor"} is
+supported.
+
 % Automatically created from the ccl files 
 % Do not worry for now.
 
diff --git a/param.ccl b/param.ccl
index d657f32..1a37fe4 100644
--- a/param.ccl
+++ b/param.ccl
@@ -6,6 +6,7 @@ shares: ADMBase
 EXTENDS KEYWORD initial_data 
 {
    "brilldata"   :: "Brill wave initial data"
+   "brilldata2D" :: "Brill wave initial data assuming axisymmetry"
 }
 
 USES KEYWORD metric_type
@@ -55,78 +56,124 @@ REAL brill_thresh "How far (absolute norm) to go"
 
 # Brill wave parameters
 
-BOOLEAN axisym "Axisymmetric initial data?"
+KEYWORD q_function  "Form of function q [0,1,2]"
 {
-} "yes"
+  "exp" :: ""
+  "eppley" ::
+  "gundlach" ::
+} "gundlach"
 
-INT brill_q  "Form of function q [0,1,2]"
+# -------------------------------------------------------------
+
+REAL exp_a "Exp Brill wave: Amplitude"
+{
+  : :: "Anything"
+} 0.0
+
+INT exp_b "Exp Brill wave: used in exponent in rho: rho^(2+b)"
 {
-  0:2 :: "Only cases 0,1,2 defined at the moment"
+  : :: "Anything"
 } 2
 
-REAL brill_a "Brill wave: Amplitude"
+REAL exp_rho0 "Exp Brill wave: radius of torus in rho"
+{
+  0:* :: "Positive"
+} 0.0
+
+REAL exp_sigmaz "Exp Brill wave: sigma in z"
+{
+  (0:* :: "Positive"
+} 1.0
+
+# -------------------------------------------------------------
+
+REAL eppley_a "Eppley Brill wave: Amplitude"
 {
   : :: "Anything"
 } 0.0
 
-INT brill_b "Brill wave: used in exponent in rho"
+INT eppley_b "Eppley Brill wave: used in exponent in rho: rho^b"
 {
   : :: "Anything"
 } 2
 
-INT brill_c "Brill wave: (r^2 - r0^2)^(c/2)"
+INT eppley_c "Eppley Brill wave: (r^2 - r0^2)^(c/2)"
 {  
   : :: "Anything"
 } 2
 
-REAL brill_rho0 "Brill wave: radius of torus in rho"
+REAL eppley_r0 "Eppley Brill wave: radius of torus in r"
 {
-  : :: "Anything"
+  (0:* :: "Positive" 
 } 0.0
 
-REAL brill_r0 "Brill wave: radius of torus in r"
+REAL eppley_sigmarho "Eppley Brill wave: sigma in rho"
 {
-  : :: "Anything" 
-} 0.0
+  : :: "Anything"
+} 1.0
 
-REAL brill_srho "Brill wave: sigma in rho"
+REAL eppley_sigmar "Eppley Brill wave: sigma in r"
 {
   : :: "Anything"
 } 1.0
 
-REAL brill_sr "Brill wave: sigma in r"
+# -------------------------------------------------------------
+
+REAL gundlach_a "Gundlach Brill wave: Amplitude"
+{
+  : :: "Anything"
+} 0.0
+
+INT gundlach_b "Gundlach Brill wave: used in exponent in rho: rho^b"
+{
+  : :: "Anything"
+} 2
+
+REAL gundlach_sigmarho "Gundlach Brill wave: sigma in rho"
 {
   : :: "Anything"
 } 1.0
 
-REAL brill_sz "Brill wave: sigma in z"
+REAL gundlach_r0 "Gundlach Brill wave: radius of torus in r"
+{
+  (0:* :: "Positive" 
+} 0.0
+
+INT gundlach_c "Gundlach Brill wave: (r^2 - r0^2)^(c/2)"
+{  
+  : :: "Anything"
+} 2
+
+REAL gundlach_sigmar "Gundlach Brill wave: sigma in r"
 {
   : :: "Anything"
 } 1.0
 
+# -------------------------------------------------------------
+
 # 3D Brill wave parameters
 
-REAL brill_d "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+REAL brill3d_d "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 {
   : :: "Anything"
 } 0.0
 
-REAL brill_e "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+REAL brill3d_e "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 {
   : :: "Anything"
 } 1.0
 
-REAL brill_m "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+REAL brill3d_m "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 {
   : :: "Anything"
 } 2.0
 
-REAL brill_n "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+REAL brill3d_n "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 {
   : :: "Anything"
 } 2.0
 
-REAL brill_phi0 "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+REAL brill3d_phi0 "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
 {
   : :: "Anything"
 } 0.0
@@ -143,3 +190,7 @@ INT sor_maxit "Maximum number of iterations"
 {
 0:* :: "Positive"
 } 100
+
+BOOLEAN output_coeffs "output coefficients for elliptic solve"
+{
+} "no"
\ No newline at end of file
diff --git a/schedule.ccl b/schedule.ccl
index 4d06086..0e8035e 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,7 +1,7 @@
 # Schedule definitions for thorn IDBrillData
 # $Header$
 
-if (CCTK_Equals(initial_data,"brilldata"))
+if (CCTK_Equals(initial_data,"brilldata") || CCTK_Equals(initial_data,"brilldata2D"))
 {
 
   STORAGE: brillconf
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index fe55173..7a04cec 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -52,7 +52,9 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS)
   if( ! CCTK_Equals(metric_type, "physical") && 
       ! CCTK_Equals(metric_type, "static conformal"))
   {
-    CCTK_PARAMWARN("Unknown ADMBase::metric_type - known types are \"physical\" and \"static conformal\"");
+    CCTK_VParamWarn("Unknown ADMBase::metric_type (%s): "
+		    "known types are \"physical\" and \"static conformal\"",
+		    metric_type);
   }
 
   if (CCTK_Equals(metric_type, "static conformal"))
@@ -72,6 +74,33 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS)
   {
     CCTK_INFO("  ... using physical metric");
   }
+
+  if (CCTK_Equals(initial_data,"brilldata"))
+  {
+    CCTK_INFO("  ... using full 3D solver");
+  }
+  else if (CCTK_Equals(initial_data,"brilldata2D"))
+  {
+    CCTK_INFO("  ... using reduced 2D solver");
+    if (brill3d_d != 0.0)
+    {
+      CCTK_WARN(0,"But 3D parameter brill3d_d set for the 2D solver !?!");
+    }
+  }
+
+  if (CCTK_Equals(q_function,"exp"))
+  {
+    CCTK_INFO("  ... Exponential brill function");
+  }
+  else if (CCTK_Equals(q_function,"eppley"))
+  {
+    CCTK_INFO("  ... generalised Eppley Brill function");
+  }
+  else if (CCTK_Equals(q_function,"exp"))
+  {
+    CCTK_INFO("  ... Gundlach-Holz Brill function");
+  }
+
 }
 
 /********************************************************************
diff --git a/src/Startup.c b/src/Startup.c
index 8660cbe..a6d4c13 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -18,7 +18,9 @@ void BrillData_InitSymBound(CCTK_ARGUMENTS)
   sym[1] = 1;
   sym[2] = 1;
   
-  SetCartSymVN(cctkGH, sym,"IDBrillData::brillpsi");
+  SetCartSymVN(cctkGH,sym,"IDBrillData::brillpsi");
+  SetCartSymVN(cctkGH,sym,"IDBrillData::brillMlinear");
+  SetCartSymVN(cctkGH,sym,"IDBrillData::brillNsource");
 
   return;
 } 
diff --git a/src/brilldata.F b/src/brilldata.F
index 7dc01e0..3c91b97 100644
--- a/src/brilldata.F
+++ b/src/brilldata.F
@@ -46,22 +46,33 @@ c     Get indices for grid functions.
 
       call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear")
       if (iMcoeff.lt.0) then
-	call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
+	call CCTK_WARN(0,"Grid variable index for brillMlinear not found")
       end if
 
       call CCTK_VarIndex(iNcoeff,"idbrilldata::brillNsource")
       if (iNcoeff.lt.0) then
-	call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
+	call CCTK_WARN(0,"Grid variable index for brillNsource not found")
       end if
 
 c     Set up background metric and coefficients for linear solve.
 
-      if (axisym.eq.1) then
+      if (CCTK_EQUALS(initial_data,"brilldata2D")) then
          call setupbrilldata2D(CCTK_ARGUMENTS)
       else
          call setupbrilldata3D(CCTK_ARGUMENTS)
       end if
 
+      if (output_coeffs .eq. 1) then
+         call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillMlinear","Debug_brillMlinear")
+         call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillNsource","Debug_brillNsource")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxx","Debug_gxx")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxy","Debug_gxy")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxz","Debug_gxz")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyy","Debug_gyy")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyz","Debug_gyz")
+         call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gzz","Debug_gzz")
+      end if
+
 c     Tolerances for elliptic solve.
 
       AbsTol(1)= brill_thresh
@@ -74,14 +85,17 @@ c     Tolerances for elliptic solve.
 
 c     Boundaries.
 
-      call Ell_SetStrKey(ierr,"yes","EllLinConfMetric::Bnd::Robin")
-      call Ell_SetRealKey(ierr,1.0D0,"EllLinConfMetric::Bnd::Robin::inf")
-      call Ell_SetIntKey(ierr,1,"EllLinConfMetric::Bnd::Robin::falloff")
+      call Ell_SetStrKey(ierr,"yes","EllLinMetric::Bnd::Robin")
+      call Ell_SetRealKey(ierr,1.0D0,"EllLinMetric::Bnd::Robin::inf")
+      call Ell_SetIntKey(ierr,1,"EllLinMetric::Bnd::Robin::falloff")
+
+c     Elliptic solver 
 
-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
-         call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit");
+         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
diff --git a/src/brillq.F b/src/brillq.F
index bbd9220..df6a65d 100644
--- a/src/brillq.F
+++ b/src/brillq.F
@@ -50,93 +50,86 @@ c         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
       implicit none
 
       DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
 
       logical firstcall
 
       integer qtype
-      integer b,c
 
       CCTK_REAL brillq,rho,z,phi
-      CCTK_REAL a,r0,sr,rho0,srho,sz
-      CCTK_REAL d,e,m,n,phi0
 
       data firstcall /.true./
 
       save firstcall,qtype
-      save a,b,c,r0,sr,rho0,srho,sz
-      save d,e,m,n,phi0
 
 c     Get parameters at first call.
 
       if (firstcall) then
 
-         qtype = brill_q
-
-         a = brill_a
-         b = brill_b
-         c = brill_c
-
-         r0   = brill_r0
-         sr   = brill_sr
-         rho0 = brill_rho0
-         srho = brill_srho
-         sz   = brill_sz
-
-         if (axisym.eq.0) then
-            d = brill_d
-            e = brill_e
-            m = brill_m
-            n = brill_n
-            phi0 = brill_phi0
+         if (CCTK_EQUALS(q_function,"exp")) then
+            qtype = 0
+         else if (CCTK_EQUALS(q_function,"eppley")) then
+            qtype = 1
+         else if (CCTK_EQUALS(q_function,"gundlach")) then
+            qtype = 2
+         else
+            call CCTK_WARN(0,"Brill wave data type not recognised")
          end if
 
          firstcall = .false.
 
       end if
 
+      if (rho < 0) then
+         rho = -rho
+      end if
+
 c     Calculate q(rho,z) from a choice of forms.
 
-c     brill_q = 0.
+c     "exp" brill data
 
       if (qtype.eq.0) then
 
-         brillq = a*dabs(rho)**(2 + b)
-     .          / dexp((rho - rho0)**2)/(rho**2 + z**2)
+         brillq = exp_a*dabs(rho)**(2 + exp_b)
+     .          / dexp((rho - exp_rho0)**2)/(rho**2 + z**2)
 
-         if (sz.ne.0.0D0) then
-            brillq = brillq/dexp((z/sz)**2)
+         if (exp_sigmaz.ne.0.0D0) then
+            brillq = brillq/dexp((z/exp_sigmaz)**2)
          end if
 
       else if (qtype.eq.1) then
 
-c     brill_q = 1.  This includes Eppleys choice of q. 
+c     "Eppley" brill data.  This includes Eppleys choice of q. 
 c     But note that q(Eppley) = 2q(Cactus).
 
-         brillq = a*(rho/srho)**b 
-     .          / ( 1.0D0 + ((rho**2 + z**2 - r0**2)/sr**2)**(c/2) )
+         brillq = eppley_a*(rho/eppley_sigmarho)**eppley_b 
+     &          / ( 1.0D0 + ((rho**2 + z**2 - eppley_r0**2)/
+     &          eppley_sigmar**2)**(eppley_c/2) )
 
       else if (qtype.eq.2) then
 
-c     brill_q = 2.   This includes my (Carstens) notion of what a
+c     "Gundlach" brill data   This includes my (Carstens) notion of what a
 c     smooth "pure quadrupole" q should be, plus the choice of
 c     Holz et al.
 
-         brillq = a*(rho/srho)**b
-     .          / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2))
+         brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b
+     &          / dexp(((rho**2 + z**2 - gundlach_r0**2)/
+     &            gundlach_sigmar**2)**(gundlach_c/2))
 
       else
 
 c     Unknown type for q function.
 
-         call CCTK_WARN(0,"Unknown type of Brill function q")
+         call CCTK_WARN(0,"brillq: Unknown type of Brill function q")
 
       end if
 
 c     If desired, multiply with a phi dependent factor.
 
-      if (axisym.eq.0) then
-         brillq = brillq*(1.0D0 + d*rho**m*cos(n*(phi-phi0))**2
-     .          / (1.0D0 + e*rho**m))
+      if (CCTK_EQUALS(initial_data,"brilldata")) then
+         brillq = brillq*(1.0D0 + brill3d_d*rho**brill3d_m*
+     &        cos(brill3D_n*(phi-brill3d_phi0))**2
+     &        / (1.0D0 + brill3d_e*rho**brill3d_m))
       end if
 
       return
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index 9af9fbf..467c28d 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -55,7 +55,7 @@ c     Brill metric calculated from q and Psi.
                y1 = y(i,j,k)
                z1 = z(i,j,k)
 
-               rho2 = x1**2 + y1**2
+               rho2 = x1*x1 + y1*y1
                rho1 = dsqrt(rho2)
 
                phi = phif(x1,y1)
@@ -74,8 +74,8 @@ c              The individual coefficients can be read off as
 
                if (rho1.gt.rhofudge) then
 
-                  gxx(i,j,k) = (e2q + (one - e2q)*y1**2/rho2)
-                  gyy(i,j,k) = (e2q + (one - e2q)*x1**2/rho2)
+                  gxx(i,j,k) = (e2q + (one - e2q)*y1*y1/rho2)
+                  gyy(i,j,k) = (e2q + (one - e2q)*x1*x1/rho2)
                   gzz(i,j,k) = e2q
                   gxy(i,j,k) = - (one - e2q)*x1*y1/rho2
 
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
index f84e825..f2e5f13 100644
--- a/src/setupbrilldata2D.F
+++ b/src/setupbrilldata2D.F
@@ -33,11 +33,8 @@ c       f
       DECLARE_CCTK_PARAMETERS
 
       integer i,j,k
-      integer nx,ny,nz
       integer ierr
 
-      integer, dimension(3) :: sym
-
       CCTK_REAL x1,y1,z1,rho1
       CCTK_REAL brillq,eps
       CCTK_REAL zp,zm,rhop,rhom
@@ -50,53 +47,34 @@ 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     Epsilon for finite differencing.
 
       eps = cctk_delta_space(1)
 
-c     Initialize psi.
-
-      brillpsi = one
-
-c     Initialize metric.
-
-      gxx = one
-      gyy = one
-      gzz = one
-
-      gxy = zero
-      gxz = zero
-      gyz = zero
-
-c     Define the symmetries for the brill GFs.
-
-      sym(1) = 1
-      sym(2) = 1
-      sym(3) = 1
-
-      call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillMlinear')
-      call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillNsource')
-
 c     Set up coefficient arrays for elliptic solve.
 c     Notice that the Cactus conventions are:
 c     __
 c     \/ psi +  Mlinear*psi  +  Nsource  =  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)
 
                x1 = x(i,j,k)
                y1 = y(i,j,k)
                z1 = z(i,j,k)
 
-               rho1 = dsqrt(x1**2 + y1**2)
+               rho1 = dsqrt(x1*x1 + y1*y1)
+
+               brillpsi(i,j,k) = one
+
+               gxx(i,j,k) = one
+               gyy(i,j,k) = one
+               gzz(i,j,k) = one
+
+               gxy(i,j,k) = zero
+               gxz(i,j,k) = zero
+               gyz(i,j,k) = zero
 
 c              Centered derivatives. Note that here we may be calling brillq
 c              with a small negative rho, but that should be ok as long as
@@ -116,17 +94,15 @@ c              will not be regular on the axis.
      .              + brillq(rhom,z1,zero)
      .              - 4.0D0*brillq(rho1,z1,zero))/eps**2
 
+               brillNsource(i,j,k) = zero
+
             end do
          end do
       end do
 
-c     Set coefficient Nsource = 0.
-
-      brillNsource = zero
-
 c     Synchronization and boundaries.
 
-      call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
+c      call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
       call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
 
       return
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
index ecd0fd3..97cf28b 100644
--- a/src/setupbrilldata3D.F
+++ b/src/setupbrilldata3D.F
@@ -43,22 +43,13 @@ c      c        \         z      rho                    phi          phi    /
       DECLARE_CCTK_FUNCTIONS
 
       integer i,j,k
-      integer nx,ny,nz
       integer ierr
 
-      integer, dimension(3) :: sym
-
       CCTK_REAL x1,y1,z1,rho1,rho2
       CCTK_REAL phi,phif,e2q
       CCTK_REAL brillq,rhofudge,eps
       CCTK_REAL zero,one,two,three
 
-c     Set up grid size.
-
-      nx = cctk_lsh(1)
-      ny = cctk_lsh(2)
-      nz = cctk_lsh(3)
-
 c     Define numbers
 
       zero  = 0.0D0
@@ -74,78 +65,44 @@ c     Epsilon for finite differencing.
 
       eps = cctk_delta_space(1)
 
-c     Initialize psi.
-
-      brillpsi = one
-
-      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)
                z1 = z(i,j,k)
 
-               rho2 = x1**2 + y1**2
+               rho2 = x1*x1 + y1*y1
                rho1 = dsqrt(rho2)
-
                phi = phif(x1,y1)
 
                e2q  = dexp(two*brillq(rho1,z1,phi))
 
-               if (rho1.gt.rhofudge) then
-                  gxx(i,j,k) = e2q + (one - e2q)*y1**2/rho2
-                  gyy(i,j,k) = e2q + (one - e2q)*x1**2/rho2
-                  gzz(i,j,k) = e2q
-                  gxy(i,j,k) = - (one - e2q)*x1*y1/rho2
-               else
-                  gxx(i,j,k) = one 
-                  gyy(i,j,k) = one 
-                  gzz(i,j,k) = one
-                  gxy(i,j,k) = zero
-               end if
-
-            end do
-         end do
-      end do
-
-      gxz = zero
-      gyz = zero
-
-c     Define the symmetries for the brill GFs.
-
-      sym(1) = 1
-      sym(2) = 1
-      sym(3) = 1
-
-      call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillMlinear')
-      call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillNsource')
-
-c     Set up coefficient arrays for elliptic solve.
-c     Notice that the Cactus conventions are:
-c     __
-c     \/ psi +  Mlinear*psi  +  Nsource  =  0
-
-      do k=1,nz
-         do j=1,ny
-            do i=1,nx
-
-               x1 = x(i,j,k)
-               y1 = y(i,j,k)
-               z1 = z(i,j,k)
-
-               rho2 = x1**2 + y1**2
-               rho1 = dsqrt(rho2)
+c              Initialise psi
 
-               phi  = phif(x1,y1)
+               brillpsi(i,j,k) = one
 
-               e2q = dexp(two*brillq(rho1,z1,phi))
+c              Set up metric and coefficient arrays for elliptic solve.
+c              Notice that the Cactus conventions are:
+c              __
+c              \/ psi +  Mlinear*psi  +  Nsource  =  0
 
 c              Find M using centered differences over a small
 c              interval.
 
+c              Here we assume that for very small rho, the
+c              phi derivatives are essentially zero.  This
+c              must always be true otherwise the function
+c              is not regular on the axis.
+
                if (rho1.gt.rhofudge) then
 
+                  gxx(i,j,k) = e2q + (one - e2q)*y1*y1/rho2
+                  gyy(i,j,k) = e2q + (one - e2q)*x1*x1/rho2
+                  gzz(i,j,k) = e2q
+                  gxy(i,j,k) = - (one - e2q)*x1*y1/rho2
+
                   brillMlinear(i,j,k) = 0.25D0/e2q
      .                 *(brillq(rho1,z1+eps,phi) 
      .                 + brillq(rho1,z1-eps,phi) 
@@ -154,8 +111,20 @@ c              interval.
      .                 - 4.0D0*brillq(rho1,z1,phi))
      .                 / eps**2
 
+                  brillMlinear(i,j,k) = brillMlinear(i,j,k) + 0.25D0/rho2
+     .               *(three*0.25D0*(brillq(rho1,z1,phi+eps)
+     .               - brillq(rho1,z1,phi-eps))**2
+     .               + two*(brillq(rho1,z1,phi+eps)
+     .               - two*brillq(rho1,z1,phi)
+     .               + brillq(rho1,z1,phi-eps)))/eps**2
+
                else
 
+                  gxx(i,j,k) = one 
+                  gyy(i,j,k) = one 
+                  gzz(i,j,k) = one
+                  gxy(i,j,k) = zero
+
                   brillMlinear(i,j,k) = 0.25D0/e2q
      .                 *(brillq(rho1,z1+eps,phi) 
      .                 + brillq(rho1,z1-eps,phi) 
@@ -165,32 +134,21 @@ c              interval.
 
                end if
 
-c              Here I assume that for very small rho, the
-c              phi derivatives are essentially zero.  This
-c              must always be true otherwise the function
-c              is not regular on the axis.
+               gxz(i,j,k) = zero
+               gyz(i,j,k) = zero
 
-               if (rho1.gt.rhofudge) then
-                  brillMlinear(i,j,k) = brillMlinear(i,j,k) + 0.25D0/rho2
-     .               *(three*0.25D0*(brillq(rho1,z1,phi+eps)
-     .               - brillq(rho1,z1,phi-eps))**2
-     .               + two*(brillq(rho1,z1,phi+eps)
-     .               - two*brillq(rho1,z1,phi)
-     .               + brillq(rho1,z1,phi-eps)))/eps**2
-               end if
+c              Set coefficient Nsource = 0
+               brillNsource(i,j,k) = zero
 
             end do
          end do
       end do
 
-c     Set coefficient Nsource = 0.
-
-      brillNsource = zero
 
 c     Synchronization and boundaries.
 
-      call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
-      call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
+c      call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
+c      call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
 
       return
       end
diff --git a/test/test_brilldata_1.par b/test/test_brilldata_1.par
index 6dc6543..5d673ba 100644
--- a/test/test_brilldata_1.par
+++ b/test/test_brilldata_1.par
@@ -17,13 +17,13 @@ cactus::cctk_itlast 				= 0
 
 # Brill wave initial data
 
-ADMBase::initial_data 				= "brilldata"
+ADMBase::initial_data 				= "brilldata2D"
 ADMBase::metric_type				= "physical"
 
 # Brill wave parameters
 
-idbrilldata::brill_q 				= 2
-idbrilldata::brill_a 				= 5.0
+idbrilldata::q_function 			= "gundlach"	
+idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
diff --git a/test/test_brilldata_2.par b/test/test_brilldata_2.par
index 858a256..630cd64 100644
--- a/test/test_brilldata_2.par
+++ b/test/test_brilldata_2.par
@@ -17,14 +17,14 @@ cactus::cctk_itlast 				= 0
 
 # Brill wave initial data
 
-ADMBase::initial_data 				= "brilldata"
+ADMBase::initial_data 				= "brilldata2D"
 ADMBase::metric_type				= "static conformal"
 StaticConformal::conformal_storage              = "factor"
 
 # Brill wave parameters
 
-idbrilldata::brill_q 				= 2
-idbrilldata::brill_a 				= 5.0
+idbrilldata::q_function 			= "gundlach"
+idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
diff --git a/test/test_wsor.par b/test/test_wsor.par
index f3fa563..a9b0337 100644
--- a/test/test_wsor.par
+++ b/test/test_wsor.par
@@ -1,4 +1,4 @@
-!DESC "Brill wave initial data (Holz type), amplitude 5, SOR, robin boundaries, physical metric"
+!DESC "Brill wave initial data (Physical, SOR)"
 
 ActiveThorns = "admbase staticconformal ellbase ellsor boundary pugh pughslab pughreduce cartgrid3d idbrilldata ioascii iobasic ioutil time"
 
@@ -16,13 +16,13 @@ cactus::cctk_itlast 				= 0
 
 # Brill wave initial data
 
-ADMBase::initial_data 				= "brilldata"
+ADMBase::initial_data 				= "brilldata2D"
 ADMBase::metric_type				= "physical"
 
 # Brill wave parameters
 
-idbrilldata::brill_q 				= 2
-idbrilldata::brill_a 				= 5.0
+idbrilldata::q_function 			= "gundlach"	
+idbrilldata::gundlach_a 			= 5.0
 
 # Elliptic solver.
 
@@ -50,4 +50,3 @@ IOASCII::out1D_vars = "admbase::gxx admbase::gyy admbase::gzz admbase::gxy idbri
 
 IO::new_filename_scheme = "no"
 
-pugh::enable_all_storage = "yes" # fixes a seg fault for now