Working on SOR. Seems to work, but get different answers with different

numbers of processors, so can't commit a testsuite.


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

4 files changed, 89 insertions, 76 deletions

diff --git a/src/Startup.c b/src/Startup.c
index 80449f3..dfe7d5f 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -10,6 +10,10 @@
 #include "CactusElliptic/EllBase/src/EllBase.h"
 #include "CactusElliptic/EllBase/src/Ell_DBstructure.h"
 
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusElliptic_EllSOR_Startup_c)
+
 void EllSOR_Register(CCTK_ARGUMENTS);
 void sor_confmetric(cGH *GH, 
 	            int *MetricPsiI, 
@@ -55,8 +59,8 @@ void EllSOR_Register(CCTK_ARGUMENTS)
      the routines that sets them ! */
 
   /* Register boundary SOR can handle */
-  err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllSOR::Bnd::Robin");
-  err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllSOR::Bnd::Const");
+  err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Robin");
+  err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Const");
 
   
 
diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c
index 6fb677b..15c4a9d 100644
--- a/src/sor_confmetric.c
+++ b/src/sor_confmetric.c
@@ -21,6 +21,10 @@
 #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,
@@ -167,18 +171,10 @@ int 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   */
-  timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
+  timelevel = 0;
   var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);
 
-  timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[0]) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
+  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]);
@@ -188,11 +184,7 @@ int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
 
   if (conformal)
   {
-    timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[6]) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
+    timelevel = 0;
     psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[6]);
   }
 
@@ -201,21 +193,13 @@ int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
      set flag Mstorage=1, dito for N */
   if (MIndex>=0)  
   { 
-    timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
+    timelevel = 0;
     Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
     Mstorage = 1;
   }
   if (NIndex>=0) 
   {
-    timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
+    timelevel = 0;
     Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
     Nstorage = 1;
   }
diff --git a/src/sor_flat.c b/src/sor_flat.c
index 752277b..1b30f1d 100644
--- a/src/sor_flat.c
+++ b/src/sor_flat.c
@@ -21,6 +21,10 @@
 #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);
@@ -48,13 +52,13 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   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,kstep;
-  int nxyz;
+  int k,ks,ke,step;
   int timelevel;
 
   /* stencil index */
@@ -70,7 +74,8 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   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) 
@@ -81,17 +86,20 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   /* 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");
-  }
+  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"))
@@ -128,39 +136,32 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      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 */
-  timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
-  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);  
+
+  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex);  
 
   if (MIndex>=0)  
   { 
-    timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
-    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
+    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
     Mstorage = 1;
   }
+  else
+  {
+    CCTK_WARN(0,"Stop");
+  }
   if (NIndex>=0) 
   {
-    timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
-    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
+    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];
-  nxyz = GH->cctk_lsh[0]*GH->cctk_lsh[1]*GH->cctk_lsh[2];
   
   dx2rec = 1.0/(dx*dx);
   dy2rec = 1.0/(dy*dy);
@@ -175,12 +176,10 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 
   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++) 
@@ -200,21 +199,26 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
       omega = 1.8;
     }
     
-    resnorm = 0.0;
+    resnorm = 0.0;    
     
-    ks = (sorit%2)+1;
-    if (GH->cctk_lsh[2]==3)
-    {
-      ks = 2;
-    }
+    is = 1;
+    js = 1;
+    ks = 1;
 
-    for (k=ks;k<ke;k+=kstep) 
+    /*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++)     
+      for (j=js;j<je;j=j+step)     
       {
-	for (i=is;i<ie;i++)    
+	for (i=is;i<ie;i=i+step)    
         {
-	  
+
 	  ac = ac_orig;
 
 	  ijk   = CCTK_GFINDEX3D(GH,i  ,j  ,k  );
@@ -248,7 +252,6 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
       }
     }
 
-    /* reduction operation on processor-local residual values */
     ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
 				&resnorm, &glob_residual, CCTK_VARIABLE_REAL);
     if (ierr<0) 
@@ -256,6 +259,10 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
       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]);
 
@@ -263,19 +270,25 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 #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,"CartSymVI failed in EllSOR loop");
+      CCTK_WARN(1,"sor_flat: CartSymVI failed in EllSOR loop");
       break;
     }
 
+    BndScalarVI(GH,sw,zero,FieldIndex);
     /* apply boundary conditions within loop */
-    if (CCTK_EQUALS(sor_bound,"robin"))
-    {
-      ierr = BndRobinVI(GH, sw, finf, npow,  FieldIndex);
+    /* 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);
@@ -287,6 +300,13 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     }
     
   }
+
+  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) 
   {
@@ -294,6 +314,7 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     retval = ELL_NOCONVERGENCE;
   }
 
+  if (ans) free(ans);
   return retval;
 }
  
diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c
index 062a846..3df20c5 100644
--- a/src/sor_wrapper.c
+++ b/src/sor_wrapper.c
@@ -24,6 +24,10 @@
 
 #include "CactusElliptic/EllBase/src/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,