Untabified. Please try to adhere to the Cactus coding guidelines when

modifying or adding files.


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

4 files changed, 200 insertions, 202 deletions

diff --git a/src/Startup.c b/src/Startup.c
index 05212b1..bdadc2c 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -12,13 +12,13 @@ static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusElliptic_EllPETSc_Startup_c)
 
 void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI,
-		      int *NI, int *AbsTol, int *RelTol);
+                      int *NI, int *AbsTol, int *RelTol);
 
 void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI,
-		  int *NI, int *AbsTol, int *RelTol);
+                  int *NI, int *AbsTol, int *RelTol);
 
 void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, 
-		int *AbsTol, int *RelTol);
+                int *AbsTol, int *RelTol);
 
 /* Registration of the petsc solvers with the Elliptic solver registry.
    This routine registers petsc_confmetric under the name "petsc" for 
diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c
index c4bfb5c..7600fd8 100644
--- a/src/petsc_confmetric_solver.c
+++ b/src/petsc_confmetric_solver.c
@@ -80,8 +80,8 @@ static SLES    sles;         /* linear solver context */
 
 
 int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 
-			     int FieldIndex, int MIndex, int NIndex, 
-			     CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+                             int FieldIndex, int MIndex, int NIndex, 
+                             CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
 void *GetDataPtr_NextTL(cGH *GH, const char *field);
 
 
@@ -105,8 +105,8 @@ void *GetDataPtr_NextTL(cGH *GH, const char *field) {
    elliptic registration routine see LinearElliptic.h*/
 
 int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 
-			     int FieldIndex, int MIndex, int NIndex, 
-			     CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
+                             int FieldIndex, int MIndex, int NIndex, 
+                             CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
 
   DECLARE_CCTK_PARAMETERS    /* CCTK passed parameters */
 
@@ -284,19 +284,19 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   /* We fix the lower and upper boundary indices, that actually "active" in 
      the sense that they are no ghostzones: */
   imin=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
-	 GH->cctk_nghostzones[0]);
+         GH->cctk_nghostzones[0]);
   imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : 
-	 pEx->lnsize[0]-GH->cctk_nghostzones[0]);
+         pEx->lnsize[0]-GH->cctk_nghostzones[0]);
 
   jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
-	GH->cctk_nghostzones[1]);
+        GH->cctk_nghostzones[1]);
   jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : 
-	 pEx->lnsize[1]-GH->cctk_nghostzones[1]);
+         pEx->lnsize[1]-GH->cctk_nghostzones[1]);
   
   kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
-	 GH->cctk_nghostzones[2]);
+         GH->cctk_nghostzones[2]);
   kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : 
-	 pEx->lnsize[2]-GH->cctk_nghostzones[2]);
+         pEx->lnsize[2]-GH->cctk_nghostzones[2]);
 
   /* We need to get the global index of gridpoints (gps) owned on my proc. 
      For that we have to know how many gps are there before my processors
@@ -388,12 +388,12 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   if (verbose) CCTK_INFO("Creating Matrices and Vectors ....");
 
   ierr = MatCreateMPIAIJ(pughGH->PUGH_COMM_WORLD,     
-			 (endpoint-startpoint),   /* # of rows */
-			 (endpoint-startpoint),   /* This is confusing */
-			 rank,rank,               /* Matrix size */
-			 19,PETSC_NULL,           /* Diagonal  */
-			 19,PETSC_NULL,           /* Off diagonal */
-			 &A[0]);                  /* The output  */
+                         (endpoint-startpoint),   /* # of rows */
+                         (endpoint-startpoint),   /* This is confusing */
+                         rank,rank,               /* Matrix size */
+                         19,PETSC_NULL,           /* Diagonal  */
+                         19,PETSC_NULL,           /* Off diagonal */
+                         &A[0]);                  /* The output  */
   CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
   CHKERRQ(ierr);
@@ -421,7 +421,7 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   else 
     if  (CCTK_EQUALS(petsc_nablaform,"up")) {
       if (verbose) 
-	CCTK_INFO("Forming nabla with upper g and finite difference of dg \n");
+        CCTK_INFO("Forming nabla with upper g and finite difference of dg \n");
       nabla_form = 3;
       uxx3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
       uxy3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
@@ -430,42 +430,42 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
       uyz3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
       uzz3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
       for (i=0;i<pEx->npoints;i++) {
-	CCTK_REAL det;
-	CCTK_REAL p12;
-	det= -(SQR(gxz[i])*gyy[i]) + 
-	  2*gxy[i]*gxz[i]*gyz[i] - 
-	  gxx[i]*SQR(gyz[i])  -
-	  SQR(gxy[i])*gzz[i] + 
-	  gxx[i]*gyy[i]*gzz[i];
-	
-	if (conformal) {
-	  pm4 = 1./pow(psi[i],4.0);
-	  p12 = pow(psi[i],12.0);
-	} else {
-	  pm4 = 1.0;
-	  p12 = 1.0;
-	}
+        CCTK_REAL det;
+        CCTK_REAL p12;
+        det= -(SQR(gxz[i])*gyy[i]) + 
+          2*gxy[i]*gxz[i]*gyz[i] - 
+          gxx[i]*SQR(gyz[i])  -
+          SQR(gxy[i])*gzz[i] + 
+          gxx[i]*gyy[i]*gzz[i];
+        
+        if (conformal) {
+          pm4 = 1./pow(psi[i],4.0);
+          p12 = pow(psi[i],12.0);
+        } else {
+          pm4 = 1.0;
+          p12 = 1.0;
+        }
       
-	/*invert metric. This is the conformal upper metric. */
-	uxx3[i]=(-SQR(gyz[i]) + gyy[i]*gzz[i])/det*pm4;
-	uxy3[i]=(gxz[i]*gyz[i] - gxy[i]*gzz[i])/det*pm4;
-	uyy3[i]=(-SQR(gxz[i]) + gxx[i]*gzz[i])/det*pm4;
-	uxz3[i]=(-gxz[i]*gyy[i] + gxy[i]*gyz[i])/det*pm4;
-	uyz3[i]=(gxy[i]*gxz[i] - gxx[i]*gyz[i])/det*pm4;
-	uzz3[i]=(-SQR(gxy[i]) + gxx[i]*gyy[i])/det*pm4;
-	
-	det = det*p12;
+        /*invert metric. This is the conformal upper metric. */
+        uxx3[i]=(-SQR(gyz[i]) + gyy[i]*gzz[i])/det*pm4;
+        uxy3[i]=(gxz[i]*gyz[i] - gxy[i]*gzz[i])/det*pm4;
+        uyy3[i]=(-SQR(gxz[i]) + gxx[i]*gzz[i])/det*pm4;
+        uxz3[i]=(-gxz[i]*gyy[i] + gxy[i]*gyz[i])/det*pm4;
+        uyz3[i]=(gxy[i]*gxz[i] - gxx[i]*gyz[i])/det*pm4;
+        uzz3[i]=(-SQR(gxy[i]) + gxx[i]*gyy[i])/det*pm4;
+        
+        det = det*p12;
       
-	/* Rescaling for the uppermetric solver */
-	if (Mstorage) Mlin[i] = Mlin[i]*sqrt(det);
-	if (Nstorage) Nlin[i] = Nlin[i]*sqrt(det);
-
-	uxx3[i]=uxx3[i]/(2.*dx*dx)*sqrt(det);
-	uyy3[i]=uyy3[i]/(2.*dy*dy)*sqrt(det);
-	uzz3[i]=uzz3[i]/(2.*dz*dz)*sqrt(det);
-	uxy3[i]=uxy3[i]/(4.*dx*dy)*sqrt(det);
-	uxz3[i]=uxz3[i]/(4.*dx*dz)*sqrt(det);
-	uyz3[i]=uyz3[i]/(4.*dy*dz)*sqrt(det);
+        /* Rescaling for the uppermetric solver */
+        if (Mstorage) Mlin[i] = Mlin[i]*sqrt(det);
+        if (Nstorage) Nlin[i] = Nlin[i]*sqrt(det);
+
+        uxx3[i]=uxx3[i]/(2.*dx*dx)*sqrt(det);
+        uyy3[i]=uyy3[i]/(2.*dy*dy)*sqrt(det);
+        uzz3[i]=uzz3[i]/(2.*dz*dz)*sqrt(det);
+        uxy3[i]=uxy3[i]/(4.*dx*dy)*sqrt(det);
+        uxz3[i]=uxz3[i]/(4.*dx*dz)*sqrt(det);
+        uyz3[i]=uyz3[i]/(4.*dy*dz)*sqrt(det);
       }
     } 
     else  CCTK_WARN(0,"Don't know how to form nabla!\n");
@@ -498,22 +498,22 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
           ig  = i + GH->cctk_lbnd[0]; 
           jg  = j + GH->cctk_lbnd[1];
           kg  = k + GH->cctk_lbnd[2];
-	
+        
           /* Get the row we are working on */
           I = wsp[ijk];
 
           /* Setup Temporaries / Psi derivatives on psi */
-	  if (conformal) {
-	    pm4   = 1./pow(psi[ijk],4.0);
-	    psixp = psix[ijk];
-	    psiyp = psiy[ijk];
-	    psizp = psiz[ijk];
-	  } else {
-	    pm4   = 1.0;
-	    psixp = 0.0;
-	    psiyp = 0.0;
-	    psizp = 0.0;
-	  }
+          if (conformal) {
+            pm4   = 1./pow(psi[ijk],4.0);
+            psixp = psix[ijk];
+            psiyp = psiy[ijk];
+            psizp = psiz[ijk];
+          } else {
+            pm4   = 1.0;
+            psixp = 0.0;
+            psiyp = 0.0;
+            psizp = 0.0;
+          }
           if (nabla_form == 2) {
             /* Use finite differences of g for the d's */
             int ijkp, ijkm;
@@ -551,9 +551,9 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
             tdzgyy = (gyy[ijkp] - gyy[ijkm])/(2.0*dz);
             tdzgyz = (gyz[ijkp] - gyz[ijkm])/(2.0*dz);
             tdzgzz = (gzz[ijkp] - gzz[ijkm])/(2.0*dz);
-	    
-	    /* great ... so start hacking away at the coefficients.
-	       Form upper metric - compute determinant */
+            
+            /* great ... so start hacking away at the coefficients.
+               Form upper metric - compute determinant */
             det= -(SQR(gxz[ijk])*gyy[ijk]) + 
               2*gxy[ijk]*gxz[ijk]*gyz[ijk] - 
               gxx[ijk]*SQR(gyz[ijk])  -
@@ -707,64 +707,64 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
               two * uyz * Gzyz;
             a(0,0,1)  = a(0,0,1)  - pm4 * tmp / (two*dz);
             a(0,0,-1) = a(0,0,-1) + pm4 * tmp / (two*dz);
-	      
+              
           } /* end if nable_form==2 */ 
 
-	  else { 
-	    /* nabla_form==3: Upper Metric Nabla Form */
-	    a(0,0,0) = -Uxx(i+1,j,k) -2.*Uxx(i,j,k) -Uxx(i-1,j,k) 
-	      -Uyy(i,j+1,k) -2.*Uyy(i,j,k) -Uyy(i,j-1,k)
-	      -Uzz(i,j,k+1) -2.*Uzz(i,j,k) -Uzz(i,j,k-1);
-	    
-	    /*$ae = uxx(i+1,j,k)+uxx(i,j,k)$*/
-	    a(1,0,0) = Uxx(i+1,j,k) + Uxx(i,j,k);
-	    /*$aw = uxx(i-1,j,k)+uxx(i,j,k)$*/
-	    a(-1,0,0) = Uxx(i-1,j,k)+Uxx(i,j,k);
-		/*$an = uyy(i,j+1,k)+uyy(i,j,k)$*/              
-	    a(0,1,0) = Uyy(i,j+1,k)+Uyy(i,j,k);
-	    /*$as = uyy(i,j-1,k)+uyy(i,j,k)$*/              
-	    a(0,-1,0) = Uyy(i,j-1,k)+Uyy(i,j,k);
-	    /*$at = uzz(i,j,k+1)+uzz(i,j,k)$*/               
-	    a(0,0,1) = Uzz(i,j,k+1)+Uzz(i,j,k);               
-	    /*$ab = uzz(i,j,k-1)+uzz(i,j,k)$*/               
-	    a(0,0,-1) = Uzz(i,j,k-1)+Uzz(i,j,k);
-	    
-	    /*$ane = uxy(i,j+1,k)+uxy(i+1,j,k)$*/              
-	    a(1,1,0) = Uxy(i,j+1,k)+Uxy(i+1,j,k);
-	    /*$anw = - uxy(i-1,j,k)-uxy(i,j+1,k)$*/               
-	    a(-1,1,0) = - Uxy(i-1,j,k)-Uxy(i,j+1,k);
-	    /*$ase = - uxy(i+1,j,k)-uxy(i,j-1,k)$*/               
-	    a(1,-1,0) = - Uxy(i+1,j,k)-Uxy(i,j-1,k);
-	    /*$asw = uxy(i-1,j,k)+uxy(i,j-1,k)$*/               
-	    a(-1,-1,0) = Uxy(i-1,j,k)+Uxy(i,j-1,k);
-	    /*$ate = uxz(i,j,k+1)+uxz(i+1,j,k)$*/               
-	    a(1,0,1) = Uxz(i,j,k+1)+Uxz(i+1,j,k);
-	    /*$atw = - uxz(i-1,j,k)-uxz(i,j,k+1)$*/               
-	    a(-1,0,1) = - Uxz(i-1,j,k)-Uxz(i,j,k+1);
-	    /*$abe = - uxz(i+1,j,k)-uxz(i,j,k-1)$*/               
-	    a(1,0,-1) = - Uxz(i+1,j,k)-Uxz(i,j,k-1);
-	    /*$abw = uxz(i-1,j,k)+uxz(i,j,k-1)$*/               
-	    a(-1,0,-1) = Uxz(i-1,j,k)+Uxz(i,j,k-1);         
-	    /*$atn = uyz(i,j+1,k)+uyz(i,j,k+1)$*/               
-	    a(0,1,1) = Uyz(i,j+1,k)+Uyz(i,j,k+1);
-	    /*$ats = - uyz(i,j,k+1)-uyz(i,j-1,k)$*/               
-	    a(0,-1,1) = - Uyz(i,j,k+1)-Uyz(i,j-1,k);
-	    /*$abn = - uyz(i,j,k-1)-uyz(i,j+1,k)$*/               
-	    a(0,1,-1) = - Uyz(i,j,k-1)-Uyz(i,j+1,k);
-	    /*$asb = uyz(i,j,k-1)+uyz(i,j-1,k)$*/
-	    a(0,-1,-1) = Uyz(i,j,k-1)+Uyz(i,j-1,k);
-	  } /* end nabla_form=3 */
+          else { 
+            /* nabla_form==3: Upper Metric Nabla Form */
+            a(0,0,0) = -Uxx(i+1,j,k) -2.*Uxx(i,j,k) -Uxx(i-1,j,k) 
+              -Uyy(i,j+1,k) -2.*Uyy(i,j,k) -Uyy(i,j-1,k)
+              -Uzz(i,j,k+1) -2.*Uzz(i,j,k) -Uzz(i,j,k-1);
+            
+            /*$ae = uxx(i+1,j,k)+uxx(i,j,k)$*/
+            a(1,0,0) = Uxx(i+1,j,k) + Uxx(i,j,k);
+            /*$aw = uxx(i-1,j,k)+uxx(i,j,k)$*/
+            a(-1,0,0) = Uxx(i-1,j,k)+Uxx(i,j,k);
+                /*$an = uyy(i,j+1,k)+uyy(i,j,k)$*/              
+            a(0,1,0) = Uyy(i,j+1,k)+Uyy(i,j,k);
+            /*$as = uyy(i,j-1,k)+uyy(i,j,k)$*/              
+            a(0,-1,0) = Uyy(i,j-1,k)+Uyy(i,j,k);
+            /*$at = uzz(i,j,k+1)+uzz(i,j,k)$*/               
+            a(0,0,1) = Uzz(i,j,k+1)+Uzz(i,j,k);               
+            /*$ab = uzz(i,j,k-1)+uzz(i,j,k)$*/               
+            a(0,0,-1) = Uzz(i,j,k-1)+Uzz(i,j,k);
+            
+            /*$ane = uxy(i,j+1,k)+uxy(i+1,j,k)$*/              
+            a(1,1,0) = Uxy(i,j+1,k)+Uxy(i+1,j,k);
+            /*$anw = - uxy(i-1,j,k)-uxy(i,j+1,k)$*/               
+            a(-1,1,0) = - Uxy(i-1,j,k)-Uxy(i,j+1,k);
+            /*$ase = - uxy(i+1,j,k)-uxy(i,j-1,k)$*/               
+            a(1,-1,0) = - Uxy(i+1,j,k)-Uxy(i,j-1,k);
+            /*$asw = uxy(i-1,j,k)+uxy(i,j-1,k)$*/               
+            a(-1,-1,0) = Uxy(i-1,j,k)+Uxy(i,j-1,k);
+            /*$ate = uxz(i,j,k+1)+uxz(i+1,j,k)$*/               
+            a(1,0,1) = Uxz(i,j,k+1)+Uxz(i+1,j,k);
+            /*$atw = - uxz(i-1,j,k)-uxz(i,j,k+1)$*/               
+            a(-1,0,1) = - Uxz(i-1,j,k)-Uxz(i,j,k+1);
+            /*$abe = - uxz(i+1,j,k)-uxz(i,j,k-1)$*/               
+            a(1,0,-1) = - Uxz(i+1,j,k)-Uxz(i,j,k-1);
+            /*$abw = uxz(i-1,j,k)+uxz(i,j,k-1)$*/               
+            a(-1,0,-1) = Uxz(i-1,j,k)+Uxz(i,j,k-1);         
+            /*$atn = uyz(i,j+1,k)+uyz(i,j,k+1)$*/               
+            a(0,1,1) = Uyz(i,j+1,k)+Uyz(i,j,k+1);
+            /*$ats = - uyz(i,j,k+1)-uyz(i,j-1,k)$*/               
+            a(0,-1,1) = - Uyz(i,j,k+1)-Uyz(i,j-1,k);
+            /*$abn = - uyz(i,j,k-1)-uyz(i,j+1,k)$*/               
+            a(0,1,-1) = - Uyz(i,j,k-1)-Uyz(i,j+1,k);
+            /*$asb = uyz(i,j,k-1)+uyz(i,j-1,k)$*/
+            a(0,-1,-1) = Uyz(i,j,k-1)+Uyz(i,j-1,k);
+          } /* end nabla_form=3 */
 
           /* 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).  */
-	  VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
+          VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
           
-	  /* Put in the octant boundary conditions.
+          /* 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
@@ -789,12 +789,12 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
              
             */
 #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");
-	  }
+          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)
@@ -839,19 +839,19 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
                     /* This is a boundary. */
                     rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)];
                   } 
-		  else {
+                  else {
                     J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRQ(ierr);
+                    CHKERRQ(ierr);
                   }
                 }
               }
 
-	  if (Nstorage)
-	    rhsval = -rhsval - Nlin[ijk] / ac;
+          if (Nstorage)
+            rhsval = -rhsval - Nlin[ijk] / ac;
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRQ(ierr);
+          CHKERRQ(ierr);
           
         }
       }
@@ -869,7 +869,7 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 
   if (nabla_form == 3) {
     if (verbose)
-	printf ("Freeing upper metric storage\n");
+        printf ("Freeing upper metric storage\n");
     free(uxx3);
     free(uxy3);
     free(uxz3);
diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c
index dc09aa3..645502a 100644
--- a/src/petsc_flat_solver.c
+++ b/src/petsc_flat_solver.c
@@ -64,10 +64,10 @@ static Vec     soln, b;      /* approx solution, RHS */
 static SLES    sles;         /* linear solver context */
 
 int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 
-		int *AbsTol, int *RelTol);
+                int *AbsTol, int *RelTol);
 
 int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 
-		int *AbsTol, int *RelTol) 
+                int *AbsTol, int *RelTol) 
 {
 
   DECLARE_CCTK_PARAMETERS    /* CCTK passed parameters */
@@ -113,7 +113,7 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   int octant;           /* Apply octant BCs inside */
   int matnormalize;     /* Normalize the central mat value to one? */
   int I,J;              /* The PETSc col and row in the matrix */
-	    
+            
   int PetscTolStyle;
 
   
@@ -206,19 +206,19 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   /* We fix the lower and upper boundary indices, that actually "active" in 
      the sense that they are no ghostzones: */
   imin=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
-	 GH->cctk_nghostzones[0]);
+         GH->cctk_nghostzones[0]);
   imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : 
-	 pEx->lnsize[0]-GH->cctk_nghostzones[0]);
+         pEx->lnsize[0]-GH->cctk_nghostzones[0]);
 
   jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
-	GH->cctk_nghostzones[1]);
+        GH->cctk_nghostzones[1]);
   jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : 
-	 pEx->lnsize[1]-GH->cctk_nghostzones[1]);
+         pEx->lnsize[1]-GH->cctk_nghostzones[1]);
   
   kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
-	 GH->cctk_nghostzones[2]);
+         GH->cctk_nghostzones[2]);
   kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : 
-	 pEx->lnsize[2]-GH->cctk_nghostzones[2]);
+         pEx->lnsize[2]-GH->cctk_nghostzones[2]);
 
   /* We need to get the global index of gridpoints (gps) owned on my proc. 
      For that we have to know how many gps are there before my processors
@@ -317,12 +317,12 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
  
   if (verbose) CCTK_INFO("Creating Matrices and Vectors ....");
   ierr = MatCreateMPIAIJ(pughGH->PUGH_COMM_WORLD,     
-			 (endpoint-startpoint),   /* # of rows */
-			 (endpoint-startpoint),   /* This is confusing */
-			 rank,rank,               /* Matrix size */
-			 19,PETSC_NULL,           /* Diagonal  */
-			 19,PETSC_NULL,           /* Off diagonal */
-			 &A[0]);                  /* The output  */
+                         (endpoint-startpoint),   /* # of rows */
+                         (endpoint-startpoint),   /* This is confusing */
+                         rank,rank,               /* Matrix size */
+                         19,PETSC_NULL,           /* Diagonal  */
+                         19,PETSC_NULL,           /* Off diagonal */
+                         &A[0]);                  /* The output  */
   CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
   CHKERRQ(ierr);
@@ -374,34 +374,34 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
           ig  = i + GH->cctk_lbnd[0]; 
           jg  = j + GH->cctk_lbnd[1];
           kg  = k + GH->cctk_lbnd[2];
-	
+        
           /* Get the row we are working on */
           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];
+          }
 
           /* Now set the values of the matrix. Here we force dirichlet. */
-	  VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
+          VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
           
 #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");
-	  }
+          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)
@@ -425,12 +425,12 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
       
           /* renormalize */
           if (matnormalize) 
-	  {
+          {
             ac = a(0,0,0);
             for (J=0;J<27;J++) a[J] = a[J] / ac;
           } 
           else 
-	    ac = 1;
+            ac = 1;
           
 
           /* This is the new way-clever look. Note it relies
@@ -439,40 +439,40 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
              19-way banded matrix. */
 
           for (l=-1;l<=1;l++)
-	  {
+          {
             for (m=-1;m<=1;m++)
-	    {
+            {
               for (n=-1;n<=1;n++) 
-	      {
+              {
                 if (l*m*n == 0) 
-		{       /* This is the 19 point ... if none are
-			 * zero, then we have no stencil here.
-			 */
+                {       /* This is the 19 point ... if none are
+                         * zero, then we have no stencil here.
+                         */
                   if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) 
-		  {
+                  {
                     /* This is a boundary. */
                     rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)];
                   } 
-		  else 
-		  {
+                  else 
+                  {
                     J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRQ(ierr);
+                    CHKERRQ(ierr);
                   }
                 }
               }
-	    }
-	  }
-	  if (Nstorage) 
-	  {
-	    rhsval = -rhsval - Nlin[ijk] / ac;
-	  }
+            }
+          }
+          if (Nstorage) 
+          {
+            rhsval = -rhsval - Nlin[ijk] / ac;
+          }
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRQ(ierr);
+          CHKERRQ(ierr);
           
-	  /* reset the central value*/
-	  a( 0, 0, 0) = -6.0;
+          /* reset the central value*/
+          a( 0, 0, 0) = -6.0;
   
         }   /* if (wsp>=0) */
       }   /*for i */
@@ -667,7 +667,7 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
         int ijk;                /* The data point for the array */
         int I;                  /* The col in the matrix */
 
-	/* these guys are easy */
+        /* these guys are easy */
         ijk = DATINDEX(pEx,i,j,k);
         if (wsp[ijk] >= 0) 
         {
@@ -707,5 +707,3 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   
   return (0);
 }
-
-
diff --git a/src/petsc_wrapper.c b/src/petsc_wrapper.c
index 8d4fbcf..0f6fb65 100644
--- a/src/petsc_wrapper.c
+++ b/src/petsc_wrapper.c
@@ -8,8 +8,8 @@ static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_wrapper_c)
 
 void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex, 
-		      int MIndex, int NIndex, 
-		      CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+                      int MIndex, int NIndex, 
+                      CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
 
 /* The wrapper functions for the core PETSc solver, that performs the 
    actual solve. One wrapper function for each class, because these 
@@ -18,14 +18,14 @@ void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex,
 
 /* Wrapper function for the class of elliptic equations that needs a metric */
 void petsc_metric(cGH *GH, int *MetricI, int FieldIndex, 
-		  int MIndex, int NIndex, 
-		  CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
+                  int MIndex, int NIndex, 
+                  CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
 
 /* Arrays for the M/N size info */
   
   int petsc_confmetric_solver(cGH *GH, int *MetricI, int MetricISize, 
-			       int FieldIndex, int MIndex, int NIndex, 
-			       CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
+                               int FieldIndex, int MIndex, int NIndex, 
+                               CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
   
   
   if (GH->cctk_dim>3) 
@@ -34,22 +34,22 @@ void petsc_metric(cGH *GH, int *MetricI, int FieldIndex,
   }
 
   petsc_confmetric_solver(GH, MetricI, 6, FieldIndex, 
-			  MIndex, NIndex, 
-			  AbsTol, RelTol);
+                          MIndex, NIndex, 
+                          AbsTol, RelTol);
 }
   
 /* wrapper function for the class of elliptic equations, that needs a conf.
    factor */
 void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex, 
-		      int MIndex, int NIndex, 
-		      CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
+                      int MIndex, int NIndex, 
+                      CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
   
   int petsc_confmetric_solver(cGH *, int *, int, int, int, int, 
-			       CCTK_REAL *, CCTK_REAL *);
+                               CCTK_REAL *, CCTK_REAL *);
 
  /* petsc_confmetric_solver expects an metric array, in the case of this 
     equation class,  needs conf.factor as last entry -> size 7 */
   petsc_confmetric_solver(GH, MetricPsiI, 7,
-			  FieldIndex, MIndex, NIndex, 
-			  AbsTol, RelTol);
+                          FieldIndex, MIndex, NIndex, 
+                          AbsTol, RelTol);
 }