less beta

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@10 1d96b42b-98df-4a6a-9d84-1b24288d4588

1 files changed, 84 insertions, 171 deletions

diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c
index 890289a..1bbefc0 100644
--- a/src/petsc_flat_solver.c
+++ b/src/petsc_flat_solver.c
@@ -1,15 +1,12 @@
 
 #include "cctk.h"
 #include "cctk_parameters.h"
-#include "cctk_Flesh.h"
 
-#include "assert.h"
 #include "mpi.h"
 #include "sles.h"
 
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
-
 #define DEBUG
 
 /*Don't know what these actually mean! FIXME */
@@ -17,47 +14,36 @@
 #define ELLCONV_RELATIVE 1
 #define ELLCONV_EITHER   2
 
-/* A few convenient macros */
+/* A convenient macros for the stencil molecule*/
 #define a(i,j,k) a[i+1 + 3*(j+1 + 3*(k+1))]
-#define SQR(a) ((a)*(a))
-
   
 
-/* Place the matrix in global space so we can keep it around
-   with the memory stripped out. Thanks, Barry!  */
-
 static int trips=0;
 static Mat     *A;           /* linear system matrix */
 static Vec     soln, b;      /* approx solution, RHS */
 static SLES    sles;         /* linear solver context */
 
 
-
 void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, 
-		       int *AbsTol, int *RelTol) {
-
+		int *AbsTol, int *RelTol) {
 
   DECLARE_CCTK_PARAMETERS    /* CCTK passed parameters */
 
-
-  void *GetDataPtr_NextTL(cGH *GH, const char *field);
-  int  CCTK_GetArgc(void);        /* PETSc wnats to know the prog. name */
-  void CCTK_GetArgv(char **argv); /* and parameters, don't think they can be passed through */
-
-
   PC     pc;            /* preconditioner context */
   KSP    ksp;           /* Krylov subspace method context */
   int    num_A;         /* Number of A-arrays as needed ny Dan's MG */
   
 
-  double norm;          /* norm of solution error */
   int    ierr;          /* Check the return status of petsc */
   int    retcode;       /* Check the return status of CCTK */
 
-  double  a[27];        /* The stencil array */
-  int     idx[27];      /* The column of the stencil array */
-  int     rank;         /* Rank of the matrix/vector */
-  int     its;          /* Number of iterations */
+  int    idx[27];      /* The column of the stencil array */
+  int    rank;         /* Rank of the matrix/vector */
+  int    its;          /* Number of iterations */
+  
+  CCTK_REAL norm;       /* norm of solution error */
+  CCTK_REAL a[27];      /* The stencil array */
+  CCTK_REAL ac;         /* Storage for a(0,0,0) for renorm */
   
   /* Loopers */
   int i,j,k,l,m,n,ll;
@@ -65,25 +51,18 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   /* loop limits put in for stencil_w !=1     Ed Evans 1/19/98 */
   int imin,imax,jmin,jmax,kmin,kmax,interior;  
    
-  pGH *pughGH;                 /* The pugh Extension handle */
-  int myproc;                  /* out processor */
-
+  pGH *pughGH;          /* The pugh Extension handle */
+ 
   /* Tolerances */
   CCTK_REAL rtol, atol, tolerance;
-    
+
   /* Values to assemble the matrix */
-  CCTK_REAL two=2.0, four=4.0, tmp;
-  CCTK_REAL pm4, psixp,psiyp,psizp,det;
+  CCTK_REAL tmp, *values;
+
+  /* Misc */
   CCTK_REAL dx,dy,dz;
-  CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz;
-  CCTK_REAL Gxxx,Gxxy,Gxxz,Gxyy,Gxyz,Gxzz; /* Christoffels */
-  CCTK_REAL Gyxx,Gyxy,Gyxz,Gyyy,Gyyz,Gyzz;
-  CCTK_REAL Gzxx,Gzxy,Gzxz,Gzyy,Gzyz,Gzzz;
-  CCTK_REAL dxxx,dxxy,dxxz,dxyy,dxyz,dxzz;
-  CCTK_REAL dyxx,dyxy,dyxz,dyyy,dyyz,dyzz;
-  CCTK_REAL dzxx,dzxy,dzxz,dzyy,dzyz,dzzz;
-  CCTK_REAL *values;
-  
+
+  int myproc; 
   int nxs,nys,nzs;      /* Size of the grid with stencils off... */
   int startpoint=0;     /* My starting index (per proc) */
   int endpoint;         /* One more than my end */
@@ -96,24 +75,25 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   int nabla_form;       /* Which form of the nable */
   int reuse;            /* Do PETSc reuse tricks? */
   int matnormalize;     /* Normalize the central mat value to one? */
-  CCTK_REAL ac;         /* Storage for a(0,0,0) for renorm */
   int pctype, ksptype;  /* Which PC and KSP to use. See below for
                            monster nested if statement. */
+  int I,J;              /* The PETSc col and row in the matrix */
+	    
   int PetscTolStyle;
 
-  /* For the upper metric form of nabla */
-  CCTK_REAL *uxx3, *uxy3, *uxz3, *uyy3, *uyz3, *uzz3;
-
   
   /* Pointers to the data of : petsc workspace/ the GF to solve /
-     metric / psi / derivs(psi) / Mlinear / Nlinear(source) */
+     Mlinear / Nlinear(source) */
   CCTK_REAL *wsp =NULL, *ell_field=NULL; 
- 
   CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
 
+  /* flags to signal if storage is ON or OFF (used/not-used) for M/N*/
   int Mstorage=0, Nstorage=0;
 
 
+
+
+
   octant       = CCTK_Equals(domain,"octant");
   verbose      = CCTK_Equals(petsc_verbose,"yes")||
                  CCTK_Equals(petsc_verbose,"debug");
@@ -121,6 +101,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   reuse        = 0; 
   matnormalize = 0; 
   
+
   /* FIXME, the TolAbs/TolRel will be evaluated here */
   PetscTolStyle=0;
   tolerance    =AbsTol[0];
@@ -149,15 +130,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
 
   trips++;
 
+  /* Allocate storage for system matrix, will have more of them later */
   num_A = 1;
   A=(Mat*) malloc (sizeof(Mat)*num_A);
-  
-  /* 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 */
-  ell_field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex);
+
+  /* Get the data pointers */
   
   /* if we have a negative index, this GF is not needed, 
      there for don't even look for it. when index positive,
@@ -170,10 +148,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
     Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex);
     Nstorage = 1;
   }
-  wsp   = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,CCTK_VarIndex("ellpetsc::wsp"));
+  ell_field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex);
+  wsp       = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,CCTK_VarIndex("ellpetsc::wsp"));
   
-  /* initialize the linear index lookup table, after it's 
-     filled up (below), the -1 indicates a boundary */
+
+  /* initialize the linear index lookup table with -1; 
+     this will indicate a boudary later */
   for (i=0;i<pughGH->npoints;i++) wsp[i] = -1.0;
 
   /* Get myproc from the CCTK, get the gridspacing */
@@ -238,14 +218,14 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
     for (j=jmin;j<jmax;j++)
       for (i=imin;i<imax;i++)
         wsp[DATINDEX(pughGH,i,j,k)] = endpoint++;
+
 #ifdef DEBUG  
   printf("ENDPOINT:  %d \n",endpoint);
 #endif DEBUG
 
   /* So now each point has a unique index of its own. If we
      do a sync that means each processor knows the indiced
-     in the ghost zones (eg, on its neighbors) in a FD stencil
-     of 1 */
+     in the ghost zones (eg, on its neighbors) in a FD stencil of 1 */
   retcode = pugh_SyncGroup(GH,"ellpetsc::petscworkspace");
   
   
@@ -289,6 +269,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
      "official", yet */
 
   /*$if (trips == 0 || reuse == 0) {$*/
+  
   if (verbose) CCTK_INFO("Creating Matrices and Vectors ....");
   ierr = MatCreateMPIAIJ(pughGH->PUGH_COMM_WORLD,     
 			 (endpoint-startpoint),   /* # of rows */
@@ -302,6 +283,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   CHKERRA(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
   CHKERRA(ierr);
+  
   /*$}$*/ 
   /*$else {$*/
   /*  ierr = MatRestoreValuesMemory(A);  */
@@ -310,16 +292,27 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   }$*/
 
 
-
-
   /* Compare the PETSc layout to Cactus, this better be a match */
   ierr = VecGetOwnershipRange(soln,&pvstart,&pvend);
   ierr = MatGetOwnershipRange(A[0],&pstart,&pend);
-  printf("CAC M-Layout: %d %d \n",startpoint, endpoint);
-  printf("PET M-Layout: %d %d \n",pstart, pend);
-  printf("PET V-Layout: %d %d \n",pvstart,pvend);
+  if (debug) { 
+    printf("CAC M-Layout: %d %d \n",startpoint, endpoint);
+    printf("PET M-Layout: %d %d \n",pstart, pend);
+    printf("PET V-Layout: %d %d \n",pvstart,pvend);
+  }
   if (pstart != startpoint && pend != endpoint) 
     CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why??)\n");
+         
+
+  /* Initialize the stencil array */
+  for (I=0;I<27;I++) a[I] = 0.0;
+  a( 0, 0, 0) = -10.38;
+  a( 1, 0, 0) = 1.56;  /* ae */
+  a(-1, 0, 0) = 1.56;  /* aw */
+  a( 0, 1, 0) = 1.56;  /* an */       
+  a( 0,-1, 0) = 1.56;  /* as */      
+  a( 0, 0, 1) = 1.56;  /* at */              
+  a( 0, 0,-1) = 1.56;  /* ab */
 
 
   for (k=kmin;k<kmax;k++) {
@@ -329,22 +322,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
         if (wsp[DATINDEX(pughGH,i,j,k)] >= 0) {
 
           /* Set up indices */
-          int ijk;              /* The data point for the array */
+          int ijk;              /* The lineax index for the array */
           int ig, jg, kg;       /* The position in global space */
-          int I,J;              /* The col and row in the matrix */
           CCTK_REAL rhsval = 0;
-          
-          for (I=0;I<27;I++) a[I] = 0.0;
-	  a( 0, 0, 0) = -10.38;
-	  a( 1, 0, 0) = 1.56;  /* ae */
-	  a(-1, 0, 0) = 1.56;  /* aw */
-	  a( 0, 1, 0) = 1.56;  /* an */       
-	  a( 0,-1, 0) = 1.56;  /* as */      
-	  a( 0, 0, 1) = -1.56;  /* at */              
-	  a( 0, 0,-1) = -1.56;  /* ab */
-  
-          /* these guys are easy */
-          ijk = DATINDEX(pughGH,i,j,k);               /* get linear index */
+
+          /* local linear index / global 3d index */
+          ijk = DATINDEX(pughGH,i,j,k);    
           ig  = i + GH->cctk_lbnd[0]; 
           jg  = j + GH->cctk_lbnd[1];
           kg  = k + GH->cctk_lbnd[2];
@@ -353,45 +336,19 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
           I = wsp[ijk];
        
           /* M phi */
-	  if (Mstorage)
-	    a(0,0,0) = a(0,0,0) - Mlin[ijk];
+	  if (Mstorage) a(0,0,0) = a(0,0,0) - Mlin[ijk];
 
-          /* Great now set the values of the matrix. This is
-             really painful due to the boundaries (here we force
-             dirichlet).  */
+          /* Great now set the values of the matrix. Here we force dirichlet. */
 	  VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
           
-	  /* Put in the octant boundary conditions.
-           
-             We do these by reflecting the -1 stencil if we are at a
-             boundary at -1. That is, imagine the 1D laplace equation
-             with the boundary condition a(0) = a(1). This means
-             the point 1 stencil (which is our first stencil in
-             the matrix) becomes
-             
-             a(0) + a(2) - 2 a(1) = rho(i) * deltax  ->
-             a(2) - 2 a(1) + a(1) = rho(i) * deltax
-             
-             OK so think about this with a general stenci
-             
-             S_0 a_0 + S_1 a_1 + S_2 a_2 = rho_1 ->
-             S_0 a_1 + S_1 a_1 + S_2 a_2 = rho_1 ->
-             (S_0 + S_1) a_1 + S_2 a_2 = rho_1
-             
-             eg, S_0 -> 0 and S_1 -> S_1 + S_0
-             
-             Great, so implement this in 3D. It is a bit trickier,
-             since we don't want to zero neighbors in the wrong
-             direction, but not impossible.
-             
-            */
-
+#ifdef DEBUG
 	  if (ijk==1918) {
 	    for (l=-1;l<=1;l++)
 	      for (m=-1;m<=1;m++)
 		for (n=-1;n<=1;n++) printf(" (%d %d %d): %f \n",l,m,n,a(l,m,n));
 	    printf("\n");
 	  }
+#endif
 
           if (octant)
             for (l=-1;l<=0;l++)
@@ -417,7 +374,8 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
             ac = a(0,0,0);
             for (J=0;J<27;J++) a[J] = a[J] / ac;
           } 
-          else ac = 1;
+          else 
+	    ac = 1;
           
 
           /* This is the new way-clever look. Note it relies
@@ -442,17 +400,19 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
                   }
                 }
               }
-
-	  if (Nstorage)
-	    rhsval = -rhsval + Nlin[ijk] / ac;
+	  
+	  if (Nstorage) rhsval = -rhsval + Nlin[ijk] / ac;
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
 	  CHKERRA(ierr);
           
-        }
-      }
-    }
-  }
+	  /* reset the central value*/
+	  a( 0, 0, 0) = -10.38;
+  
+        }   /* if (wsp>=0) */
+      }   /*for i */
+    }   /* for j */
+  }   /* for k */
   
   if (verbose) CCTK_INFO ("Assembling the vectors");
   ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
@@ -462,21 +422,11 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   ierr = VecAssemblyEnd(soln);                      CHKERRA(ierr);
   ierr = VecAssemblyEnd(b);                         CHKERRA(ierr);
   ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(ierr);
-
-  if (nabla_form == 3) {
-    if (verbose)
-	printf ("Freeing upper metric storage\n");
-    free(uxx3);
-    free(uxy3);
-    free(uxz3);
-    free(uyy3);
-    free(uyz3);
-    free(uzz3);
-  }
+  
   if (trips==0) 
     OptionsSetValue("-ksp_monitor","");
   
-  CCTK_INFO("CREATING SLES");
+  if (verbose) CCTK_INFO("CREATING SLES");
   ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRA(ierr);
   
   
@@ -489,14 +439,14 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
      using A to precondition itself. Since I don't know what this
      means, we'll leave it for now.  */
 
-  CCTK_INFO("CREATING SLES OPERATOR");
+  if (verbose) CCTK_INFO("CREATING SLES OPERATOR");
   ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr);
 
   /* Set linear solver defaults for this problem. Later this
      should be replaced/modified with appropriate parsing from the 
      parameter parser. For now it is not. These defaults are
      reasonable, I hope.  */
-  CCTK_INFO("SLESGet KSP/PC");
+  if (verbose) CCTK_INFO("SLESGet KSP/PC");
   ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr);
   ierr = SLESGetPC(sles,&pc);   CHKERRA(ierr);
   
@@ -567,48 +517,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
             PetscTolStyle);
   }
   
-  ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT,
-                          PETSC_DEFAULT);
+  ierr = KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, PETSC_DEFAULT);
   CHKERRA(ierr);
 
-
-    /* We are warned in the manual that 
-
-     The default technique for orthogonalization of the Hessenberg matrix
-     in GMRES is the modified Gram-Schmidt method, which employs many
-     VecDot() operations and can thus be slow in parallel. A fast approach
-     is to use the unmodified Gram-Schmidt method, which can be set with
-
-     ierr = KSPGMRESSetOrthogonalization(KSP ksp, 
-        KSPGMRESUnmodifiedGramSchmidtOrthogonalization); 
-
-     or the options database command
-     -ksp_gmres_unmodifiedgramschmidt. Note that this algorithm is
-     numerically unstable, but may deliver much better speed
-     performance. One can also use unmodifed Gram-Schmidt with
-     iterative refinement, by setting the orthogonalization routine,
-     KSPGMRESIROrthog(), by using the command line option
-     -ksp_gmres_irorthog.
-
-     So this my not be a smart thing to do. But lets put it here 
-     so we can turn it on or off later.
-
-   */
-  /*$ierr = KSPGMRESSetOrthogonalization(ksp, 
-                KSPGMRESUnmodifiedGramSchmidtOrthogonalization); 
-  CHKERRA(ierr);$*/
-
-  /* We also learn that  ...
-
-     By default, KSP assumes an initial guess of zero by zeroing the
-     initial value for the solution vector that is given. To use a
-     nonzero initial guess, the user must call
-
-     ierr = KSPSetInitialGuessNonzero(KSP ksp); 
-
-   */
-  CCTK_INFO("KSPSetInitialGuess\n");
-  ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr);
+  if (verbose) CCTK_INFO("KSPSetInitialGuess\n");
+  ierr = KSPSetInitialGuessNonzero(ksp); 
+  CHKERRA(ierr);
   
   /* 
     Set runtime options. Since we don't use PETSC runtime options
@@ -624,7 +538,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-  CCTK_INFO("SLES solve...");
+  if (verbose) CCTK_INFO("SLES solve...");
   ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); 
   
   /* Here we can form a "res = Ax  - b" and find its norm to get
@@ -659,7 +573,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
           /* Now this one. "Fortran-order" the matrix. But remember
              we have stripped off the ghost zones. Hence ig-1 and nxs...
           */
-          I =wsp[ijk]-startpoint;
+          I = wsp[ijk]-startpoint;
           ell_field[ijk] = values[I];
         }
       }
@@ -668,11 +582,10 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex,
   
   VecRestoreArray(soln,&values);   
 
-  /* Sync var in CCTK speak. */ 
   pugh_SyncGroup(GH,"ellpetsc::petscworkspace");
   
-  /* And finally free up the matrix memory   FIXME 
-  ierr = MatReleaseValuesMemory(A); CHKERRA(ierr); */
+  /* And finally free up the matrix memory */  
+  /*$ierr = MatReleaseValuesMemory(A); CHKERRA(ierr);$*/ 
 
   /* This code is not used anymore */
   if (verbose) CCTK_INFO("Destroying Matrices");