From 258ae09276f7737e50c5d40449679aa62f4e50aa Mon Sep 17 00:00:00 2001
From: lanfer <lanfer@fa3da13c-9f13-4301-a575-cf5b8c5e1907>
Date: Tue, 26 Sep 2000 10:55:23 +0000
Subject: SOR solvers for flat and conformal metric

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@59 fa3da13c-9f13-4301-a575-cf5b8c5e1907
---
 src/sor_confmetric.c | 413 +++++++++++++++++++++++++++++++++++++++++++++++++++
 src/sor_flat.c       | 219 +++++++++++++++++++++++++++
 2 files changed, 632 insertions(+)
 create mode 100644 src/sor_confmetric.c
 create mode 100644 src/sor_flat.c

diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c
new file mode 100644
index 0000000..b0df762
--- /dev/null
+++ b/src/sor_confmetric.c
@@ -0,0 +1,413 @@
+ /*@@
+   @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 "CactusBase/Boundary/src/Boundary.h"
+#include "CactusBase/CartGrid3D/src/Symmetry.h"
+#include "CactusElliptic/EllBase/src/Ell_DBstructure.h"
+
+
+#define SQR(a) ((a)*(a))
+
+ /*@@
+   @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 
+
+@@*/
+
+
+
+void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
+		       int FieldIndex, int MIndex, int NIndex,
+		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
+{
+  DECLARE_CCTK_PARAMETERS  
+
+  /* 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;
+  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, abs;
+
+  /* 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   */
+  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]);
+  
+
+  if (conformal)
+    psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, 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);
+    Mstorage = 1;
+  }
+  if (NIndex>=0) {
+    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,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];
+	    abs = 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] + abs*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) 
+    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;
+}
+
+  
diff --git a/src/sor_flat.c b/src/sor_flat.c
new file mode 100644
index 0000000..0117894
--- /dev/null
+++ b/src/sor_flat.c
@@ -0,0 +1,219 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "cctk.h"
+#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"
+ 
+
+void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
+		 CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
+{
+
+  DECLARE_CCTK_PARAMETERS
+
+  /* The pointer to the data fields */
+  CCTK_REAL *Mlin=NULL, *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;
+  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;
+
+  /* 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; 
+  
+
+  /* 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; 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);
+    Mstorage = 1;
+  }
+  if (NIndex>=0) {
+    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,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];
+  
+  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;
+
+  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 = 2;
+    
+    for (k=ks;k<ke;k+=kstep) {
+      for (j=js;j<je;j++)     {
+	for (i=is;i<ie;i++)    {
+	  
+	  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; 
+
+	  printf(" %d %d %d  %f \n",i,j,k,var[ijk]);
+
+	}
+      }
+    }
+
+    /* 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 Norm 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;
+    
+  }
+  
+  if (glob_residual>tol) 
+    CCTK_WARN(2,"SOR SOLVER DID NOT CONVERGE");
+
+  return;
+}
+ 
+  
-- 
cgit v1.2.3