Cleaning up and fixing the sign convention. It works!

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

6 files changed, 190 insertions, 101 deletions

diff --git a/src/brilldata.F b/src/brilldata.F
index b811ee4..04c5eeb 100644
--- a/src/brilldata.F
+++ b/src/brilldata.F
@@ -1,3 +1,12 @@
+/*@@
+  @file      brilldata.F
+  @date
+  @author    Carsten Gundlach (Cactus 4, Miguel Alcubierre)
+  @desc
+             Construct Brill wave initial data.
+  @enddesc
+  @version   $Header$
+@@*/
 
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -15,24 +24,35 @@
       DECLARE_CCTK_FUNCTIONS
 
       integer ipsi,iMcoeff,iNcoeff
+      integer metpsi_index(7)
       integer ierr
 
       CCTK_REAL AbsTol(3),RelTol(3)
 
+c     Get indices for metric.
+
+      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")
+
 c     Get indices for grid functions.
 
-      call CCTK_VarIndex(ipsi   , "idbrilldata::brillpsi")
-      if (ipsi .lt. 0) then
+      call CCTK_VarIndex(ipsi,"idbrilldata::brillpsi")
+      if (ipsi.lt.0) then
 	call CCTK_WARN(0,"Grid variable index for iphi not found")
       end if
 
-      call CCTK_VarIndex(iMcoeff, "idbrilldata::brillMlinear")
-      if (iMcoeff .lt. 0) then
+      call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear")
+      if (iMcoeff.lt.0) then
 	call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
       end if
 
-      call CCTK_VarIndex(iNcoeff, "idbrilldata::brillNsource")
-      if (iNcoeff .lt. 0) then
+      call CCTK_VarIndex(iNcoeff,"idbrilldata::brillNsource")
+      if (iNcoeff.lt.0) then
 	call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
       end if
 
@@ -54,39 +74,31 @@ c     Tolerances for elliptic solve.
       RelTol(2)= -1
       RelTol(3)= -1
 
-c     Elliptic solver.
-
-      if (axisym.eq.1) then
-
-         call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf")
-         call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")
-
-         if (CCTK_EQUALS(brill_solver,"sor")) then
-            call Ell_LinFlatSolver(ierr,cctkGH,ipsi,
-     .      iMcoeff,iNcoeff,AbsTol,RelTol,"sor")
-         end if
-
-         if (CCTK_EQUALS(brill_solver,"petsc")) then
-            call Ell_LinFlatSolver(ierr,cctkGH,ipsi,
-     .      iMcoeff,iNcoeff,AbsTol,RelTol,"petsc")
-         end if
+c     Boundaries.
 
-         if (CCTK_EQUALS(brill_solver,"bam")) then
-            call Ell_LinFlatSolver(ierr,cctkGH,ipsi,
-     .      iMcoeff,iNcoeff,AbsTol,RelTol,"bam")
-         end if
+      call Ell_SetRealKey(ierr,1.0d0,"EllLinConfMetric::Bnd::Robin::inf")
+      call Ell_SetIntKey(ierr,1,"EllLinConfMetric::Bnd::Robin::falloff")
 
-      else
+c     Elliptic solver.
 
+      if (CCTK_EQUALS(brill_solver,"sor")) then
+         call Ell_LinConfMetricSolver(ierr,cctkGH, 
+     .   metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor")
+      else if (CCTK_EQUALS(brill_solver,"petsc")) then
+         call Ell_LinConfMetricSolver(ierr,cctkGH, 
+     .   metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"petsc")
+      else if (CCTK_EQUALS(brill_solver,"bam")) then
+         call Ell_LinConfMetricSolver(ierr,cctkGH, 
+     .   metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"bam")
       end if
 
 c     Check for errors.
 
-      if (ierr .eq. ELL_SUCCESS) then
+      if (ierr.eq.ELL_SUCCESS) then
         call CCTK_INFO("Leaving elliptic solver: solve successful")
-      else if (ierr .eq. ELL_NOCONVERGENCE) then
+      else if (ierr.eq.ELL_NOCONVERGENCE) then
 	call CCTK_INFO("Leaving elliptic solver: solver failed to converge")
-      else if (ierr .eq. ELL_NOSOLVER) then
+      else if (ierr.eq.ELL_NOSOLVER) then
 	call CCTK_INFO("Elliptic solver not found")
       else 
         call CCTK_INFO("Leaving elliptic solver: solve failed")	
diff --git a/src/brillq.F b/src/brillq.F
index c7010f6..8d62064 100644
--- a/src/brillq.F
+++ b/src/brillq.F
@@ -1,11 +1,19 @@
+/*@@
+  @file      brilldata.F
+  @date
+  @author    Carsten Gundlach, Miguel Alcubierre.
+  @desc
+             Construct Brill wave q function.
+  @enddesc
+  @version   $Header$
+@@*/
+
 #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
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index ee49b57..36060e9 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -1,3 +1,12 @@
+/*@@
+  @file      finishbrilldata.F
+  @date
+  @author    Carsten Gundlach (Cactus 4, Miguel Alcubierre)
+  @desc
+             Reconstruct metric from conformal factor.
+  @enddesc
+  @version   $Header$
+@@*/
 
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -22,22 +31,22 @@
 
       external brillq
 
-C     Numbers.
+c     Numbers.
 
       zero = 0.0D0
 
-C     Set up grid size.
+c     Set up grid size.
 
       nx = cctk_lsh(1)
       ny = cctk_lsh(2)
       nz = cctk_lsh(3)
 
-C     Parameters.
+c     Parameters.
 
       rhofudge = brill_rhofudge
       
-C     Replace flat metric left over from elliptic solve by
-C     Brill metric calculated from q and Psi.
+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
@@ -55,31 +64,31 @@ C     Brill metric calculated from q and Psi.
                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
+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)
+                  gzz(i,j,k) = psi4*e2q
                   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.
+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 
+                  gxx(i,j,k) = psi4
+                  gyy(i,j,k) = psi4
+                  gzz(i,j,k) = psi4
                   gxy(i,j,k) = zero
 
                end if
@@ -88,12 +97,12 @@ C                 it should be, or the data will be singular.
          end do
       end do
 
-C     In any case,
+c     In any case,
 
       gxz = zero
       gyz = zero
 
-C     Vanishing extrinsic curvature completes the Cauchy data.
+c     Vanishing extrinsic curvature completes the Cauchy data.
 
       kxx = zero
       kyy = zero
diff --git a/src/phif.F b/src/phif.F
index c24f2ca..3d9bf79 100644
--- a/src/phif.F
+++ b/src/phif.F
@@ -1,17 +1,23 @@
+/*@@
+  @file      phif.F
+  @date      October 1998
+  @author    Miguel Alcubierre
+  @desc
+             Function to find angle phi from coordinates {x,y}.
+  @enddesc
+  @version   $Header$
+@@*/
+
 #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
+      CCTK_REAL zero,half,one,two,pi
 
-c     Define numbers.
+c     Numbers.
 
       zero = 0.0D0
       half = 0.5D0
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
index 709370d..a7b35ff 100644
--- a/src/setupbrilldata2D.F
+++ b/src/setupbrilldata2D.F
@@ -1,3 +1,12 @@
+/*@@
+  @file      setupbrilldata2D.F
+  @date
+  @author    Carsten Gundlach (Cactus 4, Miguel Alcubierre)
+  @desc
+             Set up axisymmetric Brill data for elliptic solve.
+  @enddesc
+  @version   $Header$
+@@*/
 
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -7,20 +16,18 @@
 
       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
+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  =  0
+c       f                z      rho
+c
+c     where:
+c
+c     __
+c     \/  =  Flat space Laplacian.
+c       f
 
       implicit none
 
@@ -29,6 +36,9 @@ C       f
 
       integer i,j,k
       integer nx,ny,nz
+      integer ierr
+
+      integer, dimension(3) :: sym
 
       CCTK_REAL x1,y1,z1,rho1
       CCTK_REAL brillq,eps
@@ -36,26 +46,26 @@ C       f
 
       external brillq
 
-C     Numbers.
+c     Numbers.
 
       zero = 0.0D0
       one  = 1.0D0
 
-C     Set up grid size.
+c     Set up grid size.
 
       nx = cctk_lsh(1)
       ny = cctk_lsh(2)
       nz = cctk_lsh(3)
 
-C     Parameters.
+c     Parameters.
 
       eps = brill_eps
 
-C     Initialize psi.
+c     Initialize psi.
 
       brillpsi = one
 
-C     Initialize metric.
+c     Initialize metric.
 
       psi = one
 
@@ -67,7 +77,19 @@ C     Initialize metric.
       gxz = zero
       gyz = zero
 
-C     Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho).
+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
@@ -79,25 +101,30 @@ C     Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho).
 
                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.
+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.
 
-               brillMlinear(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
+               brillMlinear(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.
+c     Set coefficient Nsource = 0.
 
       brillNsource = zero
 
+c     Synchronization and boundaries.
+
+      call CCTK_SyncGroup(cctkGH,'idbrilldata::brillelliptic')
+      call ADM_Boundary(cctkGH,'idbrilldata::brillelliptic')
+
       return
       end
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
index 0f6abfa..7b289b1 100644
--- a/src/setupbrilldata3D.F
+++ b/src/setupbrilldata3D.F
@@ -1,3 +1,12 @@
+/*@@
+  @file      setupbrilldata3D.F
+  @date      October 1998
+  @author    Miguel Alcubierre
+  @desc
+             Set up non-axisymmetric Brill data for elliptic solve.
+  @enddesc
+  @version   $Header$
+@@*/
 
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -7,13 +16,11 @@
 
       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     \/  psi  -  psi R  /  8  =  0
 c       c              c
 c
 c     where:
@@ -26,9 +33,9 @@ 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    /
+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
@@ -39,6 +46,9 @@ c      c      \         z      rho                    phi          phi    /
       integer CCTK_Equals
       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
@@ -105,7 +115,19 @@ c     Set up conformal metric.
       gxz = zero
       gyz = zero
 
-c     Set up coefficient Mlinear = - (1/8) Rc.
+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
@@ -126,7 +148,7 @@ c              Find M using centered differences over a small
 c              interval.
 
                if (rho1.gt.rhofudge) then
-                  brillMlinear(i,j,k) = - 0.25D0/e2q
+                  brillMlinear(i,j,k) = 0.25D0/e2q
      .                 *(brillq(rho1,z1+eps,phi) 
      .                 + brillq(rho1,z1-eps,phi) 
      .                 + brillq(rho1+eps,z1,phi) 
@@ -134,7 +156,7 @@ c              interval.
      .                 - 4.0D0*brillq(rho1,z1,phi))
      .                 / eps**2
                else
-                  brillMlinear(i,j,k) = - 0.25D0/e2q
+                  brillMlinear(i,j,k) = 0.25D0/e2q
      .                 *(brillq(rho1,z1+eps,phi) 
      .                 + brillq(rho1,z1-eps,phi) 
      .                 + two*brillq(eps,z1,phi) 
@@ -148,7 +170,7 @@ c              must always be true otherwise the function
 c              is not regular on the axis.
 
                if (rho1.gt.rhofudge) then
-                  brillMlinear(i,j,k) = brillMlinear(i,j,k) - 0.25D0/rho2
+                  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)
@@ -164,6 +186,11 @@ c     Set coefficient Nsource = 0.
 
       brillNsource = zero
 
+c     Synchronization and boundaries.
+
+      call CCTK_SyncGroup(cctkGH,'idbrilldata::brillelliptic')
+      call ADM_Boundary(cctkGH,'idbrilldata::brillelliptic')
+
       return
       end