Developing

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

1 files changed, 33 insertions, 16 deletions

diff --git a/src/SourceData.F b/src/SourceData.F
index d16c0de..2913efb 100644
--- a/src/SourceData.F
+++ b/src/SourceData.F
@@ -16,19 +16,17 @@ c     where Q is the total charge and R is the sphere radius
       DECLARE_CCTK_FARGUMENTS
       DECLARE_CCTK_PARAMETERS
 
-      CCTK_REAL tolerance,pi
+      CCTK_REAL pi
       CCTK_REAL AbsTol(3), RelTol(3)
 
       integer iphi,iMcoeff,iNcoeff
-      integer i,j,k
+      integer i,j,k,ierr
       CCTK_REAL charge_factor
 
       pi = 4.0*atan(1.0)
 
       charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
 
-      tolerance = 1.0d-5
-
       do k=1, cctk_lsh(3)
          do j=1, cctk_lsh(2)
             do i=1, cctk_lsh(1)
@@ -51,8 +49,6 @@ c     where Q is the total charge and R is the sphere radius
       call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
       call CCTK_VarIndex (iphi,"wavetoy::phi")
       
-      write (*,*) iMcoeff, iNcoeff, iphi
-
       AbsTol(1)=1.0d-5
       AbsTol(2)=1.0d-5
       AbsTol(3)=1.0d-5
@@ -62,30 +58,51 @@ c     where Q is the total charge and R is the sphere radius
       RelTol(3)=-1
 
 c     Call to elliptic solver will go here
-c      call Ell_LinConfMetricSolver(cctkGH, 
-c     $     metpsi_index, afield_index, Mlin_index, Nsrc_index, AbsTol, RelTol,
-c     $     "petsc")
+      call Ell_LinFlatSolver(cctkGH, 
+     &     iphi, 
+     &     iMcoeff, iNcoeff, 
+     &     AbsTol, RelTol,
+     &     "sor")
+
+c     Set up last timestep
+      do k=1, cctk_lsh(3)
+         do j=1, cctk_lsh(2)
+            do i=1, cctk_lsh(1)
+                     
+               phi_old(i,j,k) = phi(i,j,k)
+
+            end do
+         end do
+      end do
+
 
-c     Cheat, and use exact solution
+c     Output exact solution if required
       do k=1, cctk_lsh(3)
          do j=1, cctk_lsh(2)
             do i=1, cctk_lsh(1)
                      
                if (r(i,j,k) >= radius) then
-                  phi(i,j,k) = charge/r(i,j,k)
+                  temp(i,j,k) = charge/r(i,j,k)
                else
-                  phi(i,j,k) = charge/(2.0d0*radius**3)*
+                  temp(i,j,k) = charge/(2.0d0*radius**3)*
      &                 (3.0d0*radius**2-r(i,j,k)**2)
                end if
 
-c              Should put in an option for time symmetric first step
-               phi_old(i,j,k) = phi(i,j,k)
-
             end do
          end do
       end do
 
-      
+      if (output_tmp==1) then
+         call CCTK_OutputVarAsByMethod(ierr,cctkGH,
+     &        "idscalarwaveelliptic::temp",
+     &        "IOASCII_1D","phi_exact")
+      end if      
 
       return
       end
+
+
+
+
+
+