From d8f42b55c08565f64087d7c1fe54cff04207ddde Mon Sep 17 00:00:00 2001
From: allen <allen@fa3da13c-9f13-4301-a575-cf5b8c5e1907>
Date: Fri, 13 Apr 2001 14:30:23 +0000
Subject: Fixed a number of problems, notably that timelevels weren't being
 used properly, and error codes not returned.

Still some testing to do, then I'll add a testsuite for SOR.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@69 fa3da13c-9f13-4301-a575-cf5b8c5e1907
---
 src/sor_confmetric.c | 310 +++++++++++++++++++++++++++++++--------------------
 src/sor_flat.c       | 145 ++++++++++++++++++------
 src/sor_wrapper.c    |  76 +++++++------
 3 files changed, 346 insertions(+), 185 deletions(-)

diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c
index 4f5a2e4..6fb677b 100644
--- a/src/sor_confmetric.c
+++ b/src/sor_confmetric.c
@@ -16,13 +16,14 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#include "CactusBase/Boundary/src/Boundary.h"
-#include "CactusBase/CartGrid3D/src/Symmetry.h"
-#include "CactusElliptic/EllBase/src/Ell_DBstructure.h"
+#include "Boundary.h"
+#include "Symmetry.h"
+#include "Ell_DBstructure.h"
+#include "EllBase.h"
 
 #define SQR(a) ((a)*(a))
 
-void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
+int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
 		       int FieldIndex, int MIndex, int NIndex,
 		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
 
@@ -42,14 +43,15 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
 
 @@*/
 
-
-
-void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
+int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
 		       int FieldIndex, int MIndex, int NIndex,
 		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
 {
   DECLARE_CCTK_PARAMETERS  
 
+  int retval = ELL_SUCCESS;
+  int timelevel;
+
   /* The pointer to the metric fields */
   CCTK_REAL *gxx =NULL, *gxy =NULL; 
   CCTK_REAL *gxz =NULL, *gyy =NULL; 
@@ -112,12 +114,15 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
   /* Get the reduction handle */
   sum_handle = CCTK_ReductionArrayHandle("sum");
   if (sum_handle<0) 
+  {
     CCTK_WARN(1,"Cannot get reduction handle for operation >sum<");
-  
+  }
+
   /* IF Robin BCs are set, prepare for a boundary call:
      setup stencil width and get Robin constants (set by the routine
      which is calling the solver interface) */
-  if (CCTK_EQUALS(sor_bound,"robin")) { 
+  if (CCTK_EQUALS(sor_bound,"robin")) 
+  { 
     sw[0]=1; 
     sw[1]=1; 
     sw[2]=1;
@@ -129,21 +134,32 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
   /* Only supports absolute tolerance */
   tol   = AbsTol[0];
   if (CCTK_EQUALS(sor_accel,"const"))
+  {
     accel_const = 1;
+  }
   else if (CCTK_EQUALS(sor_accel,"cheb"))
+  {
     accel_cheb  = 1;
+  }
 
   /* Things to do only once! */
-  if (firstcall==1) {
+  if (firstcall==1) 
+  {
     if (CCTK_Equals(elliptic_verbose, "yes"))
+    {
+      if (accel_cheb)
+      {
+        printf("SOR with Chebyshev acceleration with radius of 1\n");
+      }
+      else if (accel_const)
+      {
+        printf("SOR with hardcoded omega = 1.8\n");
+      }
+      else
       {
-	if (accel_cheb)
-	  printf("SOR with Chebyshev acceleration with radius of 1\n");
-	else if (accel_const)
-	  printf("SOR with hardcoded omega = 1.8\n");
-	else
-	  printf("SOR with unaccelearted relaxation (omega = 1)\n");
+        printf("SOR with unaccelearted relaxation (omega = 1)\n");
       }
+    }
     firstcall = 0;
   }
    
@@ -151,27 +167,56 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
   /* Get the data ptr of these GFs, They all have to be
      on the same timelevel; derive the metric data pointer from 
      the index array. Note the ordering in the metric   */
-  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex);
-  gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[0]);
-  gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[1]);
-  gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[2]);
-  gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[3]);
-  gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[4]);
-  gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[5]);
-  
+  timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[0]) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+  gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[0]);
+  gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[1]);
+  gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[2]);
+  gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[3]);
+  gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[4]);
+  gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[5]);
 
   if (conformal)
-    psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[6]);
+  {
+    timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[6]) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[6]);
+  }
 
   /* if we have a negative index for M/N, this GF is not needed, 
      there for don't even look for it. when index positive,
      set flag Mstorage=1, dito for N */
-  if (MIndex>=0)  { 
-    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
+  if (MIndex>=0)  
+  { 
+    timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
     Mstorage = 1;
   }
-  if (NIndex>=0) {
-    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex);
+  if (NIndex>=0) 
+  {
+    timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
     Nstorage = 1;
   }
 
@@ -201,7 +246,8 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
 
   /* Calculate the inverse metric */
 
-  for (ijk=0;ijk<nxyz;ijk++) {
+  for (ijk=0;ijk<nxyz;ijk++) 
+  {
 
     /* determinant */
     det = -(SQR(gxz[ijk])*gyy[ijk]) + 
@@ -210,10 +256,13 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
       SQR(gxy[ijk])*gzz[ijk] + 
       gxx[ijk]*gyy[ijk]*gzz[ijk];
     
-    if (conformal) {
+    if (conformal) 
+    {
       pm4 = 1.0/pow(psi[ijk],4.0);
       p12 = pow(psi[ijk],12.0);
-    } else {
+    } 
+    else 
+    {
       pm4 = 1.0;
       p12 = 1.0;
     }
@@ -247,8 +296,10 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
   }	
 
   /*$for (k=1;k<GH->cctk_lsh[2]-1;k++) {
-    for (j=1;j<GH->cctk_lsh[1]-1;j++)   {
-      for (i=1;i<GH->cctk_lsh[0]-1;i++)   {
+    for (j=1;j<GH->cctk_lsh[1]-1;j++)   
+    {
+      for (i=1;i<GH->cctk_lsh[0]-1;i++)   
+      {
 	ijk   = CCTK_GFINDEX3D(GH,i,j,k); 
 	printf("U G : %d %d %d   %7.8g  %7.8g  %7.8g \n",
 	       i,j,k,uxx[ijk],uyy[ijk],gxx[ijk]);
@@ -272,99 +323,111 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
     rjacobian =  1.0;
 
     if (accel_cheb) 
+    {
       for (chebit=2;chebit<sorit;chebit++)
+      {
 	omega = 1.0/(1.0 - 0.25*rjacobian*rjacobian*omega);
+      }
+    }
     if (accel_const)
+    {
       omega = 1.8;
-    
+    }
+
     resnorm = 0.0;
 
     ks = (sorit%2)+1;
     if (GH->cctk_lsh[2]==3)
       ks = 1;
 
-    for (k=ks;k<ke;k+=kstep) {
-      for (j=js;j<je;j++)     {
-	for (i=is;i<ie;i++)    {
-
-	    ijk   = CCTK_GFINDEX3D(GH,i  ,j  ,k  );
-
-	    ipjk  = CCTK_GFINDEX3D(GH,i+1,j  ,k  );
-	    imjk  = CCTK_GFINDEX3D(GH,i-1,j  ,k  );
-	    ijpk  = CCTK_GFINDEX3D(GH,i  ,j+1,k  );
-	    ijmk  = CCTK_GFINDEX3D(GH,i  ,j-1,k  );
-	    ijkp  = CCTK_GFINDEX3D(GH,i  ,j  ,k+1);
-	    ijkm  = CCTK_GFINDEX3D(GH,i  ,j  ,k-1);
-
-	    ipjpk = CCTK_GFINDEX3D(GH,i+1,j+1,k  );
-	    ipjmk = CCTK_GFINDEX3D(GH,i+1,j-1,k  );
-	    imjpk = CCTK_GFINDEX3D(GH,i-1,j+1,k  );
-	    imjmk = CCTK_GFINDEX3D(GH,i-1,j-1,k  );
-
-	    ijpkp = CCTK_GFINDEX3D(GH,i  ,j+1,k+1);
-	    ijpkm = CCTK_GFINDEX3D(GH,i  ,j+1,k-1);
-	    ijmkp = CCTK_GFINDEX3D(GH,i  ,j-1,k+1);
-	    ijmkm = CCTK_GFINDEX3D(GH,i  ,j-1,k-1);
-
-	    ipjkp = CCTK_GFINDEX3D(GH,i+1,j  ,k+1);
-	    ipjkm = CCTK_GFINDEX3D(GH,i+1,j  ,k-1);
-	    imjkp = CCTK_GFINDEX3D(GH,i-1,j  ,k+1);
-	    imjkm = CCTK_GFINDEX3D(GH,i-1,j  ,k-1);
-	    
-
-	    ac = -1.0*uxx[ipjk] - 2.0*uxx[ijk] - 1.0*uxx[imjk]
-	         -1.0*uyy[ijpk] - 2.0*uyy[ijk] - 1.0*uyy[ijmk]
-	         -1.0*uzz[ijkp] - 2.0*uzz[ijk] - 1.0*uzz[ijkm];
-
-
-	    if (Mstorage) 
-	      ac += Mlin[ijk];
-
-
-	    ae  = uxx[ipjk]+uxx[ijk];
-	    aw  = uxx[imjk]+uxx[ijk];
-	    an  = uyy[ijpk]+uyy[ijk];
-	    as  = uyy[ijmk]+uyy[ijk];
-	    at  = uzz[ijkp]+uzz[ijk];
-	    ab  = uzz[ijkm]+uzz[ijk];
-
-	    ane = uxy[ijpk]+uxy[ipjk];
-	    anw =-uxy[imjk]-uxy[ijpk];
-	    ase =-uxy[ipjk]-uxy[ijmk];
-	    asw = uxy[imjk]+uxy[ijmk];
-
-	    ate = uxz[ijkp]+uxz[ipjk];
-	    atw =-uxz[imjk]-uxz[ijkp];
-	    abe =-uxz[ipjk]-uxz[ijkm];
-	    abw = uxz[imjk]+uxz[ijkm];
-	    
-	    atn = uyz[ijpk]+uyz[ijkp];               
-	    ats =-uyz[ijkp]-uyz[ijmk];
-	    abn =-uyz[ijkm]-uyz[ijpk];
-	    absol = uyz[ijkm]+uyz[ijmk];
-
-	    residual = ac * var[ijk]
-	      + ae *var[ipjk]  +  aw*var[imjk]
-	      + an *var[ijpk]  +  as*var[ijmk]
-	      + at *var[ijkp]  +  ab*var[ijkm];
-
-	    residual = residual 
-	      + ane*var[ipjpk] + anw*var[imjpk]; 
-	    
-	    residual = residual 
-	      + ase*var[ipjmk] + asw*var[imjmk];
-
-	    residual = residual    
-	      + ate*var[ipjkp] + atw*var[imjkp] 
-	      + abe*var[ipjkm] + abw*var[imjkm]
-	      + atn*var[ijpkp] + ats*var[ijmkp]
-	      + abn*var[ijpkm] + absol*var[ijmkm];
-	    
-	    if (Nstorage)
-	      residual +=Nlin[ijk];
+    for (k=ks;k<ke;k+=kstep) 
+    {
+      for (j=js;j<je;j++)     
+      {
+	for (i=is;i<ie;i++)    
+        {
+
+	  ijk   = CCTK_GFINDEX3D(GH,i  ,j  ,k  );
+	  
+	  ipjk  = CCTK_GFINDEX3D(GH,i+1,j  ,k  );
+	  imjk  = CCTK_GFINDEX3D(GH,i-1,j  ,k  );
+	  ijpk  = CCTK_GFINDEX3D(GH,i  ,j+1,k  );
+	  ijmk  = CCTK_GFINDEX3D(GH,i  ,j-1,k  );
+	  ijkp  = CCTK_GFINDEX3D(GH,i  ,j  ,k+1);
+	  ijkm  = CCTK_GFINDEX3D(GH,i  ,j  ,k-1);
+	  
+	  ipjpk = CCTK_GFINDEX3D(GH,i+1,j+1,k  );
+	  ipjmk = CCTK_GFINDEX3D(GH,i+1,j-1,k  );
+	  imjpk = CCTK_GFINDEX3D(GH,i-1,j+1,k  );
+	  imjmk = CCTK_GFINDEX3D(GH,i-1,j-1,k  );
+	  
+	  ijpkp = CCTK_GFINDEX3D(GH,i  ,j+1,k+1);
+	  ijpkm = CCTK_GFINDEX3D(GH,i  ,j+1,k-1);
+	  ijmkp = CCTK_GFINDEX3D(GH,i  ,j-1,k+1);
+	  ijmkm = CCTK_GFINDEX3D(GH,i  ,j-1,k-1);
+	  
+	  ipjkp = CCTK_GFINDEX3D(GH,i+1,j  ,k+1);
+	  ipjkm = CCTK_GFINDEX3D(GH,i+1,j  ,k-1);
+	  imjkp = CCTK_GFINDEX3D(GH,i-1,j  ,k+1);
+	  imjkm = CCTK_GFINDEX3D(GH,i-1,j  ,k-1);
+	  
+
+	  ac = -1.0*uxx[ipjk] - 2.0*uxx[ijk] - 1.0*uxx[imjk]
+	    -1.0*uyy[ijpk] - 2.0*uyy[ijk] - 1.0*uyy[ijmk]
+	    -1.0*uzz[ijkp] - 2.0*uzz[ijk] - 1.0*uzz[ijkm];
+	  
+	  
+	  if (Mstorage) 
+	  {
+	    ac += Mlin[ijk];
+	  }
+
+	  ae  = uxx[ipjk]+uxx[ijk];
+	  aw  = uxx[imjk]+uxx[ijk];
+	  an  = uyy[ijpk]+uyy[ijk];
+	  as  = uyy[ijmk]+uyy[ijk];
+	  at  = uzz[ijkp]+uzz[ijk];
+	  ab  = uzz[ijkm]+uzz[ijk];
+	  
+	  ane = uxy[ijpk]+uxy[ipjk];
+	  anw =-uxy[imjk]-uxy[ijpk];
+	  ase =-uxy[ipjk]-uxy[ijmk];
+	  asw = uxy[imjk]+uxy[ijmk];
+	  
+	  ate = uxz[ijkp]+uxz[ipjk];
+	  atw =-uxz[imjk]-uxz[ijkp];
+	  abe =-uxz[ipjk]-uxz[ijkm];
+	  abw = uxz[imjk]+uxz[ijkm];
+	  
+	  atn = uyz[ijpk]+uyz[ijkp];               
+	  ats =-uyz[ijkp]-uyz[ijmk];
+	  abn =-uyz[ijkm]-uyz[ijpk];
+	  absol = uyz[ijkm]+uyz[ijmk];
+	  
+	  residual = ac * var[ijk]
+	    + ae *var[ipjk]  +  aw*var[imjk]
+	    + an *var[ijpk]  +  as*var[ijmk]
+	    + at *var[ijkp]  +  ab*var[ijkm];
+	  
+	  residual = residual 
+	    + ane*var[ipjpk] + anw*var[imjpk]; 
+	  
+	  residual = residual 
+	    + ase*var[ipjmk] + asw*var[imjmk];
+	  
+	  residual = residual    
+	    + ate*var[ipjkp] + atw*var[imjkp] 
+	    + abe*var[ipjkm] + abw*var[imjkm]
+	    + atn*var[ijpkp] + ats*var[ijmkp]
+	    + abn*var[ijpkm] + absol*var[ijmkm];
+	  
+	  if (Nstorage)
+	  {
+	    residual +=Nlin[ijk];
+	  }
     
-	    resnorm  = resnorm + fabs(residual);
-	    var[ijk] = var[ijk] - omega*residual/ac; 
+	  resnorm  = resnorm + fabs(residual);
+	  var[ijk] = var[ijk] - omega*residual/ac; 
 	}
       }
     }
@@ -373,35 +436,46 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
     ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
 				&resnorm, &glob_residual, CCTK_VARIABLE_REAL); 
     if (ierr<0) 
+    {
       CCTK_WARN(1,"Reduction of residual  failed");
+    }
 
     glob_residual = glob_residual /  
       (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]);
 
     /* apply symmetry boundary conditions within loop */    
     if (CartSymVI(GH,FieldIndex)<0) 
+    {
       CCTK_WARN(1,"CartSymVI failed in EllSOR loop");
-  
+    }
+
     /* apply boundary conditions within loop */
     if (CCTK_EQUALS(sor_bound,"robin"))
+    {
       ierr = BndRobinVI(GH, sw, finf, npow,  FieldIndex);
+    }
 
     /* synchronization of grid variable */
     CCTK_SyncGroupWithVarI(GH, FieldIndex);
 
     /* Leave iteration loop if tolerance criterium is met */
     if (glob_residual<tol)
+    {
       break;
+    }
   }
 
   /* Information for the user if the solve did not converge within 
      the given constraints of max.iteration and tolerance */
   if (residual>tol) 
+  {
+    retval = ELL_NOCONVERGENCE;
     CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
 	       "SOR SOLVER DID NOT CONVERGE within given "
 	       "tolerance/max.number of iterations.\n "
 	       "max. iteration %d   residual (tolerance): %g (%g)\n",
 	       maxit, glob_residual, tol );
+  }
 
   if (uxx) free(uxx); 
   if (uyy) free(uyy); 
@@ -410,7 +484,7 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
   if (uxz) free(uxz); 
   if (uyz) free(uyz);
 
-  return;
+  return retval;
 }
 
   
diff --git a/src/sor_flat.c b/src/sor_flat.c
index 67f0463..752277b 100644
--- a/src/sor_flat.c
+++ b/src/sor_flat.c
@@ -1,3 +1,14 @@
+ /*@@
+   @file      sor_flat.c
+   @date      Tue Aug 24 12:50:07 1999
+   @author    Gerd Lanfermann
+   @desc 
+   SOR solver for 3D flat equation
+   @enddesc 
+ @@*/
+
+/*#define DEBUG_ELLIPTIC*/
+
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
@@ -6,21 +17,25 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#include "CactusBase/Boundary/src/Boundary.h"
-#include "CactusBase/CartGrid3D/src/Symmetry.h"
-#include "CactusElliptic/EllBase/src/Ell_DBstructure.h"
+#include "Boundary.h"
+#include "Symmetry.h"
+#include "Ell_DBstructure.h"
+#include "EllBase.h"
  
-void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 		 CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
 
-void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 		 CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
 {
 
   DECLARE_CCTK_PARAMETERS
 
+  int retval = ELL_SUCCESS;
+
   /* The pointer to the data fields */
-  CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
+  CCTK_REAL *Mlin=NULL;
+  CCTK_REAL *Nlin=NULL;   
   CCTK_REAL *var =NULL; 
 
   /* shortcuts for deltas,etc. */
@@ -40,6 +55,7 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   int j,js,je;
   int k,ks,ke,kstep;
   int nxyz;
+  int timelevel;
 
   /* stencil index */
   int ijk;
@@ -55,55 +71,88 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   static int firstcall = 1;
   CCTK_REAL  dx2rec, dy2rec, dz2rec; 
   
-
   /* Get the reduction handle */
   sum_handle = CCTK_ReductionArrayHandle("sum");
   if (sum_handle<0) 
-    CCTK_WARN(1,"Cannot get reduction handle for operation >sum<");
-  
+  {
+    CCTK_WARN(1,"sor_flat_3d: Cannot get reduction handle for sum operation");
+  }
+
   /* IF Robin BCs are set, prepare for a boundary call:
      setup stencil width and get Robin constants (set by the routine
      which is calling the solver interface) */
-  if (CCTK_EQUALS(sor_bound,"robin")) { 
+  if (CCTK_EQUALS(sor_bound,"robin")) 
+  { 
     sw[0]=1; 
     sw[1]=1; 
     sw[2]=1;
     
     ierr = Ell_GetRealKey(&finf, "EllLinConfMetric::Bnd::Robin::inf");
+
     ierr = Ell_GetIntKey (&npow, "EllLinConfMetric::Bnd::Robin::falloff");
   }
 
   /* Only supports absolute tolerance */
   tol   = AbsTol[0];
   if (CCTK_EQUALS(sor_accel,"const"))
+  {
     accel_const = 1;
+  }
   else if (CCTK_EQUALS(sor_accel,"cheb"))
+  {
     accel_cheb  = 1;
+  }
 
   /* Things to do only once! */
-  if (firstcall==1) {
+  if (firstcall==1) 
+  {
     if (CCTK_Equals(elliptic_verbose, "yes"))
+    {
+      if (accel_cheb)
+      {
+	CCTK_INFO("SOR with Chebyshev acceleration with radius of 1");
+      }
+      else if (accel_const)
+      {
+	CCTK_INFO("SOR with hardcoded omega = 1.8");
+      }
+      else
       {
-	if (accel_cheb)
-	  printf("SOR with Chebyshev acceleration with radius of 1\n");
-	else if (accel_const)
-	  printf("SOR with hardcoded omega = 1.8\n");
-	else
-	  printf("SOR with unaccelearted relaxation (omega = 1)\n");
+	CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)");
       }
+    }
     firstcall = 0;
   }
 
   /* Get the data ptr of these GFs, They all have to be
      on the same timelevel; if we have a negative index for M/N, 
      this GF is not set,  there for don't even look for it and flag it  */
-  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex);  
-  if (MIndex>=0)  { 
-    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
+  /* get the current time level */
+  timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);  
+
+  if (MIndex>=0)  
+  { 
+    timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
     Mstorage = 1;
   }
-  if (NIndex>=0) {
-    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex);
+  if (NIndex>=0) 
+  {
+    timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
     Nstorage = 1;
   }
 
@@ -134,26 +183,37 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   kstep = 2;
 
   /* start at 1 for historic (Fortran) reasons */
-  for (sorit=1; sorit<=maxit; sorit++) {
-    
+  for (sorit=1; sorit<=maxit; sorit++) 
+  {
     omega     =  1.0;
     rjacobian =  1.0;
     
     if (accel_cheb) 
+    {
       for (chebit=2;chebit<sorit;chebit++)
+      {
 	omega = 1.0/(1.0 - 0.25*rjacobian*rjacobian*omega);
+      }
+    }
     if (accel_const)
+    {
       omega = 1.8;
+    }
     
     resnorm = 0.0;
     
     ks = (sorit%2)+1;
     if (GH->cctk_lsh[2]==3)
+    {
       ks = 2;
-    
-    for (k=ks;k<ke;k+=kstep) {
-      for (j=js;j<je;j++)     {
-	for (i=is;i<ie;i++)    {
+    }
+
+    for (k=ks;k<ke;k+=kstep) 
+    {
+      for (j=js;j<je;j++)     
+      {
+	for (i=is;i<ie;i++)    
+        {
 	  
 	  ac = ac_orig;
 
@@ -166,7 +226,9 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 	  ijkm  = CCTK_GFINDEX3D(GH,i  ,j  ,k-1);
 	  
 	  if (Mstorage)
+	  {
 	    ac += Mlin[ijk];
+	  }
 
 	  residual = ac * var[ijk]
 	      + ae *var[ipjk]  +  aw*var[imjk]
@@ -174,14 +236,14 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 	      + at *var[ijkp]  +  ab*var[ijkm];
 
 	  if (Nstorage)
+	  {
 	    residual +=Nlin[ijk];
+	  }
 
 	  resnorm  = resnorm + fabs(residual);
 
 	  var[ijk] = var[ijk] - omega*residual/ac; 
 
-	  printf(" %d %d %d  %f \n",i,j,k,var[ijk]);
-
 	}
       }
     }
@@ -190,32 +252,49 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
 				&resnorm, &glob_residual, CCTK_VARIABLE_REAL);
     if (ierr<0) 
+    {
       CCTK_WARN(1,"Reduction of Norm failed");
+    }
 
     glob_residual = glob_residual / 
       (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]);
 
+
+#ifdef DEBUG_ELLIPTIC
+    printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
+#endif
+
     /* apply symmetry boundary conditions within loop */    
-    if (CartSymVI(GH,FieldIndex)<0) 
+    if (CartSymVI(GH,FieldIndex)<0)
+    { 
       CCTK_WARN(1,"CartSymVI failed in EllSOR loop");
+      break;
+    }
 
     /* apply boundary conditions within loop */
     if (CCTK_EQUALS(sor_bound,"robin"))
+    {
       ierr = BndRobinVI(GH, sw, finf, npow,  FieldIndex);
+    }
 
     /* synchronization of grid variable */
     CCTK_SyncGroupWithVarI(GH, FieldIndex);
 
     /* Leave iteration loop if tolerance criterium is met */
     if (glob_residual<tol)
+    {
       break;
+    }
     
   }
   
   if (glob_residual>tol) 
-    CCTK_WARN(2,"SOR SOLVER DID NOT CONVERGE");
+  {
+    CCTK_WARN(1,"sor_flat: SOR Solver did not converge");
+    retval = ELL_NOCONVERGENCE;
+  }
 
-  return;
+  return retval;
 }
  
   
diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c
index 4f2a8b0..062a846 100644
--- a/src/sor_wrapper.c
+++ b/src/sor_wrapper.c
@@ -22,19 +22,21 @@
 #include "cctk_Parameters.h"
 #include "cctk_FortranString.h"
 
-void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
-		 CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
-void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
-		       int FieldIndex, int MIndex, int NIndex,
-		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
-void sor_confmetric(cGH *GH, 
-		    int *MetricPsiI, 
-		    int FieldIndex, 
-		    int MIndex, 
-		    int NIndex, 
-		    CCTK_REAL *AbsTol,
-		    CCTK_REAL *RelTol);
-void sor_flat(cGH *GH, 
+#include "CactusElliptic/EllBase/src/EllBase.h"
+
+int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+		CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
+		      int FieldIndex, int MIndex, int NIndex,
+		      CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+int sor_confmetric(cGH *GH, 
+		   int *MetricPsiI, 
+		   int FieldIndex, 
+		   int MIndex, 
+		   int NIndex, 
+		   CCTK_REAL *AbsTol,
+		   CCTK_REAL *RelTol);
+int sor_flat(cGH *GH, 
 	     int FieldIndex, 
 	     int MIndex, 
 	     int NIndex, 
@@ -65,32 +67,36 @@ void sor_flat(cGH *GH,
 
 @@*/
 
-void sor_confmetric(cGH *GH, 
-		    int *MetricPsiI, 
-		    int FieldIndex, 
-		    int MIndex, 
-		    int NIndex, 
-		    CCTK_REAL *AbsTol,
-		    CCTK_REAL *RelTol) 
+int sor_confmetric(cGH *GH, 
+		   int *MetricPsiI, 
+		   int FieldIndex, 
+		   int MIndex, 
+		   int NIndex, 
+		   CCTK_REAL *AbsTol,
+		   CCTK_REAL *RelTol) 
 {
+  int retval = ELL_NOSOLVER;
   
   switch (CCTK_GroupDimFromVarI(FieldIndex))
   {
     case 1:  
-      CCTK_WARN(0,"No 1D SOR solver implemented");
+      CCTK_WARN(0,"sor_confmetric: No 1D SOR solver implemented");
       break;
     case 2:
-      CCTK_WARN(0,"No 2D SOR solver implemented");
+      CCTK_WARN(0,"sor_confmetric: No 2D SOR solver implemented");
       break;
     case 3:
-      sor_confmetric_3d(GH, MetricPsiI, 1, FieldIndex, MIndex, NIndex,
-			AbsTol, RelTol);
+      retval = sor_confmetric_3d(GH, MetricPsiI, 1, 
+				 FieldIndex, MIndex, NIndex,
+				 AbsTol, RelTol);
       break;
     default:
-      CCTK_WARN(1,"Requested SOR solver for dimension equal "
+      CCTK_WARN(1,"sor_confmetric: Requested SOR solver for dimension equal "
 		  "zero or gt. three not implemented!");
       break;
   }
+
+  return retval;
 }
  
 /*@@
@@ -105,9 +111,9 @@ void sor_confmetric(cGH *GH,
      This wrapper is registered and if it is being called with 
      a n-dim. grid function, it goes of and picks the correct solver.
 
-     We pass in the arguments that are neccessary for this class of elliptic eq. 
-     this solver is intended to solve. See ./CactusElliptic/EllBase/src/ for the
-     classes of elliptic eq. 
+     We pass in the arguments that are necessary for this class of 
+     elliptic equations this solver is intended to solve. 
+     See ./CactusElliptic/EllBase/src/ for the classes of elliptic eq. 
    @enddesc 
    @calls     
    @calledby   
@@ -117,30 +123,32 @@ void sor_confmetric(cGH *GH,
 
 @@*/
 
-void sor_flat(cGH *GH, 
+int sor_flat(cGH *GH, 
 	     int FieldIndex, 
 	     int MIndex, 
 	     int NIndex, 
 	     CCTK_REAL *AbsTol, 
 	     CCTK_REAL *RelTol) 
 {
+  int retval = ELL_NOSOLVER;
   switch (CCTK_GroupDimFromVarI(FieldIndex))
   {
     case 1:  
-      CCTK_WARN(1,"No 1D SOR solver implemented");
+      CCTK_WARN(1,"sor_flat: No 1D SOR solver implemented");
       break;
     case 2:
-      CCTK_WARN(1,"No 2D SOR solver implemented");
+      CCTK_WARN(1,"sor_flat: No 2D SOR solver implemented");
       break;
     case 3:
-      printf(" wrapper: %d %d %d \n",FieldIndex, MIndex, NIndex);
-      sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol);
+      retval = sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol);
       break;
   default:
-      CCTK_WARN(1,"Requested SOR solver for dimension equal" 
+      CCTK_WARN(1,"sor_flat: Requested SOR solver for dimension equal" 
 		  "zero or gt. three not implemented!");
       break;
   }
 
+  return retval;
+
 }
  
-- 
cgit v1.2.3