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

7 files changed, 144 insertions, 172 deletions

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