Tidying so that the conformal solver isn't used to solve for the conformal

factor, in fact the conformal factor is now only set at the end of
the thorn.

Now if the conformal metric is used, psi will be set to the conformal factor
which was actually solved for. This means that right now that you can't use IDBrillData with conformal derivatives, although it would presumably be
easy enough to numerically calculate them, or to return 1 and 0's instead
just for this case.


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

5 files changed, 73 insertions, 65 deletions

diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 50a7188..fe55173 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -49,18 +49,26 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
 
   /* Do we know how to deal with this type of metric ? */
-  if( ! CCTK_EQUALS(metric_type, "physical") && 
-      ! CCTK_EQUALS(metric_type, "static conformal"))
+  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\"");
   }
 
+  if (CCTK_Equals(metric_type, "static conformal"))
+  {
+    if (!CCTK_Equals(conformal_storage,"factor"))
+    {
+      CCTK_PARAMWARN("BrillData only sets the conformal factor (this could easily be changed ... please ask)");
+    }
+  }
+
   CCTK_INFO("Setting up Brill data");
-  if (CCTK_EQUALS(metric_type, "static conformal"))
+  if (CCTK_Equals(metric_type, "static conformal"))
   {
      CCTK_VInfo(CCTK_THORNSTRING,"  ... using trivial conformal %s",conformal_storage);
   }
-  else if (CCTK_EQUALS(metric_type, "physical"))
+  else if (CCTK_Equals(metric_type, "physical"))
   {
     CCTK_INFO("  ... using physical metric");
   }
diff --git a/src/brilldata.F b/src/brilldata.F
index 596da0d..7dc01e0 100644
--- a/src/brilldata.F
+++ b/src/brilldata.F
@@ -23,26 +23,25 @@
       DECLARE_CCTK_FUNCTIONS
 
       integer ipsi,iMcoeff,iNcoeff
-      integer metpsi_index(7)
+      integer metric_index(6)
       integer ierr
 
       CCTK_REAL AbsTol(3),RelTol(3)
 
 c     Get indices for metric.
 
-      call CCTK_VarIndex(metpsi_index(1), "admbase::gxx")
-      call CCTK_VarIndex(metpsi_index(2), "admbase::gxy")
-      call CCTK_VarIndex(metpsi_index(3), "admbase::gxz")
-      call CCTK_VarIndex(metpsi_index(4), "admbase::gyy")
-      call CCTK_VarIndex(metpsi_index(5), "admbase::gyz")
-      call CCTK_VarIndex(metpsi_index(6), "admbase::gzz")
-      call CCTK_VarIndex(metpsi_index(7), "staticconformal::psi")
+      call CCTK_VarIndex(metric_index(1), "admbase::gxx")
+      call CCTK_VarIndex(metric_index(2), "admbase::gxy")
+      call CCTK_VarIndex(metric_index(3), "admbase::gxz")
+      call CCTK_VarIndex(metric_index(4), "admbase::gyy")
+      call CCTK_VarIndex(metric_index(5), "admbase::gyz")
+      call CCTK_VarIndex(metric_index(6), "admbase::gzz")
 
 c     Get indices for grid functions.
 
       call CCTK_VarIndex(ipsi,"idbrilldata::brillpsi")
       if (ipsi.lt.0) then
-	call CCTK_WARN(0,"Grid variable index for iphi not found")
+	call CCTK_WARN(0,"Grid variable index for ipsi not found")
       end if
 
       call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear")
@@ -83,14 +82,14 @@ c     Elliptic solver.
 
       if (CCTK_EQUALS(brill_solver,"sor")) then
          call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit");
-         call Ell_LinConfMetricSolver(ierr,cctkGH, 
-     .   metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor")
+         call Ell_LinMetricSolver(ierr,cctkGH, 
+     .   metric_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")
+         call Ell_LinMetricSolver(ierr,cctkGH, 
+     .   metric_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")
+         call Ell_LinMetricSolver(ierr,cctkGH, 
+     .   metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"bam")
       end if
 
 c     Check for errors.
@@ -110,9 +109,9 @@ c     Synchronization and symmetry boundaries.
       call CCTK_SyncGroup(ierr,cctkGH,"idbrilldata::brillconf")
       call CartSymGN(ierr,cctkGH,"idbrilldata::brillconf")
 
-c     Reconstruct physical metric.
+c     Construct final metric.
 
-      call finishbrilldata(CCTK_ARGUMENTS)
+      call IDBrillData_Finish(CCTK_ARGUMENTS)
 
       return
       end
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index 1505d90..9af9fbf 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -12,12 +12,13 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
-      subroutine finishbrilldata(CCTK_ARGUMENTS)
+      subroutine IDBrillData_Finish(CCTK_ARGUMENTS)
 
       implicit none
 
       DECLARE_CCTK_ARGUMENTS
       DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
 
       integer i,j,k
       integer nx,ny,nz
@@ -59,7 +60,6 @@ c     Brill metric calculated from q and Psi.
 
                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,
@@ -74,20 +74,20 @@ c              The individual coefficients can be read off as
 
                if (rho1.gt.rhofudge) then
 
-                  gxx(i,j,k) = psi4*(e2q + (one - e2q)*y1**2/rho2)
-                  gyy(i,j,k) = psi4*(e2q + (one - e2q)*x1**2/rho2)
-                  gzz(i,j,k) = psi4*e2q
-                  gxy(i,j,k) = - psi4*(one - e2q)*x1*y1/rho2
+                  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
 
 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) = psi4
-                  gyy(i,j,k) = psi4
-                  gzz(i,j,k) = psi4
-                  gxy(i,j,k) = zero
+                  gxx(i,j,k) = 0.0d0
+                  gyy(i,j,k) = 0.0d0
+                  gzz(i,j,k) = 0.0d0
+                  gxy(i,j,k) = 0.0d0
 
                end if
 
@@ -95,6 +95,41 @@ c                 it should be, or the data will be singular.
          end do
       end do
 
+      if (CCTK_EQUALS(metric_type,"static conformal")) then
+
+         conformal_state = 1
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  psi(i,j,k) = brillpsi(i,j,k)
+                  
+               end do
+            end do
+         end do
+
+      else
+
+         conformal_state = 0
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  psi4  = brillpsi(i,j,k)**4
+
+                  gxx(i,j,k) = psi4*gxx(i,j,k)
+                  gyy(i,j,k) = psi4*gyy(i,j,k)
+                  gzz(i,j,k) = psi4*gzz(i,j,k)
+                  gxy(i,j,k) = psi4*gxy(i,j,k)
+                  
+               end do
+            end do
+         end do
+
+      end if
+
 c     In any case,
 
       gxz = zero
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
index 98ffe80..f84e825 100644
--- a/src/setupbrilldata2D.F
+++ b/src/setupbrilldata2D.F
@@ -66,8 +66,6 @@ c     Initialize psi.
 
 c     Initialize metric.
 
-      psi = one
-
       gxx = one
       gyy = one
       gzz = one
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
index 9daa65f..ecd0fd3 100644
--- a/src/setupbrilldata3D.F
+++ b/src/setupbrilldata3D.F
@@ -78,38 +78,6 @@ c     Initialize psi.
 
       brillpsi = one
 
-c     Set up conformal metric.
-
-      if (CCTK_EQUALS(metric_type,"static conformal")) then
-      
-         conformal_state = 1
-         do k=1,nz
-            do j=1,ny
-               do i=1,nx
-                  psi(i,j,k) = one
-                  if (CCTK_EQUALS(conformal_storage,"factor+derivs") .or. 
-     &                 CCTK_EQUALS(conformal_storage,"factor+derivs+2nd derivs")) then
-                     psix(i,j,k) = zero
-                     psiy(i,j,k) = zero
-                     psiz(i,j,k) = zero
-                     conformal_state = 2
-                  end if
-                  if (CCTK_EQUALS(conformal_storage,"factor+derivs+2nd derivs")) then
-                     psixx(i,j,k) = zero
-                     psiyy(i,j,k) = zero
-                     psizz(i,j,k) = zero
-                     psixy(i,j,k) = zero
-                     psixz(i,j,k) = zero
-                     psiyz(i,j,k) = zero
-                     conformal_state = 3
-                  end if
-               end do
-            end do
-         end do
-
-      end if
-         
-
       do k=1,nz
          do j=1,ny
             do i=1,nx