bug fixes to SOR, still testing

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@33 fa3da13c-9f13-4301-a575-cf5b8c5e1907

1 files changed, 12 insertions, 30 deletions

diff --git a/src/sor_flat.F b/src/sor_flat.F
index d02272f..e7722e5 100644
--- a/src/sor_flat.F
+++ b/src/sor_flat.F
@@ -150,27 +150,14 @@ c      none = contains("sor_accel","none").ne.0
       two     = 2.0D0
       four    = 4.0D0
       resnorm = 0
-
-      ae = 1.
-      aw = 1.
-      an = 1.
-      as = 1.
-      at = 1.
-      ab = 1.
-      ac =-6.
-
-c      ane = 0.            
-c      anw = 0.               
-c      ase = 0.                
-c      asw = 0.               
-c      ate = 0.               
-c      atw = 0.               
-c      abe = 0.               
-c      abw = 0.             
-c      atn = 0.                
-c      ats = 0.               
-c      abn = 0.              
-c      asb = 0. 
+      
+      ae = 1.0d0/dx**2.
+      aw = 1.0d0/dx**2.
+      an = 1.0d0/dy**2.
+      as = 1.0d0/dy**2.
+      at = 1.0d0/dz**2.
+      ab = 1.0d0/dz**2.
+      ac_orig = -2.0d0/dx**2 - 2.0d0/dy**2 - 2.0d0/dz**2
 
       do sor_iteration=1,maxit
 
@@ -203,9 +190,9 @@ c     Start loop with Red Black
          do k=ks,cctk_lsh(3)-1,kstep
             do j=2,cctk_lsh(2)-1
                do i=2,cctk_lsh(1)-1
-                  
-                  ac = -6.0
 
+                  ac = ac_orig                  
+             
                   if (Mlinear_storage.eq.1) then
                      ac = ac - Mlinear(i,j,k)
                   endif
@@ -215,13 +202,6 @@ c     Start loop with Red Black
      &                 + an*var(i,j+1,k) + as*var(i,j-1,k)
      &                 + at*var(i,j,k+1) + ab*var(i,j,k-1)
 
-c     &                 + ane*var(i+1,j+1,k) + anw*var(i-1,j+1,k)
-c     &                 + ase*var(i+1,j-1,k) + asw*var(i-1,j-1,k)
-c     &                 + ate*var(i+1,j,k+1) + atw*var(i-1,j,k+1)
-c     &                 + abe*var(i+1,j,k-1) + abw*var(i-1,j,k-1)
-c     &                 + atn*var(i,j+1,k+1) + ats*var(i,j-1,k+1)
-c     &                 + abn*var(i,j+1,k-1) + asb*var(i,j-1,k-1)
-
                   if (Nsource_storage.eq.1) then
                      residual = residual - Nsource(i,j,k)
                   endif
@@ -241,6 +221,7 @@ c     Reduce the norm
             call CCTK_WARN(1,"Reduction of norm failed!");
          endif
          residual = resnorm
+         residual = residual / (cctk_gsh(1)*cctk_gsh(2)*cctk_gsh(3))
 
 c     Synchronize the variables
          call CCTK_SyncGroupWithVarI(cctkGH, var_idx)
@@ -251,6 +232,7 @@ c     Synchronize the variables
 
 c     Apply boundary conditions
 c         call Ell_GetStringKey(nchar, mybound,"EllLinFlat::Bnd")
+
          if (CCTK_EQUALS(sor_bound,"robin")) then
             call RobinBCVarI(ierr, cctkGH, finf, npow, sw, var_idx);
             if (ierr.ne.0) then