New tidied version of SOR, Linear solver uses grid sweeping which gives

consistent results for multiprocessors


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

7 files changed, 1001 insertions, 860 deletions

diff --git a/src/ConfMetric.c b/src/ConfMetric.c
new file mode 100644
index 0000000..be9b0d2
--- /dev/null
+++ b/src/ConfMetric.c
@@ -0,0 +1,546 @@
+ /*@@
+   @file      ConfMetric.c
+   @date      Tue Sep 26 11:29:18 2000
+   @author    Gerd Lanfermann
+   @desc 
+     Provides sor solver engine routines
+   @enddesc 
+ @@*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "Boundary.h"
+#include "Symmetry.h"
+#include "Ell_DBstructure.h"
+#include "EllBase.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusElliptic_EllSOR_ConfMetric_c)
+
+#define SQR(a) ((a)*(a))
+
+static int firstcall = 1;
+
+int SORConfMetric3D(cGH *GH, int *MetricPsiI, int conformal,
+		    int FieldIndex, int MIndex, int NIndex,
+		    CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+
+ /*@@
+   @routine    SORConfMetric3D
+   @date       Tue Sep 26 11:28:08 2000
+   @author     Joan Masso, Paul Walker, Gerd Lanfermann
+   @desc 
+   This is a standalone sor solver engine, 
+   called by the wrapper functions
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int SORConfMetric3D(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;
+  int origin_sign;
+
+  /* The pointer to the metric fields */
+  CCTK_REAL *gxx =NULL, *gxy =NULL; 
+  CCTK_REAL *gxz =NULL, *gyy =NULL; 
+  CCTK_REAL *gyz =NULL, *gzz =NULL; 
+  CCTK_REAL *psi =NULL;
+  CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
+  CCTK_REAL *var =NULL; 
+
+  /* The inverse metric, allocated here */
+  CCTK_REAL *uxx=NULL, *uyy=NULL, 
+            *uzz=NULL, *uxz=NULL,
+            *uxy=NULL, *uyz=NULL;
+
+  /* shortcuts for metric, psi, deltas, etc. */
+  CCTK_REAL pm4, p12, det;
+
+  CCTK_REAL dx,dy,dz;
+  CCTK_REAL dxdx, dydy, dzdz, 
+            dxdy, dxdz, dydz;
+
+  /* Some physical variables */
+  CCTK_REAL omega, resnorm=0.0, residual=0.0;
+  CCTK_REAL glob_residual=0.0, rjacobian=0.0; 
+  CCTK_REAL finf;
+  CCTK_INT npow;
+  CCTK_REAL tol;
+
+  /* Iteration / stepping  variables */
+  int sorit, maxit; 
+  int i,j,k;
+  int nxyz;
+  
+  /* 19 point stencil index */
+  int ijk;
+  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
+  int ipjpk, ipjmk, imjpk, imjmk;
+  int ipjkp, ipjkm, imjkp, imjkm;
+  int ijpkp, ijpkm, ijmkp, ijmkm;
+ 
+  /* 19 point stencil coefficients */
+  CCTK_REAL ac;
+  CCTK_REAL ae,aw,an,as,at,ab;
+  CCTK_REAL ane, anw, ase, asw, ate, atw, abe, abw;
+  CCTK_REAL atn, ats, abn, absol;
+
+  /* Miscellaneous */
+  int sum_handle=-1;
+  int sw[3], ierr;
+  int Mstorage=0, Nstorage=0;
+  size_t varsize;
+  CCTK_REAL detrec_pm4, sqrtdet;
+  char *robin=NULL;
+  char *sor_accel=NULL;
+
+
+  /* Get the reduction handle */
+  sum_handle = CCTK_ReductionArrayHandle("sum");
+  if (sum_handle<0) 
+  {
+    CCTK_WARN(1,"SORConfMetric3D: "
+	      "Cannot get handle for sum reduction");
+    retval = ELL_FAILURE;
+  }
+
+  /* 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 (Ell_GetStrKey(&robin, "EllLinConfMetric::Bnd::Robin")< 0)
+  {
+    CCTK_WARN(1,"SORConfMetric3D: Robin keys not set");
+  }
+  
+  if (CCTK_EQUALS(robin,"yes")) 
+  { 
+    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");
+  }
+
+  /* Get the maximum number of iterations allowed. */
+  if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET)
+  {
+    CCTK_WARN(1,"SORConfMetric3D: Ell::SORmaxit not set. "
+	      "Maximum allowed iterations being set to 100.");
+    maxit = 100;
+  }
+
+  /* Only supports absolute tolerance */
+  tol   = AbsTol[0];
+
+  /* Get the type of acceleration to use. */
+  if (Ell_GetStrKey(&sor_accel, "Ell::SORaccel") == ELLGET_NOTSET)
+  {
+    const char tmpstr[6] = "const";
+    CCTK_WARN(1,"SORConfMetric3D: Ell::SORaccel not set. "
+	      "Omega being set to a constant value of 1.8.");
+    sor_accel = strdup(tmpstr);
+  }
+
+  /* Things to do only once! */
+  if (firstcall==1) 
+  {
+    if (CCTK_Equals(elliptic_verbose, "yes"))
+    {
+      if (CCTK_Equals(sor_accel, "cheb"))
+      {
+        CCTK_INFO("SOR with Chebyshev acceleration");
+      }
+      else if (CCTK_Equals(sor_accel, "const"))
+      {
+        CCTK_INFO("SOR with hardcoded omega = 1.8");
+      }
+      else if (CCTK_Equals(sor_accel, "none"))
+      {
+        CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)");
+      }
+      else
+      {
+        CCTK_WARN(0, "SORConfMetric3D: sor_accel not set");
+      }
+    }
+    firstcall = 0;
+  }
+
+  /* 
+     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   
+  */
+
+  timelevel = 0;
+  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);
+
+  timelevel = 0;
+  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)
+  {
+    timelevel = 0;
+    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)  
+  { 
+    timelevel = 0;
+    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
+    if (Mlin)
+    {
+      Mstorage = 1;
+    }
+    else
+    {
+      CCTK_WARN(0, "SORConfMetric3D: Unable to get pointer to M.");
+    }
+  }
+  if (NIndex>=0) 
+  {
+    timelevel = 0;
+    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
+    if (Nlin)
+    {
+      Nstorage = 1;
+    }
+    else
+    {
+      CCTK_WARN(0, "SORConfMetric3D: Unable to get pointer to N.");
+    }
+  }
+
+  /* Shortcuts */
+  dx   = GH->cctk_delta_space[0];
+  dy   = GH->cctk_delta_space[1];
+  dz   = GH->cctk_delta_space[2];
+  nxyz = GH->cctk_lsh[0]*GH->cctk_lsh[1]*GH->cctk_lsh[2];
+
+  /* Allocate the inverse metric */
+  varsize = (size_t)CCTK_VarTypeSize(CCTK_VarTypeI(FieldIndex))*nxyz;
+  uxx = (CCTK_REAL*) malloc(varsize);
+  uxy = (CCTK_REAL*) malloc(varsize);
+  uxz = (CCTK_REAL*) malloc(varsize);
+  uyy = (CCTK_REAL*) malloc(varsize);
+  uyz = (CCTK_REAL*) malloc(varsize);
+  uzz = (CCTK_REAL*) malloc(varsize);
+
+  if (!uxx || !uxy || !uxz || !uyy || !uyz || !uzz)
+  {
+    CCTK_WARN(0,"SORConfMetric3D: Memory allocation failed for upper metric");
+  }
+  
+  /* calculate the differential coefficient */
+  dxdx = 0.5/(dx*dx);
+  dydy = 0.5/(dy*dy);
+  dzdz = 0.5/(dz*dz);
+  dxdy = 0.25/(dx*dy);
+  dxdz = 0.25/(dx*dz);
+  dydz = 0.25/(dy*dz);
+
+  /* Calculate the inverse metric */
+  for (ijk=0; ijk<nxyz; ijk++) 
+  {
+    det = -(SQR(gxz[ijk])*gyy[ijk]) + 
+      2*gxy[ijk]*gxz[ijk]*gyz[ijk] - 
+      gxx[ijk]*SQR(gyz[ijk])  -
+      SQR(gxy[ijk])*gzz[ijk] + 
+      gxx[ijk]*gyy[ijk]*gzz[ijk];
+    
+    if (conformal) 
+    {
+      pm4 = 1.0/pow(psi[ijk],4.0);
+      p12 = pow(psi[ijk],12.0);
+    } 
+    else 
+    {
+      pm4 = 1.0;
+      p12 = 1.0;
+    }
+
+    /* try to avoid constly division */
+    detrec_pm4 = 1.0/det*pm4;
+    sqrtdet    = sqrt(det);
+
+    uxx[ijk]=(-SQR(gyz[ijk])     + gyy[ijk]*gzz[ijk])*detrec_pm4;
+    uxy[ijk]=( gxz[ijk]*gyz[ijk] - gxy[ijk]*gzz[ijk])*detrec_pm4;
+    uxz[ijk]=(-gxz[ijk]*gyy[ijk] + gxy[ijk]*gyz[ijk])*detrec_pm4;
+    uyy[ijk]=(-SQR(gxz[ijk])     + gxx[ijk]*gzz[ijk])*detrec_pm4;
+    uyz[ijk]=( gxy[ijk]*gxz[ijk] - gxx[ijk]*gyz[ijk])*detrec_pm4;
+    uzz[ijk]=(-SQR(gxy[ijk])     + gxx[ijk]*gyy[ijk])*detrec_pm4;
+ 	
+    det    = det*p12;
+    sqrtdet= sqrt(det);
+   
+    /* Rescaling for the uppermetric solver */
+    if (Mstorage) 
+    {
+      Mlin[ijk] = Mlin[ijk]*sqrt(det);
+    }
+    if (Nstorage)
+    {
+      Nlin[ijk] = Nlin[ijk]*sqrt(det);
+    }
+    
+    uxx[ijk]=uxx[ijk]*dxdx*sqrtdet;
+    uyy[ijk]=uyy[ijk]*dydy*sqrtdet;
+    uzz[ijk]=uzz[ijk]*dzdz*sqrtdet;
+    uxy[ijk]=uxy[ijk]*dxdy*sqrtdet;
+    uxz[ijk]=uxz[ijk]*dxdz*sqrtdet;
+    uyz[ijk]=uyz[ijk]*dydz*sqrtdet;
+
+  }	
+
+  /* Initialize omega. */
+  /* TODO: Make it so that the value of the constant omega can be set. */
+  if (CCTK_Equals(sor_accel, "const"))
+  {
+    omega = 1.8;
+  }
+  else
+  {
+    omega = 1.;
+  }
+
+  /* Set the spectral radius of the Jacobi iteration. */
+  rjacobian =  0.999;
+
+  /* Compute whether the origin of this processor's 
+     sub grid is "red" or "black". */
+
+  origin_sign = (GH->cctk_lbnd[0] + GH->cctk_lbnd[1] + GH->cctk_lbnd[2])%2;
+
+  /* start at 1 for historic (Fortran) reasons */
+  for (sorit=1; sorit<=maxit; sorit++)
+  {
+    resnorm = 0.0;
+
+    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++)    
+        {
+	  if ((origin_sign+i+j+k)%2 == sorit%2)
+          {
+            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; 
+
+          }
+	}
+      }
+    }
+
+    /* reduction operation on processor-local residual values */
+    if (CCTK_ReduceLocScalar(GH, -1, sum_handle,
+			     &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) 
+    {
+      CCTK_WARN(1,"SORConfMetric3D: 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,"SORConfMetric3D: CartSymVI failed in EllSOR loop");
+    }
+
+    /* apply Robin boundary conditions within loop */
+    if (CCTK_EQUALS(robin,"yes")) 
+    {
+      if (BndRobinVI(GH, sw, finf, npow,  FieldIndex)<0)
+      {
+        CCTK_WARN(1, "SORConfMetric3D: BndRobinVI failed in EllSOR loop");
+        break;
+      }
+    }
+
+    /* synchronization of grid variable */
+    ierr = CCTK_SyncGroupWithVarI(GH, FieldIndex);
+    if (ierr < 0)
+    {
+      CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ConfMetric3D: Synchronization of %s failed (Error: %d)",
+		 CCTK_VarName(FieldIndex),ierr);
+    }
+
+    /* Leave iteration loop if tolerance criterium is met */
+    if (glob_residual<tol)
+    {
+      break;
+    }
+
+    /* Update omega if using Chebyshev acceleration. */
+    if (CCTK_Equals(sor_accel, "cheb"))
+    {
+      if (sorit == 1)
+      {
+          omega = 1. / (1. - 0.5 * rjacobian * rjacobian);
+      }
+      else
+      {
+         omega = 1. / (1. - 0.25 * rjacobian * rjacobian * omega);
+      }
+    }
+
+  }
+
+  if (elliptic_verbose)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Required SOR tolerence was set at %f",tol);
+    CCTK_VInfo(CCTK_THORNSTRING,"Achieved SOR residual was %f",glob_residual);
+    CCTK_VInfo(CCTK_THORNSTRING,
+	       "Number of iterations was %d (max: %d)",sorit,maxit);
+  }
+
+  if (glob_residual>tol) 
+  {
+    CCTK_WARN(1, "SORConfMetric3D: SOR Solver did not converge");
+    retval = ELL_NOCONVERGENCE;
+  }
+
+  if (uxx) 
+  {
+    free(uxx); 
+  }
+  if (uyy) 
+  {
+    free(uyy); 
+  }
+  if (uzz)
+  { 
+    free(uzz);
+  }
+  if (uxy) 
+  {
+    free(uxy); 
+  }
+  if (uxz) 
+  {
+    free(uxz); 
+  }
+  if (uyz) 
+  {
+    free(uyz);
+  }
+  if (robin) 
+  {
+    free(robin);
+  }
+  if (sor_accel) 
+  {
+    free(sor_accel);
+  }
+
+  return retval;
+}
diff --git a/src/Flat.c b/src/Flat.c
new file mode 100644
index 0000000..6fcbb0e
--- /dev/null
+++ b/src/Flat.c
@@ -0,0 +1,373 @@
+ /*@@
+   @file      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>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "Boundary.h"
+#include "Symmetry.h"
+#include "Ell_DBstructure.h"
+#include "EllBase.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusElliptic_EllSOR_Flat_c)
+ 
+int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+	      CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+
+int SORFlat3D(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;
+  CCTK_REAL *Nlin=NULL;   
+  CCTK_REAL *var =NULL; 
+
+  /* shortcuts for deltas,etc. */
+  CCTK_REAL dx,dy,dz; 
+
+  /* Some physical variables */
+  CCTK_REAL omega, resnorm, residual, glob_residual, rjacobian; 
+  CCTK_REAL finf;
+  CCTK_INT npow;
+  CCTK_REAL tol;
+
+  /* Iteration and stepping  variables */
+  int sorit, maxit;
+  int i, j, k;
+  int origin_sign;
+
+  /* stencil index */
+  int ijk;
+  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
+
+  /* Coefficients for the stencil...  */
+  CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab;
+ 
+  /* Miscellaneous */
+  int sum_handle;
+  int sw[3];
+  int Mstorage=0, Nstorage=0;
+  static int firstcall = 1;
+  CCTK_REAL  dx2rec, dy2rec, dz2rec; 
+  char *robin = NULL;
+  char *sor_accel = NULL;
+
+  /* Boundary conditions */
+  int use_robin = 0;
+
+  /* Avoid compiler warnings */
+  RelTol = RelTol;
+
+  /* Get the reduction handle */
+  sum_handle = CCTK_ReductionArrayHandle("sum");
+  if (sum_handle<0) 
+  {
+    CCTK_WARN(1,"SORFlat3D: 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 (Ell_GetStrKey(&robin,"EllLinFlat::Bnd::Robin") < 0)
+  {
+    CCTK_WARN(2,"SORFlat3D: Robin keys not set");
+  }
+  else
+  {
+    if (CCTK_Equals(robin,"yes"))
+    {
+      use_robin = 1;
+
+      sw[0]=1; 
+      sw[1]=1; 
+      sw[2]=1;
+
+      if (Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf") == 
+	  ELLGET_NOTSET)
+      {
+	CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::inf not set, "
+		  "setting to 1");
+	finf = 1;
+      }
+      if (Ell_GetIntKey(&npow, "EllLinFlat::Bnd::Robin::falloff") 
+	  == ELLGET_NOTSET)
+      {
+	CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::falloff not set, "
+		  "setting to 1");
+	npow = 1;
+      }
+    }
+  }
+
+  /* Get the maximum number of iterations allowed. */
+  if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET)
+  {
+    CCTK_WARN(1,"SORFlat3D: Ell::SORmaxit not set. "
+	      "Maximum allowed iterations being set to 100.");
+    maxit = 100;
+  }
+  
+  /* Only supports absolute tolerance */
+  tol   = AbsTol[0];
+
+  /* Get the type of acceleration to use. */
+  if (Ell_GetStrKey(&sor_accel, "Ell::SORaccel") == ELLGET_NOTSET)
+  {
+    const char tmpstr[6] = "const";
+    CCTK_WARN(1, "SORFlat3D: Ell::SORaccel not set. "
+	      "Omega being set to a constant value of 1.8.");
+    sor_accel = strdup(tmpstr);
+  }
+
+  /* Things to do only once! */
+  /* TODO: Need to handle this differently, since it may be called
+           for different reasons by the same code. */
+  if (firstcall==1) 
+  {
+    if (CCTK_Equals(elliptic_verbose, "yes"))
+    {
+      if (CCTK_Equals(sor_accel, "cheb"))
+      {
+	CCTK_INFO("SOR with Chebyshev acceleration");
+      }
+      else if (CCTK_Equals(sor_accel, "const"))
+      {
+	CCTK_INFO("SOR with hardcoded omega = 1.8");
+      }
+      else if (CCTK_Equals(sor_accel, "none"))
+      {
+	CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)");
+      }
+      else
+      {
+         CCTK_WARN(0, "SORFlat3D: sor_accel not set");
+      }
+    }
+    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 (var == NULL)
+  {
+    CCTK_WARN(0, "SORFlat3D: Unable to get pointer to GF variable");
+  }
+
+  if (MIndex>=0)  
+  { 
+    Mlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, MIndex);
+    if (Mlin)
+    {
+      Mstorage = 1;
+    }
+    else
+    {
+      CCTK_WARN(0, "SORFlat3D: Unable to get pointer to M");
+    }
+  }
+
+  if (NIndex>=0) 
+  {
+    Nlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, NIndex);
+    if (Nlin)
+    {
+      Nstorage = 1;
+    }
+    else
+    {
+      CCTK_WARN(0, "SORFlat3D: Unable to get pointer to N");
+    }
+  }
+
+  /* Shortcuts */
+  dx   = GH->cctk_delta_space[0];
+  dy   = GH->cctk_delta_space[1];
+  dz   = GH->cctk_delta_space[2];
+  
+  dx2rec = 1.0/(dx*dx);
+  dy2rec = 1.0/(dy*dy);
+  dz2rec = 1.0/(dz*dz);
+
+  ae = dx2rec;
+  aw = dx2rec;
+  an = dy2rec;
+  as = dy2rec;
+  at = dz2rec;
+  ab = dz2rec;
+
+  ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec;
+
+  /* Initialize omega. */
+  /* TODO: Make it so that the value of the constant omega can be set. */
+  if (CCTK_Equals(sor_accel, "const"))
+  {
+    omega = 1.8;
+  }
+  else
+  {
+    omega = 1.;
+  }
+
+  /* Set the spectral radius of the Jacobi iteration. */
+  /* TODO: I think this can be easily computed for flat metrics? */
+  rjacobian =  0.999;
+
+  /* Compute whether the origin of this processor's sub grid is "red" or
+     "black". */
+  origin_sign = (GH->cctk_lbnd[0] + GH->cctk_lbnd[1] + GH->cctk_lbnd[2])%2;
+
+  /* start at 1 for historic (Fortran) reasons */
+  for (sorit=1; sorit<=maxit; sorit++) 
+  {
+    resnorm = 0.;
+    
+    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++)    
+        {
+          if ((origin_sign + i + j + k)%2 == sorit%2)
+          {
+	    ac = ac_orig;
+
+	    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);
+	  
+	    if (Mstorage)
+	    {
+	      ac += Mlin[ijk];
+	    }
+
+	    residual = ac * var[ijk]
+	        + ae*var[ipjk] + aw*var[imjk]
+	        + an*var[ijpk] + as*var[ijmk]
+	        + at*var[ijkp] + ab*var[ijkm];
+
+	    if (Nstorage)
+	    {
+	      residual += Nlin[ijk];
+	    }
+
+	    resnorm += fabs(residual);
+
+	    var[ijk] -= omega*residual/ac; 
+          }
+	}
+      }
+    }
+
+    /* reduction operation on processor-local residual values */
+
+    if (CCTK_ReduceLocScalar(GH, -1, sum_handle,
+			     &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) 
+    {
+      CCTK_WARN(1,"SORFlat3D: Reduction of Norm failed");
+    }
+
+#ifdef DEBUG_ELLIPTIC
+    printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
+#endif
+
+    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)
+    { 
+      CCTK_WARN(1,"SORFlat3D: CartSymVI failed in EllSOR loop");
+      break;
+    }
+
+    /* apply Robin boundary conditions within loop */
+    if (use_robin)
+    {
+      if (BndRobinVI(GH, sw, finf, npow,  FieldIndex)<0)
+      { 
+        CCTK_WARN(1,"SORFlat3D: BndRobinVI failed in EllSOR loop");
+        break;
+      }
+    }
+
+    /* synchronization of grid variable */
+    CCTK_SyncGroupWithVarI(GH, FieldIndex);
+
+    /* Leave iteration loop if tolerance criterium is met */
+    if (glob_residual<tol)
+    {
+      break;
+    }
+
+    /* Update omega if using Chebyshev acceleration. */
+    if (CCTK_Equals(sor_accel, "cheb"))
+    {
+      if (sorit == 1)
+      {
+        omega = 1. / (1. - 0.5 * rjacobian * rjacobian);
+      }
+      else
+      {
+         omega = 1. / (1. - 0.25 * rjacobian * rjacobian * omega);
+      }
+    }
+    
+  }
+
+  if (elliptic_verbose)
+  {
+    CCTK_VInfo("EllSOR","Required SOR tolerence was set at %f",tol);
+    CCTK_VInfo("EllSOR","Achieved SOR residual was %f",glob_residual);
+    CCTK_VInfo("EllSOR","Number of iterations was %d (max: %d)",sorit,maxit);
+  }
+  
+  if (glob_residual>tol) 
+  {
+    CCTK_WARN(1,"SORFlat3D: SOR Solver did not converge");
+    retval = ELL_NOCONVERGENCE;
+  }
+
+  if (use_robin) 
+  {
+    free(robin);
+  }
+  if (sor_accel)
+  { 
+    free(sor_accel);
+  }
+
+  return retval;
+}
+ 
+  
diff --git a/src/Startup.c b/src/Startup.c
index dfe7d5f..19a75bd 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -15,21 +15,21 @@ static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusElliptic_EllSOR_Startup_c)
 
 void EllSOR_Register(CCTK_ARGUMENTS);
-void sor_confmetric(cGH *GH, 
+void SORConfMetric(cGH *GH, 
 	            int *MetricPsiI, 
 	            int FieldI, 
 	            int MI,
                     int NI, 
 	            CCTK_REAL *AbsTol, 
 	            CCTK_REAL *RelTol);
-void sor_flat(cGH *GH, 
-	      int FieldI,
-	      int MI, 
-	      int NI, 
-	      CCTK_REAL *AbsTol, 
-	      CCTK_REAL *RelTol);
-
-/* routine registers the SOR solver "sor_confmetric" under the name "sor" 
+void SORFlat(cGH *GH, 
+	     int FieldI,
+	     int MI, 
+	     int NI, 
+	     CCTK_REAL *AbsTol, 
+	     CCTK_REAL *RelTol);
+
+/* routine registers the SOR solver "SORConfMetric" under the name "sor" 
    with the Elliptic Class "LinConfMetric".
 */
 
@@ -39,20 +39,20 @@ void EllSOR_Register(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS 
   DECLARE_CCTK_PARAMETERS 
 
-
   int err;
 
   /* Register the solver with the elliptic classes */
 
-  err = Ell_RegisterSolver(sor_confmetric,"sor","Ell_LinConfMetric");
+  err = Ell_RegisterSolver(SORConfMetric,"sor","Ell_LinConfMetric");
   if (!err==ELL_SUCCESS)
   {
-    CCTK_WARN(0,"Failed to register sor for Ell_LinConfMetric");
+    CCTK_WARN(0,"EllSOR_Register: Failed to register sor for Ell_LinConfMetric");
   }
-  err = Ell_RegisterSolver(sor_flat,"sor","Ell_LinFlat");
+
+  err = Ell_RegisterSolver(SORFlat,"sor","Ell_LinFlat");
   if (!err==ELL_SUCCESS)
   {
-    CCTK_WARN(0,"Failed to register sor for Ell_LinFlat");
+    CCTK_WARN(0,"EllSOR_Register: Failed to register sor for Ell_LinFlat");
   }
 
   /* These "keys" have to be same in other elliptic solver and in 
@@ -62,6 +62,22 @@ void EllSOR_Register(CCTK_ARGUMENTS)
   err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Robin");
   err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Const");
 
-  
+  /* Create a key to be associated with the maximum number of iterations
+     allowed. */
+  err = Ell_CreateKey(CCTK_VARIABLE_INT, "Ell::SORmaxit");
+  if (err != ELL_SUCCESS)
+  {
+    CCTK_WARN(0, "EllSOR_Register: "
+	      "Failed to create integer key Ell::SORmaxit");
+  } 
+
+  /* Create a key to be associated with the type of acceleration to be
+     used. */
+  err = Ell_CreateKey(CCTK_VARIABLE_STRING, "Ell::SORaccel");
+  if (err != ELL_SUCCESS)
+  {
+    CCTK_WARN(0, "EllSOR_Register: "
+	      "Failed to create integer key Ell::SORaccel");
+  } 
 
 }
diff --git a/src/sor_wrapper.c b/src/Wrapper.c
index 3df20c5..b3ca1c5 100644
--- a/src/sor_wrapper.c
+++ b/src/Wrapper.c
@@ -1,11 +1,11 @@
  /*@@
-   @file      jacobi_wrapper.c
+   @file      Wrapper.c
    @date      Tue Aug 24 12:50:07 1999
    @author    Gerd Lanfermann
    @desc 
    The C wrapper, which calles the core Fortran routine, which 
    performs the actual solve.
-   We cannot derive the pointers to the GF data from the indeces in 
+   We cannot derive the pointers to the GF data from the indices in 
    Fortran. So we do this here in C and then pass the everything 
    over to the Fortran routine.
 
@@ -22,33 +22,33 @@
 #include "cctk_Parameters.h"
 #include "cctk_FortranString.h"
 
-#include "CactusElliptic/EllBase/src/EllBase.h"
+#include "EllBase.h"
 
 static const char *rcsid = "$Header$";
 
-CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_wrapper_c)
-
-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, 
-	     CCTK_REAL *AbsTol, 
-	     CCTK_REAL *RelTol);
+CCTK_FILEVERSION(CactusElliptic_EllSOR_Wrapper_c)
+
+int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+	      CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+int SORConfMetric3D(cGH *GH, int *MetricPsiI, int conformal,
+		    int FieldIndex, int MIndex, int NIndex,
+		    CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+int SORConfMetric(cGH *GH, 
+		  int *MetricPsiI, 
+		  int FieldIndex, 
+		  int MIndex, 
+		  int NIndex, 
+		  CCTK_REAL *AbsTol,
+		  CCTK_REAL *RelTol);
+int SORFlat(cGH *GH, 
+	    int FieldIndex, 
+	    int MIndex, 
+	    int NIndex, 
+	    CCTK_REAL *AbsTol, 
+	    CCTK_REAL *RelTol);
 
 /*@@
-   @routine    sor_confmetric
+   @routine    SORConfMetric
    @date       Tue Sep 26 11:31:42 2000
    @author     Gerd Lanfermann
    @desc 
@@ -71,7 +71,7 @@ int sor_flat(cGH *GH,
 
 @@*/
 
-int sor_confmetric(cGH *GH, 
+int SORConfMetric(cGH *GH, 
 		   int *MetricPsiI, 
 		   int FieldIndex, 
 		   int MIndex, 
@@ -80,23 +80,22 @@ int sor_confmetric(cGH *GH,
 		   CCTK_REAL *RelTol) 
 {
   int retval = ELL_NOSOLVER;
-  
+
   switch (CCTK_GroupDimFromVarI(FieldIndex))
   {
     case 1:  
-      CCTK_WARN(0,"sor_confmetric: No 1D SOR solver implemented");
+      CCTK_WARN(0,"SORConfMetric: No 1D SOR solver implemented");
       break;
     case 2:
-      CCTK_WARN(0,"sor_confmetric: No 2D SOR solver implemented");
+      CCTK_WARN(0,"SORConfMetric: No 2D SOR solver implemented");
       break;
     case 3:
-      retval = sor_confmetric_3d(GH, MetricPsiI, 1, 
-				 FieldIndex, MIndex, NIndex,
-				 AbsTol, RelTol);
+      retval = SORConfMetric3D(GH, MetricPsiI, 1, 
+			       FieldIndex, MIndex, NIndex,
+			       AbsTol, RelTol);
       break;
     default:
-      CCTK_WARN(1,"sor_confmetric: Requested SOR solver for dimension equal "
-		  "zero or gt. three not implemented!");
+      CCTK_WARN(1,"SORConfMetric: Solver only implemented for 3D");
       break;
   }
 
@@ -104,20 +103,20 @@ int sor_confmetric(cGH *GH,
 }
  
 /*@@
-   @routine    sor_flat
+   @routine    SORFlat
    @date       Tue Sep 26 11:31:42 2000
    @author     Gerd Lanfermann
    @desc 
-     elliptic solver wrapper which provides a interface to the 
+     Elliptic solver wrapper which provides a interface to the 
      different n-dimensional SOR solvers (flat case), 
      of which only 3D is coded currently.
 
      This wrapper is registered and if it is being called with 
-     a n-dim. grid function, it goes of and picks the correct solver.
+     a n-dimensional grid function, it then picks the correct solver.
 
      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. 
+     See ./CactusElliptic/EllBase/src/ for the classes of elliptic equations.
    @enddesc 
    @calls     
    @calledby   
@@ -127,28 +126,30 @@ int sor_confmetric(cGH *GH,
 
 @@*/
 
-int sor_flat(cGH *GH, 
-	     int FieldIndex, 
-	     int MIndex, 
-	     int NIndex, 
-	     CCTK_REAL *AbsTol, 
-	     CCTK_REAL *RelTol) 
+int SORFlat(cGH *GH, 
+	    int FieldIndex, 
+	    int MIndex, 
+	    int NIndex, 
+	    CCTK_REAL *AbsTol, 
+	    CCTK_REAL *RelTol) 
 {
-  int retval = ELL_NOSOLVER;
+  int retval;
+
+  retval = ELL_NOSOLVER;
+
   switch (CCTK_GroupDimFromVarI(FieldIndex))
   {
     case 1:  
-      CCTK_WARN(1,"sor_flat: No 1D SOR solver implemented");
+      CCTK_WARN(1,"SORFlat: No 1D SOR solver implemented");
       break;
     case 2:
-      CCTK_WARN(1,"sor_flat: No 2D SOR solver implemented");
+      CCTK_WARN(1,"SORFlat: No 2D SOR solver implemented");
       break;
     case 3:
-      retval = sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol);
+      retval = SORFlat3D(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol);
       break;
   default:
-      CCTK_WARN(1,"sor_flat: Requested SOR solver for dimension equal" 
-		  "zero or gt. three not implemented!");
+      CCTK_WARN(1,"SORFlat: Solver only implemented for 3D");
       break;
   }
 
diff --git a/src/make.code.defn b/src/make.code.defn
index 4e86518..117e3e5 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Startup.c sor_wrapper.c sor_confmetric.c sor_flat.c
+SRCS = Startup.c Wrapper.c ConfMetric.c Flat.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c
deleted file mode 100644
index 15c4a9d..0000000
--- a/src/sor_confmetric.c
+++ /dev/null
@@ -1,474 +0,0 @@
- /*@@
-   @file      sor_confmetric.c
-   @date      Tue Sep 26 11:29:18 2000
-   @author    Gerd Lanfermann
-   @desc 
-     Provides sor solver engine routines
-   @enddesc 
- @@*/
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-#include "Boundary.h"
-#include "Symmetry.h"
-#include "Ell_DBstructure.h"
-#include "EllBase.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_confmetric_c)
-
-#define SQR(a) ((a)*(a))
-
-int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
-		       int FieldIndex, int MIndex, int NIndex,
-		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
-
- /*@@
-   @routine    sor_confmetric
-   @date       Tue Sep 26 11:28:08 2000
-   @author     Joan Masso, Paul Walker, Gerd Lanfermann
-   @desc 
-   This is a standalone sor solver engine, 
-   called by the wrapper functions
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-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; 
-  CCTK_REAL *gyz =NULL, *gzz =NULL; 
-  CCTK_REAL *psi =NULL;
-  CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
-  CCTK_REAL *var =NULL; 
-
-  /* The inverse metric, allocated here */
-  CCTK_REAL *uxx=NULL, *uyy=NULL, 
-            *uzz=NULL, *uxz=NULL,
-            *uxy=NULL, *uyz=NULL;
-
-  /* shortcuts for metric, psi, deltas, etc. */
-  CCTK_REAL pm4, p12, det;
-
-  CCTK_REAL dx,dy,dz;
-  CCTK_REAL dxdx, dydy, dzdz, 
-            dxdy, dxdz, dydz;
-
-  /* Some physical variables */
-  int accel_cheb=0, accel_const=0;
-  int chebit;
-  CCTK_REAL omega, resnorm=0.0, residual=0.0;
-  CCTK_REAL glob_residual=0.0, rjacobian=0.0; 
-  CCTK_REAL finf;
-  CCTK_INT npow;
-  CCTK_REAL tol;
-
-
-  /* Iteration / stepping  variables */
-  int sorit; 
-  int i,is,ie;
-  int j,js,je;
-  int k,ks,ke,kstep;
-  int nxyz;
-  
-  /* 19 point stencil index */
-  int ijk;
-  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
-  int ipjpk, ipjmk, imjpk, imjmk;
-  int ipjkp, ipjkm, imjkp, imjkm;
-  int ijpkp, ijpkm, ijmkp, ijmkm;
- 
-  /* 19 point stencil coefficients */
-  CCTK_REAL ac;
-  CCTK_REAL ae,aw,an,as,at,ab;
-  CCTK_REAL ane, anw, ase, asw, ate, atw, abe, abw;
-  CCTK_REAL atn, ats, abn, absol;
-
-  /* Miscellaneous */
-  int sum_handle=-1;
-  int sw[3], ierr;
-  int Mstorage=0, Nstorage=0;
-  size_t varsize;
-  static int firstcall = 1;
-  CCTK_REAL detrec_pm4, sqrtdet;
-
-
-  /* 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")) 
-  { 
-    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 (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
-      {
-        printf("SOR with unaccelearted relaxation (omega = 1)\n");
-      }
-    }
-    firstcall = 0;
-  }
-   
-
-  /* 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   */
-  timelevel = 0;
-  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);
-
-  timelevel = 0;
-  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)
-  {
-    timelevel = 0;
-    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)  
-  { 
-    timelevel = 0;
-    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
-    Mstorage = 1;
-  }
-  if (NIndex>=0) 
-  {
-    timelevel = 0;
-    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
-    Nstorage = 1;
-  }
-
-  /* Shortcuts */
-  dx   = GH->cctk_delta_space[0];
-  dy   = GH->cctk_delta_space[1];
-  dz   = GH->cctk_delta_space[2];
-  nxyz = GH->cctk_lsh[0]*GH->cctk_lsh[1]*GH->cctk_lsh[2];
-
-
-  /* Allocate the inverse metric */
-  varsize = (size_t)CCTK_VarTypeSize(CCTK_VarTypeI(FieldIndex))*nxyz;
-  if (!(uxx = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-  if (!(uyy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-  if (!(uzz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-  if (!(uxy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-  if (!(uxz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-  if (!(uyz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
-
-  /* PreCalc the differential coeff. */
-  dxdx = 1.0/(2.0*dx*dx);
-  dydy = 1.0/(2.0*dy*dy);
-  dzdz = 1.0/(2.0*dz*dz);
-  dxdy = 1.0/(4.0*dx*dy);
-  dxdz = 1.0/(4.0*dx*dz);
-  dydz = 1.0/(4.0*dy*dz);
-
-  /* Calculate the inverse metric */
-
-  for (ijk=0;ijk<nxyz;ijk++) 
-  {
-
-    /* determinant */
-    det = -(SQR(gxz[ijk])*gyy[ijk]) + 
-      2*gxy[ijk]*gxz[ijk]*gyz[ijk] - 
-      gxx[ijk]*SQR(gyz[ijk])  -
-      SQR(gxy[ijk])*gzz[ijk] + 
-      gxx[ijk]*gyy[ijk]*gzz[ijk];
-    
-    if (conformal) 
-    {
-      pm4 = 1.0/pow(psi[ijk],4.0);
-      p12 = pow(psi[ijk],12.0);
-    } 
-    else 
-    {
-      pm4 = 1.0;
-      p12 = 1.0;
-    }
-
-
-    /* try to avoid constly division */
-    detrec_pm4 = 1.0/det*pm4;
-    sqrtdet    = sqrt(det);
-
-    uxx[ijk]=(-SQR(gyz[ijk])     + gyy[ijk]*gzz[ijk])*detrec_pm4;
-    uxy[ijk]=( gxz[ijk]*gyz[ijk] - gxy[ijk]*gzz[ijk])*detrec_pm4;
-    uxz[ijk]=(-gxz[ijk]*gyy[ijk] + gxy[ijk]*gyz[ijk])*detrec_pm4;
-    uyy[ijk]=(-SQR(gxz[ijk])     + gxx[ijk]*gzz[ijk])*detrec_pm4;
-    uyz[ijk]=( gxy[ijk]*gxz[ijk] - gxx[ijk]*gyz[ijk])*detrec_pm4;
-    uzz[ijk]=(-SQR(gxy[ijk])     + gxx[ijk]*gyy[ijk])*detrec_pm4;
- 	
-    det    = det*p12;
-    sqrtdet= sqrt(det);
-   
-    /* Rescaling for the uppermetric solver */
-    if (Mstorage) Mlin[ijk] = Mlin[ijk]*sqrt(det);
-    if (Nstorage) Nlin[ijk] = Nlin[ijk]*sqrt(det);
-    
-    uxx[ijk]=uxx[ijk]*dxdx*sqrtdet;
-    uyy[ijk]=uyy[ijk]*dydy*sqrtdet;
-    uzz[ijk]=uzz[ijk]*dzdz*sqrtdet;
-    uxy[ijk]=uxy[ijk]*dxdy*sqrtdet;
-    uxz[ijk]=uxz[ijk]*dxdz*sqrtdet;
-    uyz[ijk]=uyz[ijk]*dydz*sqrtdet;
-
-  }	
-
-  /*$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++)   
-      {
-	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]);
-      }
-    }
-  }$*/
-	
-
-  /* iteration boundaries (ks set inside loop)*/
-  is = 1;
-  js = 1;
-  ie = GH->cctk_lsh[0]-1;
-  je = GH->cctk_lsh[1]-1;
-  ke = GH->cctk_lsh[2]-1;
-  kstep = 2;
-
-  /* start at 1 for historic (Fortran) reasons */
-  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 = 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];
-	  }
-    
-	  resnorm  = resnorm + fabs(residual);
-	  var[ijk] = var[ijk] - omega*residual/ac; 
-	}
-      }
-    }
-
-    /* reduction operation on processor-local residual values */
-    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); 
-  if (uzz) free(uzz);
-  if (uxy) free(uxy); 
-  if (uxz) free(uxz); 
-  if (uyz) free(uyz);
-
-  return retval;
-}
-
-  
diff --git a/src/sor_flat.c b/src/sor_flat.c
deleted file mode 100644
index 1b30f1d..0000000
--- a/src/sor_flat.c
+++ /dev/null
@@ -1,321 +0,0 @@
- /*@@
-   @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>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-#include "Boundary.h"
-#include "Symmetry.h"
-#include "Ell_DBstructure.h"
-#include "EllBase.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_flat_c)
- 
-int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
-		 CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
-
-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;
-  CCTK_REAL *Nlin=NULL;   
-  CCTK_REAL *var =NULL; 
-
-  /* shortcuts for deltas,etc. */
-  CCTK_REAL dx,dy,dz; 
-
-  /* Some physical variables */
-  int accel_cheb=0, accel_const=0;
-  int chebit;
-  CCTK_REAL omega, resnorm, residual, glob_residual, rjacobian; 
-  CCTK_REAL finf;
-  CCTK_INT npow;
-  CCTK_REAL tol;
-  CCTK_REAL zero=0.0;
-
-  /* Iteration / stepping  variables */
-  int sorit; 
-  int i,is,ie;
-  int j,js,je;
-  int k,ks,ke,step;
-  int timelevel;
-
-  /* stencil index */
-  int ijk;
-  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
-
-  /* Coeeficients for the stencil...  */
-  CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab;
- 
-  /* Miscellaneous */
-  int sum_handle=-1;
-  int sw[3], ierr;
-  int Mstorage=0, Nstorage=0;
-  static int firstcall = 1;
-  CCTK_REAL  dx2rec, dy2rec, dz2rec; 
-  char *ans;
-
-  /* Get the reduction handle */
-  sum_handle = CCTK_ReductionArrayHandle("sum");
-  if (sum_handle<0) 
-  {
-    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) */
-
-  sw[0]=1; 
-  sw[1]=1; 
-  sw[2]=1;
-
-  Ell_GetStrKey(&ans,"EllLinFlat::Bnd::Robin");
-
-  if (CCTK_Equals(ans,"yes"))
-  {
-    CCTK_WARN(1,"Robin boundary conditions are not working properly yet");
-    ierr = Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf");
-    ierr = Ell_GetIntKey (&npow, "EllLinFlat::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 (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
-      {
-	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  */
-  /* get the current time level */
-
-  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex);  
-
-  if (MIndex>=0)  
-  { 
-    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
-    Mstorage = 1;
-  }
-  else
-  {
-    CCTK_WARN(0,"Stop");
-  }
-  if (NIndex>=0) 
-  {
-    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex);
-    Nstorage = 1;
-  }
-  else
-  {
-    CCTK_WARN(0,"Stop");
-  }
-
-  /* Shortcuts */
-  dx   = GH->cctk_delta_space[0];
-  dy   = GH->cctk_delta_space[1];
-  dz   = GH->cctk_delta_space[2];
-  
-  dx2rec = 1.0/(dx*dx);
-  dy2rec = 1.0/(dy*dy);
-  dz2rec = 1.0/(dz*dz);
-
-  ae = dx2rec;
-  aw = dx2rec;
-  an = dy2rec;
-  as = dy2rec;
-  at = dz2rec;
-  ab = dz2rec;
-
-  ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec;
-
-
-  ie = GH->cctk_lsh[0]-1;
-  je = GH->cctk_lsh[1]-1;
-  ke = GH->cctk_lsh[2]-1;
-
-  /* start at 1 for historic (Fortran) reasons */
-  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;    
-    
-    is = 1;
-    js = 1;
-    ks = 1;
-
-    /*is = ((GH->cctk_lbnd[0])%2)+(sorit%2)+1;
-      js = ((GH->cctk_lbnd[1])%2)+(sorit%2)+1;
-      ks = ((GH->cctk_lbnd[2])%2)+(sorit%2)+1;*/
-
-    /*step = 2;*/
-    step = 1;
-
-    for (k=ks;k<ke;k=k+step) 
-    {
-      for (j=js;j<je;j=j+step)     
-      {
-	for (i=is;i<ie;i=i+step)    
-        {
-
-	  ac = ac_orig;
-
-	  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);
-	  
-	  if (Mstorage)
-	  {
-	    ac += Mlin[ijk];
-	  }
-
-	  residual = ac * var[ijk]
-	      + ae *var[ipjk]  +  aw*var[imjk]
-	      + an *var[ijpk]  +  as*var[ijmk]
-	      + at *var[ijkp]  +  ab*var[ijkm];
-
-	  if (Nstorage)
-	  {
-	    residual +=Nlin[ijk];
-	  }
-
-	  resnorm  = resnorm + fabs(residual);
-
-	  var[ijk] = var[ijk] - omega*residual/ac; 
-
-	}
-      }
-    }
-
-    ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
-				&resnorm, &glob_residual, CCTK_VARIABLE_REAL);
-    if (ierr<0) 
-    {
-      CCTK_WARN(1,"Reduction of Norm failed");
-    }
-
-#ifdef DEBUG_ELLIPTIC
-    printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
-#endif
-
-    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
-    /*    CCTK_WARN(0,"stop");*/
-
-    /* apply symmetry boundary conditions within loop */    
-    if (CartSymVI(GH,FieldIndex)<0)
-    { 
-      CCTK_WARN(1,"sor_flat: CartSymVI failed in EllSOR loop");
-      break;
-    }
-
-    BndScalarVI(GH,sw,zero,FieldIndex);
-    /* apply boundary conditions within loop */
-    /* FIXME */
-    /*
-    if (BndRobinVI(GH, sw, finf, npow,  FieldIndex)<0)
-    { 
-      CCTK_WARN(1,"sor_flat: BndRobinVI failed in EllSOR loop");
-      break;
-    }
-    */
-
-    /* synchronization of grid variable */
-    CCTK_SyncGroupWithVarI(GH, FieldIndex);
-
-    /* Leave iteration loop if tolerance criterium is met */
-    if (glob_residual<tol)
-    {
-      break;
-    }
-    
-  }
-
-  if (elliptic_verbose)
-  {
-    CCTK_VInfo("EllSOR","Required SOR tolerence was set at %f",tol);
-    CCTK_VInfo("EllSOR","Achieved SOR residual was %f",glob_residual);
-    CCTK_VInfo("EllSOR","Number of iterations was %d (max: %d)",sorit,maxit);
-  }
-  
-  if (glob_residual>tol) 
-  {
-    CCTK_WARN(1,"sor_flat: SOR Solver did not converge");
-    retval = ELL_NOCONVERGENCE;
-  }
-
-  if (ans) free(ans);
-  return retval;
-}
- 
-