Cleaning up. It is not working correctly yet.

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

3 files changed, 66 insertions, 56 deletions

diff --git a/src/brilldata.F b/src/brilldata.F
index 4473166..b811ee4 100644
--- a/src/brilldata.F
+++ b/src/brilldata.F
@@ -4,6 +4,7 @@
 #include "cctk_Arguments.h"
 
 #include "CactusEinstein/Einstein/src/Einstein.h"
+#include "CactusElliptic/EllBase/src/EllBase.h"
 
       subroutine brilldata(CCTK_FARGUMENTS)
 
@@ -13,12 +14,28 @@
       DECLARE_CCTK_PARAMETERS
       DECLARE_CCTK_FUNCTIONS
 
-      integer Mlin_index,Nsrc_index,field_index
-      integer metpsi_index(7)
+      integer ipsi,iMcoeff,iNcoeff
       integer ierr
 
       CCTK_REAL AbsTol(3),RelTol(3)
 
+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")
+      end if
+
+      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_WARN(0,"Grid variable index for Ncoeff not found")
+      end if
+
 c     Set up background metric and coefficients for linear solve.
 
       if (axisym.eq.1) then
@@ -27,64 +44,58 @@ c     Set up background metric and coefficients for linear solve.
          call setupbrilldata3D(CCTK_FARGUMENTS)
       end if
 
-c     Call the Linear Elliptic solver interface to find
-c     conformal factor.  First some preliminaries.
-
-      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")
+c     Tolerances for elliptic solve.
 
       AbsTol(1)= brill_thresh
-      AbsTol(2)= -1
-      AbsTol(3)= -1
+      AbsTol(2)= brill_thresh
+      AbsTol(3)= brill_thresh
  
       RelTol(1)= -1
       RelTol(2)= -1
       RelTol(3)= -1
 
-c     Boundaries.
+c     Elliptic solver.
 
-      if (CCTK_EQUALS(brill_bound,"const")) then
-         call Ell_SetRealKey(ierr,brill_const_v0,
-     .   "EllLinConfMetric::Bnd::Const::V0")
-      end if
+      if (axisym.eq.1) then
 
-      if (CCTK_EQUALS(brill_bound,"robin")) then
-         call Ell_SetIntKey(ierr,brill_robin_falloff,
-     .   "EllLinConfMetric::Bnd::Robin::falloff")
-         call Ell_SetRealKey(ierr,brill_robin_inf,
-     .   "EllLinConfMetric::Bnd::Robin::inf")
-      end if
+         call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf")
+         call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")
 
-c     Elliptic solver.
+         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,"sor")) then
-         call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index,
-     .   field_index,Mlin_index,Nsrc_index,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
+
+         if (CCTK_EQUALS(brill_solver,"bam")) then
+            call Ell_LinFlatSolver(ierr,cctkGH,ipsi,
+     .      iMcoeff,iNcoeff,AbsTol,RelTol,"bam")
+         end if
+
+      else
 
-      if (CCTK_EQUALS(brill_solver,"petsc")) then
-         call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index,
-     .   field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"petsc")
       end if
 
-      if (CCTK_EQUALS(brill_solver,"bam")) then
-         call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index,
-     .   field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"bam")
+c     Check for errors.
+
+      if (ierr .eq. ELL_SUCCESS) then
+        call CCTK_INFO("Leaving elliptic solver: solve successful")
+      else if (ierr .eq. ELL_NOCONVERGENCE) then
+	call CCTK_INFO("Leaving elliptic solver: solver failed to converge")
+      else if (ierr .eq. ELL_NOSOLVER) then
+	call CCTK_INFO("Elliptic solver not found")
+      else 
+        call CCTK_INFO("Leaving elliptic solver: solve failed")	
       end if
 
 c     Synchronization and symmetry boundaries.
 
-      call CCTK_SyncGroup(cctkGH,"IDBrillData::brillconf")
-      call CartSymGN(ierr,cctkGH,"IDBrillData::brillconf")
+      call CCTK_SyncGroup(cctkGH,"idbrilldata::brillconf")
+      call CartSymGN(ierr,cctkGH,"idbrilldata::brillconf")
 
 c     Reconstruct physical metric.
 
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F
index 5643cf0..709370d 100644
--- a/src/setupbrilldata2D.F
+++ b/src/setupbrilldata2D.F
@@ -1,7 +1,9 @@
+
 #include "cctk.h" 
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
-#include "../../../CactusEinstein/Einstein/src/Einstein.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
 
       subroutine setupbrilldata2D(CCTK_FARGUMENTS)
 
@@ -53,15 +55,10 @@ C     Initialize psi.
 
       brillpsi = one
 
-C     Set up flat metric for the elliptic solve.
-C     [Delta(g) + Mlinear] psi = Nsource.
+C     Initialize metric.
 
       psi = one
 
-      psix = zero
-      psiy = zero
-      psiz = zero
-
       gxx = one
       gyy = one
       gzz = one
@@ -87,7 +84,7 @@ 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 
+               brillMlinear(i,j,k) = - 0.25d0 
      $              *(brillq(rho1,z1+eps,zero) 
      $              + brillq(rho1,z1-eps,zero) 
      $              + brillq(rho1+eps,z1,zero) 
@@ -100,7 +97,7 @@ C              will not be regular on the axis.
 
 C     Set coefficient Nsource = 0.
 
-      Nsource = zero
+      brillNsource = zero
 
       return
       end
diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F
index 22ac6d1..0f6abfa 100644
--- a/src/setupbrilldata3D.F
+++ b/src/setupbrilldata3D.F
@@ -1,7 +1,9 @@
+
 #include "cctk.h" 
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
-#include "../../../CactusEinstein/Einstein/src/Einstein.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
 
       subroutine setupbrilldata3D(CCTK_FARGUMENTS)
 
@@ -124,7 +126,7 @@ c              Find M using centered differences over a small
 c              interval.
 
                if (rho1.gt.rhofudge) then
-                  Mlinear(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) 
@@ -132,7 +134,7 @@ c              interval.
      .                 - 4.0D0*brillq(rho1,z1,phi))
      .                 / eps**2
                else
-                  Mlinear(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) 
@@ -146,7 +148,7 @@ 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
+                  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)
@@ -160,7 +162,7 @@ c              is not regular on the axis.
 
 c     Set coefficient Nsource = 0.
 
-      Nsource = zero
+      brillNsource = zero
 
       return
       end