Updates

git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWaveElliptic/trunk@19 41e88fdd-2190-4c69-9c84-4659c8cf322e

2 files changed, 39 insertions, 16 deletions

diff --git a/par/charge.par b/par/charge.par
index 0e15d40..1224cff 100644
--- a/par/charge.par
+++ b/par/charge.par
@@ -10,7 +10,7 @@
 #   @enddesc
 # @@*/          
 
-ActiveThorns = "ellsor ellbase idscalarwaveelliptic idscalarwave time wavetoyf77 pugh cartgrid3d ioutil ioascii iobasic"
+ActiveThorns = "ellsor ellbase idscalarwaveelliptic idscalarwave time wavetoyf90 pugh cartgrid3d ioutil ioascii iobasic"
  
 time::dtfac = 0.5
 
@@ -19,16 +19,15 @@ idscalarwaveelliptic::output_tmp = "yes"
 idscalarwaveelliptic::radius = 5.5
 idscalarwaveelliptic::charge = 1
 
-wavetoyf77::bound = "radiation"
+wavetoyf90::bound = "radiation"
 
 grid::type = "BySpacing"
 grid::domain = "octant"
 grid::dxyz = 0.6
 
-driver::global_nx = 60
-driver::global_ny = 60
-driver::global_nz = 60
+driver::global_nsize = 50
 
+ellbase::elliptic_verbose="yes"
 ellsor::maxit = 100
 cactus::cctk_itlast = 120
 
diff --git a/src/SourceData.F77 b/src/SourceData.F77
index bbd786c..5ab3157 100644
--- a/src/SourceData.F77
+++ b/src/SourceData.F77
@@ -26,18 +26,40 @@ c     where Q is the total charge and R is the sphere radius
 
       CCTK_REAL pi
       CCTK_REAL AbsTol(3), RelTol(3)
+      CCTK_REAL charge_factor
 
       integer iphi,iMcoeff,iNcoeff
       integer i,j,k,ierr
-      CCTK_REAL charge_factor
 
-      pi = 4.0*atan(1.0)
 
+c     Get variable indices for all grid functions
+
+      call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
+      if (iMcoeff .lt. 0) then
+	call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
+      end if
+
+      call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
+      if (iNcoeff .lt. 0) then
+	call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
+      end if
+
+      call CCTK_VarIndex (iphi,"wavetoy::phi")
+      if (iphi .lt. 0) then
+	call CCTK_WARN(0,"Grid variable index for iphi not found")
+      end if
+
+
+c     Set up all coefficients and initial guess for solution
+
+      pi = 4.0*atan(1.0)
       charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
 
       do k=1, cctk_lsh(3)
          do j=1, cctk_lsh(2)
             do i=1, cctk_lsh(1)
+
+               phi(i,j,k) = 0.0d0
                      
                Mcoeff(i,j,k) = 0.0d0
 
@@ -50,13 +72,10 @@ c     where Q is the total charge and R is the sphere radius
             end do
          end do
       end do
+      
 
-      phi(i,j,k) = 0.0d0
+c     Set tolerance for stopping elliptic solve
 
-      call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
-      call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
-      call CCTK_VarIndex (iphi,"wavetoy::phi")
-      
       AbsTol(1)=1.0d-5
       AbsTol(2)=1.0d-5
       AbsTol(3)=1.0d-5
@@ -65,22 +84,27 @@ c     where Q is the total charge and R is the sphere radius
       RelTol(2)=-1
       RelTol(3)=-1
 
-c     Call elliptic solver to fill out phi
+
+c     Set any options needed for the elliptic solve
 
 c      call Ell_SetStringKey(ierr, "yes", "EllLinFlat::Bnd::Robin")
 c      write (*,*) "SourceData: SetString ",ierr
-
       call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf")
       call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")
 
 
+c     Call elliptic solver to fill out phi
+
       write (*,*) "charge: Going into elliptic solver"
-      call Ell_LinFlatSolver(ierr,cctkGH, 
+      call Ell_LinFlatSolver(
+     &	   ierr,
+     &     cctkGH, 
      &     iphi, 
      &     iMcoeff, iNcoeff, 
      &     AbsTol, RelTol,
      &     "sor")
-      write (*,*) "charge: Exit elliptic solver"
+      write (*,*) "charge: Exit elliptic solver ierr =",ierr
+
 
 c     Set up last timestep ... assume time symmetry