This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

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

13 files changed, 992 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..a6b13a9
--- /dev/null
+++ b/README
@@ -0,0 +1,9 @@
+Cactus Code Thorn BrillData
+Authors    : Carsten Gundlach, Miguel Alcubierre
+Managed by : C.G. <gundlach@aei-potsdam.mpg.de>
+             M.A. <miguel@aei-potsdam.mpg.de>
+--------------------------------------------------------------------------
+
+Initializes Brill wave data.
+
+
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..59a7721
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,142 @@
+% Thorn documentation template
+\documentclass{article}
+\begin{document}
+
+\title{IDBrillData}
+\author{Carsten Gundlach}
+\date{6 September 1999}
+\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.}
+
+\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
+\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},
+\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.
+
+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
+constraints are trivially satisfied.
+
+The thorn is activated by choosing the standard Cactus parameter
+``initial\_data'' in one of the following two ways:
+
+\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}
+
+
+\section{Parameters for the thorn}
+
+The thorn is controlled by the following parameters:
+
+\begin{itemize}
+
+\item brill\_q (INT):  Form of the function $q$ [0,1,2] (default 2):
+
+\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  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}}
+\]
+
+\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 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]
+\]
+
+\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 brill\_e (REAL):  $e$ in above expressions (default 1.0).
+
+\item brill\_m (REAL):  $m$ in above expressions (default 2.0).
+
+\item brill\_n (REAL):  $n$ in above expressions (default 2.0).
+
+\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).
+  
+\item savepsi (KEYWORD): Save conformal factor for output?
+  [``yes'',''no''] Normally, the conformal factor is calculated in the
+  grid function ``psi'', but this is set back to 1 at the end once the
+  physical metric has been constructed.  Setting this parameter to
+  ``yes'' copies the conformal factor to the grid function
+  ``brillpsi0'' before resetting it to 1 (default ``no'').
+
+\end{itemize}
+
+The elliptic solver is controlled by additional the parameters:
+
+\begin{itemize}
+  
+\item solver (KEYWORD): Elliptic solver used to solve the
+  hamiltonian constraint [sor/petsc/bam] (default ``sor'').
+  
+\item thresh (REAL): Threshold for elliptic solver (default
+  0.00001).
+
+\end{itemize}
+
+% Automatically created from the ccl files 
+% Do not worry for now.
+
+\include{interface}
+\include{param}
+\include{schedule}
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..f7b01c8
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,17 @@
+implements: IDBrillData
+inherits: einstein
+
+# Grid functions for linear elliptic equation.
+
+private: 
+
+real linearelliptic TYPE=GF
+{
+  Mlinear,
+  Nsource
+} "Coefficients for linear elliptic equation"
+
+real brillconf TYPE=GF
+{
+  brillpsi
+} "Conformal factor for Brill data"
diff --git a/par/brilldata.par b/par/brilldata.par
new file mode 100644
index 0000000..52a4575
--- /dev/null
+++ b/par/brilldata.par
@@ -0,0 +1,50 @@
+#################################################################
+#
+# brilldata.par
+#
+# Simply Brill wave data of the Hoz et.al. type
+#
+#################################################################
+
+ActiveThorns = "iobasic ellbase ellsor IDBrillData pugh einstein cartgrid3d ioutil ioascii"
+
+# Brill wave initial data
+
+Einstein::initial_data = "brilldata"
+
+# Brill wave parameters [Holz et.al. type data]
+
+idbrilldata::brill_q = 2                  # Form of function q [0,1,2]
+idbrilldata::brill_a = 1.0                # Amplitude (default 0)
+
+idbrilldata::savepsi = "yes"
+
+# Elliptic solver.
+
+idbrilldata::solver = "sor"               #[petsc,sor,bam]
+
+# General
+
+driver::global_nx = 32     
+driver::global_ny = 32     
+driver::global_nz = 32     
+
+grid::type = "byspacing"
+grid::dxyz = 0.2
+grid::domain = "octant" 
+
+cactus::cctk_itlast = 0.0      
+
+# Output.
+
+IO::outdir = "BrillData"
+
+IO::out_every = 1
+
+IOASCII::out1D_vars = "einstein::gxx einstein::gyy einstein::gzz
+                     einstein::gxy einstein::gxz einstein::gyz
+                     IDBrillData: brillpsi0"
+
+IOBasic::outScalar_vars = "einstein::gxx einstein::gyy einstein::gzz
+                     einstein::gxy einstein::gxz einstein::gyz
+                     IDBrillData: brillpsi0"
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..e624a2c
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,119 @@
+# Parameter definitions for thorn IDBrillData
+
+shares:einstein
+
+EXTENDS KEYWORD initial_data ""
+{
+   "brilldata"   :: "Brill wave initial data"
+} ""
+
+
+private:
+
+# Parameters for elliptic solve
+
+KEYWORD solver "Which elliptic solver to use"
+{
+  "sor"    ::  "Use SOR solver"
+  "petsc"  ::  "Use PETSc solver"
+  "bam"    ::  "Use bam solver"
+} "sor"
+
+REAL thresh "How far (absolute norm) to go"
+{
+  0.0: :: "Positive number please"
+} 0.00001
+
+
+# Brill wave parameters
+
+BOOLEAN axisym "Axisymmetric initial data?"
+{
+} "yes"
+
+INT brill_q  "Form of function q [0,1,2]"
+{
+  0:2 :: "Only cases 0,1,2 defined at the moment"
+} 1
+
+REAL brill_a "Brill wave: Amplitude"
+{
+  : ::
+} 0.0
+
+REAL brill_b "Brill wave: rho^b"
+{
+  : :: 
+} 2.0
+
+REAL brill_c "Brill wave: (r^2 - r0^2)^(c/2)"
+{  
+  : ::
+} 2.0
+
+REAL brill_rho0 "Brill wave: radius of torus in rho"
+{
+  : ::
+} 0.0
+
+REAL brill_r0 "Brill wave: radius of torus in r"
+{
+  : ::
+} 0.0
+
+REAL brill_srho "Brill wave: sigma in rho"
+{
+  : ::
+} 1.0
+
+REAL brill_sr "Brill wave: sigma in r"
+{
+  : ::
+} 1.0
+
+# 3D Brill wave parameters
+
+REAL brill_d "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+{
+  : ::
+} 0.0
+
+REAL brill_e "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+{
+  : ::
+} 1.0
+
+REAL brill_m "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+{
+  : ::
+} 2.0
+
+REAL brill_n "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+{
+  : ::
+} 2.0
+
+REAL brill_phi0 "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
+{
+  : ::
+} 0.0
+
+KEYWORD outputpsi "Output conformal factor?"
+{
+ "yes" :: "Output brillpsi"
+ "no"  :: "Don't output brillpsi"
+}"yes"
+
+
+# Additional parameters
+
+REAL brill_eps "epsilon for finite differencing"
+{
+  0: :: "Positive please"
+} 1.e-6
+
+REAL brill_rhofudge "delta rho for axis fudge"
+{
+  0: :: "Positive please"
+} 1.e-6
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..94bd666
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,12 @@
+
+if (CCTK_Equals(initial_data,"brilldata"))
+{
+  schedule brilldata at CCTK_INITIAL
+  {
+    STORAGE: linearelliptic,brillconf
+    COMMUNICATION: linearelliptic,brillconf
+    LANG: Fortran
+  } "Construct Brill wave initial data"
+}
+
+
diff --git a/src/brilldata.F b/src/brilldata.F
new file mode 100644
index 0000000..d7f8348
--- /dev/null
+++ b/src/brilldata.F
@@ -0,0 +1,67 @@
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+      subroutine brilldata(CCTK_FARGUMENTS)
+
+c     Author: Carsten Gundlach.
+c
+c     Driver routine for calculating Brill wave initial data.
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+c     Declare arrays for the linear elliptic solver call:
+
+      integer   Mlin_index,Nsrc_index,field_index,ierr
+      integer   metpsi_index(7)
+
+      CCTK_REAL AbsTol(3), RelTol(3)
+
+c     Set up background metric and coefficients for linear solve.
+
+      if (axisym.eq.1) then
+         call setupbrilldata2D(CCTK_FARGUMENTS)
+      else
+         call setupbrilldata3D(CCTK_FARGUMENTS)
+      end if
+
+c     Call the Linear Elliptic solver interface
+c     to find conformal factor.
+
+      call CCTK_VarIndex (metpsi_index(1), "einstein::gxx")
+      call CCTK_VarIndex (metpsi_index(2), "einstein::gxy")
+      call CCTK_VarIndex (metpsi_index(3), "einstein::gxz")
+      call CCTK_VarIndex (metpsi_index(4), "einstein::gyy")
+      call CCTK_VarIndex (metpsi_index(5), "einstein::gyz")
+      call CCTK_VarIndex (metpsi_index(6), "einstein::gzz")
+      call CCTK_VarIndex (metpsi_index(7), "einstein::psi")
+
+      call CCTK_VarIndex (field_index,     "IDBrillData::brillpsi")
+      call CCTK_VarIndex (Mlin_index,      "IDBrillData::Mlinear")
+      call CCTK_VarIndex (Nsrc_index,      "IDBrillData::Nsource")
+
+      AbsTol(1)= 0.0001
+      AbsTol(2)= -1
+      AbsTol(3)= -1
+ 
+      RelTol(1)= -1
+      RelTol(2)= -1
+      RelTol(3)= -1
+
+      call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index,
+     .field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"sor")
+
+c     Synchronize conformal factor.
+
+      call CCTK_SyncGroup(cctkGH,"einstein::confac")
+
+c     Reconstruct physical metric.
+
+      call finishbrilldata(CCTK_FARGUMENTS)
+
+      return
+      end
diff --git a/src/brillq.F b/src/brillq.F
new file mode 100644
index 0000000..fb28639
--- /dev/null
+++ b/src/brillq.F
@@ -0,0 +1,141 @@
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+      function brillq(rho,z,phi)
+
+C     Authors:  Carsten Gundlach, Miguel Alcubierre.
+C
+C     Calculates the function q that appear in the conformal
+C     metric for Brill waves:
+C
+C     ds^2  =  psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 )
+C
+C     There are three different choices for the form of q depending
+C     on the value of the parameter "brill_q":
+C
+C     brill_q = 0:
+C                                           2
+C                           2+b  - (rho-rho0)      2     2
+C                  q = a rho    e            (z/sz)  /  r
+C
+C
+C     brill_q = 1:  (includes Eppleys form as special case)
+C
+C
+C                                  b               2    2      2  c/2
+C                  q = a (rho/srho)  /  { 1 + [ ( r - r0 ) / sr  ]    }
+C
+C
+C     brill_q = 2:  (includes Holz et al form as special case)
+C
+C                                            2    2      2  c/2 
+C                                  b  - [ ( r - r0 ) / sr  ]
+C                  q = a (rho/srho)  e
+C
+C
+C     If we want 3D initial data, the function q is multiplied by an
+C     additional factor:
+C
+C                           m    2                              m
+C         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
+
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+
+      logical firstcall
+
+      integer CCTK_Equals
+      integer qtype
+
+      real*8 brillq,rho,z,phi
+      real*8 a,b,c,r0,sr,rho0,srho,brill_sz
+      real*8 d,e,m,n,phi0
+
+      data firstcall /.true./
+
+      save firstcall,qtype,a,b,c,r0,sr,rho0,srho,brill_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
+         brill_sz = brill_sz
+
+         if (axisym.eq.0) then
+            d = brill_d
+            e = brill_e
+            m = brill_m
+            n = brill_n
+            phi0 = brill_phi0
+         end if
+
+         firstcall = .false.
+
+      end if
+
+      if (rho.lt.0) write(*,*) "Warning: negative rho in brillq:", rho
+
+C     Calculate q(rho,z) from a choice of forms. Type 0, 1 and 2 are those
+C     of Bruegmann.
+
+C     brill_q = 0.
+
+      if (qtype .eq. 0) then
+
+         brillq = a * rho**(2.d0+b)
+     $        / dexp((rho - rho0)**2) / (rho**2 + z**2) 
+         if (brill_sz .ne. 0.d0) then
+            brillq = brillq / exp((z/brill_sz)**2)
+         end if
+
+      else if (qtype .eq. 1) then
+
+C     brill_q = 1.  This includes Eppleys choice of q. 
+C     But note that q(Eppley) = 2q(Cactus).
+
+         brillq = a * (rho/srho)**b 
+     $        / ( 1.d0 + ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) )
+
+      else if (qtype .eq. 2) then
+
+C     brill_q = 2.   This includes my (Carstens) notion of what a
+C     smooth "pure quadrupole" q should be.
+
+         brillq = a * (rho/srho)**b 
+     $        / dexp( ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) )
+
+      else
+
+C     Unknown type for q function.
+
+         write(*,*) "Unknown type of Brill function q."
+         STOP
+
+      end if
+
+C     If desired, multiply with 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))
+      end if
+
+      return
+      end
+
+
+
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
new file mode 100644
index 0000000..605abb7
--- /dev/null
+++ b/src/finishbrilldata.F
@@ -0,0 +1,112 @@
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+#include "../../../CactusEinstein/Einstein/src/Einstein.h"
+
+      subroutine finishbrilldata(CCTK_FARGUMENTS)
+
+C     Author: Carsten Gundlach.
+C
+C     Set up 3-metric, extrinsic curvature and BM variables
+C     once the conformal factor has been found.
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      integer i,j,k
+      integer nx,ny,nz
+      integer CCTK_Equals
+
+      CCTK_REAL x1,y1,z1,rho1,rho2
+      CCTK_REAL brillq,psi4,e2q,rhofudge
+      CCTK_REAL phi,phif
+      CCTK_REAL zero
+
+      external brillq
+
+C     Numbers.
+
+      zero = 0.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
+
+               x1 = x(i,j,k)
+               y1 = y(i,j,k)
+               z1 = z(i,j,k)
+
+               rho2 = x1**2 + y1**2
+               rho1 = dsqrt(rho2)
+
+               phi = phif(x1,y1)
+
+               psi4 = brillpsi(i,j,k)**4
+               e2q  = dexp(2.d0*brillq(rho1,z1,phi))
+
+C              Fudge division by rho^2 on axis. (Physically, y^/rho^2,
+C              x^2/rho^2 and xy/rho^2 are of course regular.)
+C              Transform Brills form of the physical metric to Cartesian
+C              coordinates via
+C
+C              e^2q (drho^2 + dz^2) + rho^2 dphi^2 =
+C              e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2
+C
+C               The individual coefficients can be read off as
+
+               if (rho1 .gt. rhofudge) then
+
+                  gzz(i,j,k) = psi4*e2q
+                  gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2)
+                  gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2)
+                  gxy(i,j,k) = - psi4*(1.d0 - e2q)*x1*y1/rho2
+
+               else
+
+C              This fudge assumes that q = O(rho^2) near the axis. Which
+C              it should be, or the data will be singular.
+
+                  gzz(i,j,k) = psi4 
+                  gxx(i,j,k) = psi4 
+                  gyy(i,j,k) = psi4 
+                  gxy(i,j,k) = zero
+
+               end if
+
+            end do
+         end do
+      end do
+
+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/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..ad8953d
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# -*-Makefile-*-
+# Sub-makefile for thorn IDBrillData
+
+SRCS = brilldata.F setupbrilldata2D.F setupbrilldata3D.F \
+       brillq.F phif.F finishbrilldata.F
+
+
+
diff --git a/src/phif.F b/src/phif.F
new file mode 100644
index 0000000..c24f2ca
--- /dev/null
+++ b/src/phif.F
@@ -0,0 +1,47 @@
+#include "cctk.h"
+
+      CCTK_REAL function phif(x,y)
+
+c     Author: Miguel Alcubierre, October 1998.
+c
+c     Function to find angle phi from coordinates {x,y}.
+
+      implicit none
+
+      CCTK_REAL x,y
+      real*8 zero,half,one,two,pi
+
+c     Define numbers.
+
+      zero = 0.0D0
+      half = 0.5D0
+      one  = 1.0D0
+      two  = 2.0D0
+
+      pi = acos(-one)
+
+c     Find angle between 0 and pi/2 such that tan(phi) = |y/x|.
+
+      if (dabs(x).gt.dabs(y)) then
+         phif = atan(dabs(y/x))
+      else if (dabs(x).lt.dabs(y)) then
+         phif = half*pi - atan(dabs(x/y))
+      else
+         phif = 0.25D0*pi
+      end if
+
+c     Use signs of {x,y} to move to correct quadrant.
+
+      if ((x.eq.zero).and.(y.eq.zero)) then
+         phif = zero
+      else if ((x.le.zero).and.(y.ge.zero)) then
+         phif = pi - phif
+      else if ((x.le.zero).and.(y.le.zero)) then
+         phif = pi + phif
+      else if ((x.ge.zero).and.(y.le.zero)) then
+         phif = two*pi - phif
+      end if
+
+      return
+      end
+
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
new file mode 100644
index 0000000..2c4b351
--- /dev/null
+++ b/src/setupbrilldata2D.F
@@ -0,0 +1,101 @@
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+#include "../../../CactusEinstein/Einstein/src/Einstein.h"
+
+      subroutine setupbrilldata2D(CCTK_FARGUMENTS)
+
+C     Author: Carsten Gundlach.
+C
+C     Set up axisymmetric Brill data for elliptic solve. The elliptic
+C     equation we need to solve is:
+C
+C     __                 2      2
+C     \/  psi  =  psi ( d q  + d   q ) / 4
+C       f                z      rho
+C
+C     where:
+C
+C     __
+C     \/  =  Flat space Laplacian.
+C       f
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      integer i,j,k
+      integer nx,ny,nz
+
+      CCTK_REAL x1,y1,z1,rho1
+      CCTK_REAL brillq,eps
+      CCTK_REAL zero,one
+
+      external brillq
+
+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.
+
+      eps = brill_eps
+
+C     Initialize psi.
+
+      brillpsi = one
+
+C     Set up flat metric for the elliptic solve.
+C     [Delta(g) + Mlinear] psi = Nsource.
+
+      psi = one
+
+      gxx = one
+      gyy = one
+      gzz = one
+      gxy = zero
+      gxz = zero
+      gyz = zero
+
+C     Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho).
+
+      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)
+
+               rho1 = dsqrt(x1**2 + y1**2)
+
+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
+C              brillq is even in rho - physically it must be, or the data
+C              will not be regular on the axis.
+
+               Mlinear(i,j,k) = - 0.25d0 
+     $              *(brillq(rho1,z1+eps,zero) 
+     $              + brillq(rho1,z1-eps,zero) 
+     $              + brillq(rho1+eps,z1,zero) 
+     $              + brillq(rho1-eps,z1,zero) 
+     $              - 4.d0*brillq(rho1,z1,zero))/eps**2
+
+            end do
+         end do
+      end do
+
+C     Set coefficient Nsource = 0.
+
+      Nsource = zero
+
+      return
+      end
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
new file mode 100644
index 0000000..719528b
--- /dev/null
+++ b/src/setupbrilldata3D.F
@@ -0,0 +1,167 @@
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+#include "../../../CactusEinstein/Einstein/src/Einstein.h"
+
+      subroutine setupbrilldata3D(CCTK_FARGUMENTS)
+
+c     Author: Miguel Alcubierre, October 1998.
+c
+c     Set up 3D Brill data for elliptic solve.  The elliptic
+c     equation we need to solve is:
+c
+c     __
+c     \/  psi  =  psi R  /  8
+c       c              c
+c
+c     where:
+c     __
+c     \/  =  Laplacian operator for conformal metric.
+c       c
+c
+c     R   =  Ricci scalar for conformal metric. 
+c      c
+c
+c     The Ricci scalar for the conformal metric turns out to be:
+c
+c             /  -2q    2      2             -2            2       2      \
+c     R  =  2 | e    ( d q  + d   q )  +  rho   ( 3 (d   q)  +  2 d   q ) |
+c      c      \         z      rho                    phi          phi    /
+
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      integer CCTK_Equals
+      integer i,j,k
+      integer nx,ny,nz
+
+      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
+      one   = 1.0D0
+      two   = 2.0D0
+      three = 3.0D0
+
+c     Parameters.
+
+      rhofudge = brill_rhofudge
+      eps      = brill_eps
+
+c     Initialize psi.
+
+      brillpsi = one
+
+c     Set up conformal metric.
+
+      psi = one
+
+      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)
+
+               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     Set up coefficient Mlinear = - (1/8) Rc.
+
+      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)
+
+               phi  = phif(x1,y1)
+
+               e2q = dexp(two*brillq(rho1,z1,phi))
+
+c              Find M using centered differences over a small
+c              interval.
+
+               if (rho1.gt.rhofudge) then
+                  Mlinear(i,j,k) = - 0.25D0/e2q
+     .                 *(brillq(rho1,z1+eps,phi) 
+     .                 + brillq(rho1,z1-eps,phi) 
+     .                 + brillq(rho1+eps,z1,phi) 
+     .                 + brillq(rho1-eps,z1,phi) 
+     .                 - 4.0D0*brillq(rho1,z1,phi))
+     .                 / eps**2
+               else
+                  Mlinear(i,j,k) = - 0.25D0/e2q
+     .                 *(brillq(rho1,z1+eps,phi) 
+     .                 + brillq(rho1,z1-eps,phi) 
+     .                 + two*brillq(eps,z1,phi) 
+     .                 - two*brillq(rho1,z1,phi))
+     .                 / eps**2
+               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.
+
+               if (rho1.gt.rhofudge) then
+                  Mlinear(i,j,k) = Mlinear(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
+
+            end do
+         end do
+      end do
+
+c     Set coefficient Nsource = 0.
+
+      Nsource = zero
+
+      return
+      end
+