Changes for new PUGH ... compiles but untested

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

3 files changed, 95 insertions, 93 deletions

diff --git a/src/petsc_confmetric.c b/src/petsc_confmetric.c
index ea78b3c..6dd03d3 100644
--- a/src/petsc_confmetric.c
+++ b/src/petsc_confmetric.c
@@ -11,10 +11,7 @@
 
 #include "cctk.h"
 
-#include "cctk_CommandLine.h"
-
 #include "cctk_Parameters.h"
-#include "cctk_Flesh.h"
 
 #include "assert.h"
 #include "mpi.h"
diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c
index 54b4002..b439c8f 100644
--- a/src/petsc_confmetric_solver.c
+++ b/src/petsc_confmetric_solver.c
@@ -11,11 +11,7 @@
 
 
 #include "cctk.h"
-
-#include "cctk_CommandLine.h"
-
 #include "cctk_Parameters.h"
-#include "cctk_Flesh.h"
 
 #include "assert.h"
 #include "mpi.h"
@@ -38,12 +34,12 @@ static char *rcsid = "$Header$";
 
   
 /* Some useful definitions to deal with the upper metric */
-#define Uxx(i,j,k) uxx3[DATINDEX(pughGH,(i),(j),(k))]
-#define Uxy(i,j,k) uxy3[DATINDEX(pughGH,(i),(j),(k))]
-#define Uxz(i,j,k) uxz3[DATINDEX(pughGH,(i),(j),(k))]
-#define Uyy(i,j,k) uyy3[DATINDEX(pughGH,(i),(j),(k))]
-#define Uyz(i,j,k) uyz3[DATINDEX(pughGH,(i),(j),(k))]
-#define Uzz(i,j,k) uzz3[DATINDEX(pughGH,(i),(j),(k))]
+#define Uxx(i,j,k) uxx3[DATINDEX(pEx,(i),(j),(k))]
+#define Uxy(i,j,k) uxy3[DATINDEX(pEx,(i),(j),(k))]
+#define Uxz(i,j,k) uxz3[DATINDEX(pEx,(i),(j),(k))]
+#define Uyy(i,j,k) uyy3[DATINDEX(pEx,(i),(j),(k))]
+#define Uyz(i,j,k) uyz3[DATINDEX(pEx,(i),(j),(k))]
+#define Uzz(i,j,k) uzz3[DATINDEX(pEx,(i),(j),(k))]
   
 
 /* Place the matrix in global space so we can keep it around
@@ -104,6 +100,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   int imin,imax,jmin,jmax,kmin,kmax,interior;  
    
   pGH *pughGH;                 /* The pugh Extension handle */
+  pGExtras *pEx;
+  pConnectivity *pCon;
   int myproc;                  /* out processor */
 
   /* Tolerances */
@@ -176,6 +174,10 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   pughGH = (pGH*)GH->extensions[CCTK_GHExtensionHandle("PUGH")];
   if (!pughGH) CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle\n");
 
+  /* Get the extras extension for 3D grid functions */
+  pEx = pughGH->GFExtras[2];
+  pCon = pughGH->Connectivity[2];
+
   /* Things to do on first iteration */
   if (trips==0) {
     int argc;
@@ -240,7 +242,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   
   /* initialize the linear index lookup table, after it's 
      filled up (below), the -1 indicates a boundary */
-  for (i=0;i<pughGH->npoints;i++) wsp[i] = -1.0;
+  for (i=0;i<pEx->npoints;i++) wsp[i] = -1.0;
 
   /* Get myproc from the CCTK, get the gridspacing */
   myproc = pugh_MyProc(GH);
@@ -250,20 +252,20 @@ void 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=((pughGH->neighbors[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
+  imin=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
 	 GH->cctk_nghostzones[0]);
-  imax=((pughGH->neighbors[pughGH->myproc][XDP]<0) ? pughGH->lnsize[0]-1 : 
-	 pughGH->lnsize[0]-GH->cctk_nghostzones[0]);
+  imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : 
+	 pEx->lnsize[0]-GH->cctk_nghostzones[0]);
 
-  jmin=((pughGH->neighbors[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
+  jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
 	GH->cctk_nghostzones[1]);
-  jmax=((pughGH->neighbors[pughGH->myproc][YDP]<0) ? pughGH->lnsize[1]-1 : 
-	 pughGH->lnsize[1]-GH->cctk_nghostzones[1]);
+  jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : 
+	 pEx->lnsize[1]-GH->cctk_nghostzones[1]);
   
-  kmin=((pughGH->neighbors[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
+  kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
 	 GH->cctk_nghostzones[2]);
-  kmax=((pughGH->neighbors[pughGH->myproc][ZDP]<0) ? pughGH->lnsize[2]-1 : 
-	 pughGH->lnsize[2]-GH->cctk_nghostzones[2]);
+  kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : 
+	 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
@@ -273,20 +275,20 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 
   for (i=0;i<myproc;i++) {
     
-    nxs=pughGH->rnsize[i][0];
-    nxs=((pughGH->neighbors[i][XDM]<0) && !(octant) ? nxs-1 : 
+    nxs=pEx->rnsize[i][0];
+    nxs=((pCon->neighbours[i][XDM]<0) && !(octant) ? nxs-1 : 
          nxs-GH->cctk_nghostzones[0]); 
-    nxs=((pughGH->neighbors[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]);
+    nxs=((pCon->neighbours[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]);
     
-    nys=pughGH->rnsize[i][1];
-    nys=((pughGH->neighbors[i][YDM]<0) && !(octant) ? nys-1 : 
+    nys=pEx->rnsize[i][1];
+    nys=((pCon->neighbours[i][YDM]<0) && !(octant) ? nys-1 : 
          nys-GH->cctk_nghostzones[1]);
-    nys=((pughGH->neighbors[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]);
+    nys=((pCon->neighbours[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]);
     
-    nzs=pughGH->rnsize[i][2];
-    nzs=((pughGH->neighbors[i][ZDM]<0) && !(octant) ? nzs-1 : 
+    nzs=pEx->rnsize[i][2];
+    nzs=((pCon->neighbours[i][ZDM]<0) && !(octant) ? nzs-1 : 
          nzs-GH->cctk_nghostzones[2]);
-    nzs=((pughGH->neighbors[i][ZDP]<0) ? nzs-1 : nzs-GH->cctk_nghostzones[2]);    
+    nzs=((pCon->neighbours[i][ZDP]<0) ? nzs-1 : nzs-GH->cctk_nghostzones[2]);    
     
     printf("PROC: %d  nxyzs %d %d %d \n",i,nxs,nys,nzs);
     startpoint += nxs*nys*nzs;
@@ -303,11 +305,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   for (k=kmin;k<kmax;k++)
     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
-
+        wsp[DATINDEX(pEx,i,j,k)] = endpoint++;
+  
   /* 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
@@ -320,14 +319,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
      as well as the column for the stencil of the neighbors
      (workspace->data[DATINDEX(GH,i+1,j,k)] etc...). So onwards */
 
-  nxs = pughGH->nsize[0]-2;
-  nys = pughGH->nsize[1]-2;
-  nzs = pughGH->nsize[2]-2;
+  nxs = pEx->nsize[0]-2;
+  nys = pEx->nsize[1]-2;
+  nzs = pEx->nsize[2]-2;
 
 #ifdef DEBUG
   if (debug) {
     printf("nxyzs %d %d %d \n",nxs,nys,nzs);
-    printf("lnxyz: %d %d %d \n",pughGH->lnsize[0],pughGH->lnsize[1],pughGH->lnsize[2]);
+    printf("lnxyz: %d %d %d \n",pEx->lnsize[0],pEx->lnsize[1],pEx->lnsize[2]);
   }
 #endif
 
@@ -402,13 +401,13 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
       if (verbose) 
 	CCTK_INFO("Forming nabla with upper g and finite difference of dg \n");
       nabla_form = 3;
-      uxx3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      uxy3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      uxz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      uyy3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      uyz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      uzz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL));
-      for (i=0;i<pughGH->npoints;i++) {
+      uxx3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
+      uxy3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
+      uxz3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
+      uyy3 = (CCTK_REAL*)malloc(pEx->npoints*sizeof(CCTK_REAL));
+      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]) + 
@@ -458,7 +457,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     for (j=jmin;j<jmax;j++) {
       for (i=imin;i<imax;i++) {
 
-        if (wsp[DATINDEX(pughGH,i,j,k)] >= 0) {
+        if (wsp[DATINDEX(pEx,i,j,k)] >= 0) {
 
           CCTK_REAL tdxgxx, tdxgxy, tdxgxz, tdxgyy, tdxgyz, tdxgzz;
           CCTK_REAL tdygxx, tdygxy, tdygxz, tdygyy, tdygyz, tdygzz;
@@ -473,7 +472,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
           for (I=0;I<27;I++) a[I] = 0.0;
           
           /* these guys are easy */
-          ijk = DATINDEX(pughGH,i,j,k);               /* get linear index */
+          ijk = DATINDEX(pEx,i,j,k);               /* get linear index */
           ig  = i + GH->cctk_lbnd[0]; 
           jg  = j + GH->cctk_lbnd[1];
           kg  = k + GH->cctk_lbnd[2];
@@ -498,8 +497,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
             int ijkp, ijkm;
             
             /* X derivatives */
-            ijkp = DATINDEX(pughGH,i+1,j,k);
-            ijkm = DATINDEX(pughGH,i-1,j,k);
+            ijkp = DATINDEX(pEx,i+1,j,k);
+            ijkm = DATINDEX(pEx,i-1,j,k);
             
             tdxgxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dx);
             tdxgxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dx);
@@ -510,8 +509,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
               
               
             /* Y derivatives */
-            ijkp = DATINDEX(pughGH,i,j+1,k);
-            ijkm = DATINDEX(pughGH,i,j-1,k);
+            ijkp = DATINDEX(pEx,i,j+1,k);
+            ijkm = DATINDEX(pEx,i,j-1,k);
             
             tdygxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dy);
             tdygxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dy);
@@ -521,8 +520,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
             tdygzz = (gzz[ijkp] - gzz[ijkm])/(2.0*dy);
             
             /* X derivatives */
-            ijkp = DATINDEX(pughGH,i,j,k+1);
-            ijkm = DATINDEX(pughGH,i,j,k-1);
+            ijkp = DATINDEX(pEx,i,j,k+1);
+            ijkm = DATINDEX(pEx,i,j,k-1);
             
             tdzgxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dz);
             tdzgxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dz);
@@ -781,7 +780,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
               for (m=-1;m<=0;m++)
                 for (n=-1;n<=0;n++)
                   if (l*m*n == 0) 
-                    if (wsp[DATINDEX(pughGH,i+l,j+m,k+n)] < 0 &&
+                    if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0 &&
                         (ig == imin || jg == jmin || kg == kmin))  
                       {
                         /* We are on an inner boundary point, eg, 
@@ -814,12 +813,12 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
                 if (l*m*n == 0) {       /* This is the 19 point ... if none are
                                          * zero, then we have no stencil here.
                                          */
-                  if (wsp[DATINDEX(pughGH,i+l,j+m,k+n)] < 0) {
+                  if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) {
                     /* This is a boundary. */
-                    rhsval += a(l,m,n) * ell_field[DATINDEX(pughGH,i+l,j+m,k+n)];
+                    rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)];
                   } 
 		  else {
-                    J = wsp[DATINDEX(pughGH,i+l,j+m,k+n)];
+                    J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
 		    CHKERRA(ierr);
                   }
@@ -1033,7 +1032,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
         CCTK_REAL rhsval = 0;
 
         /* these guys are easy */
-        ijk = DATINDEX(pughGH,i,j,k);
+        ijk = DATINDEX(pEx,i,j,k);
         if (wsp[ijk] >= 0) {
           ig = i + GH->cctk_lbnd[0];
           jg = j + GH->cctk_lbnd[1];
diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c
index 662c81b..2607937 100644
--- a/src/petsc_flat_solver.c
+++ b/src/petsc_flat_solver.c
@@ -66,7 +66,9 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   int imin,imax,jmin,jmax,kmin,kmax,interior;  
    
   pGH *pughGH;          /* The pugh Extension handle */
- 
+  pGExtras *pEx;
+  pConnectivity *pCon;
+
   /* Tolerances */
   CCTK_REAL rtol, atol, tolerance;
 
@@ -124,6 +126,10 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   pughGH = (pGH*)GH->extensions[CCTK_GHExtensionHandle("PUGH")];
   if (!pughGH) CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle\n");
 
+  /* Get the extras extension for 3D grid functions */
+  pEx = pughGH->GFExtras[2];
+  pCon = pughGH->Connectivity[2];
+
   /* Things to do on first iteration */
   if (trips==0) {
     int argc;
@@ -167,7 +173,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 
   /* 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;
+  for (i=0;i<pEx->npoints;i++) wsp[i] = -1.0;
 
   /* Get myproc from the CCTK, get the gridspacing */
   myproc = pugh_MyProc(GH);
@@ -177,20 +183,20 @@ void 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=((pughGH->neighbors[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
+  imin=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : 
 	 GH->cctk_nghostzones[0]);
-  imax=((pughGH->neighbors[pughGH->myproc][XDP]<0) ? pughGH->lnsize[0]-1 : 
-	 pughGH->lnsize[0]-GH->cctk_nghostzones[0]);
+  imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : 
+	 pEx->lnsize[0]-GH->cctk_nghostzones[0]);
 
-  jmin=((pughGH->neighbors[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
+  jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
 	GH->cctk_nghostzones[1]);
-  jmax=((pughGH->neighbors[pughGH->myproc][YDP]<0) ? pughGH->lnsize[1]-1 : 
-	 pughGH->lnsize[1]-GH->cctk_nghostzones[1]);
+  jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : 
+	 pEx->lnsize[1]-GH->cctk_nghostzones[1]);
   
-  kmin=((pughGH->neighbors[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
+  kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
 	 GH->cctk_nghostzones[2]);
-  kmax=((pughGH->neighbors[pughGH->myproc][ZDP]<0) ? pughGH->lnsize[2]-1 : 
-	 pughGH->lnsize[2]-GH->cctk_nghostzones[2]);
+  kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : 
+	 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
@@ -200,20 +206,20 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 
   for (i=0;i<myproc;i++) {
     
-    nxs=pughGH->rnsize[i][0];
-    nxs=((pughGH->neighbors[i][XDM]<0) && !(octant) ? nxs-1 : 
+    nxs=pEx->rnsize[i][0];
+    nxs=((pCon->neighbours[i][XDM]<0) && !(octant) ? nxs-1 : 
          nxs-GH->cctk_nghostzones[0]); 
-    nxs=((pughGH->neighbors[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]);
+    nxs=((pCon->neighbours[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]);
     
-    nys=pughGH->rnsize[i][1];
-    nys=((pughGH->neighbors[i][YDM]<0) && !(octant) ? nys-1 : 
+    nys=pEx->rnsize[i][1];
+    nys=((pCon->neighbours[i][YDM]<0) && !(octant) ? nys-1 : 
          nys-GH->cctk_nghostzones[1]);
-    nys=((pughGH->neighbors[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]);
+    nys=((pCon->neighbours[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]);
     
-    nzs=pughGH->rnsize[i][2];
-    nzs=((pughGH->neighbors[i][ZDM]<0) && !(octant) ? nzs-1 : 
+    nzs=pEx->rnsize[i][2];
+    nzs=((pCon->neighbours[i][ZDM]<0) && !(octant) ? nzs-1 : 
          nzs-GH->cctk_nghostzones[2]);
-    nzs=((pughGH->neighbors[i][ZDP]<0) ? nzs-1 : nzs-GH->cctk_nghostzones[2]);    
+    nzs=((pCon->neighbours[i][ZDP]<0) ? nzs-1 : nzs-GH->cctk_nghostzones[2]);    
     
     printf("PROC: %d  nxyzs %d %d %d \n",i,nxs,nys,nzs);
     startpoint += nxs*nys*nzs;
@@ -230,7 +236,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   for (k=kmin;k<kmax;k++)
     for (j=jmin;j<jmax;j++)
       for (i=imin;i<imax;i++)
-        wsp[DATINDEX(pughGH,i,j,k)] = endpoint++;
+        wsp[DATINDEX(pEx,i,j,k)] = endpoint++;
 
 #ifdef DEBUG  
   printf("ENDPOINT:  %d \n",endpoint);
@@ -247,14 +253,14 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      as well as the column for the stencil of the neighbors
      (workspace->data[DATINDEX(GH,i+1,j,k)] etc...). So onwards */
 
-  nxs = pughGH->nsize[0]-2;
-  nys = pughGH->nsize[1]-2;
-  nzs = pughGH->nsize[2]-2;
+  nxs = pEx->nsize[0]-2;
+  nys = pEx->nsize[1]-2;
+  nzs = pEx->nsize[2]-2;
 
 #ifdef DEBUG
   if (debug) {
     printf("nxyzs %d %d %d \n",nxs,nys,nzs);
-    printf("lnxyz: %d %d %d \n",pughGH->lnsize[0],pughGH->lnsize[1],pughGH->lnsize[2]);
+    printf("lnxyz: %d %d %d \n",pEx->lnsize[0],pEx->lnsize[1],pEx->lnsize[2]);
   }
 #endif
 
@@ -332,7 +338,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     for (j=jmin;j<jmax;j++) {
       for (i=imin;i<imax;i++) {
 
-        if (wsp[DATINDEX(pughGH,i,j,k)] >= 0) {
+        if (wsp[DATINDEX(pEx,i,j,k)] >= 0) {
 
           /* Set up indices */
           int ijk;              /* The lineax index for the array */
@@ -340,7 +346,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
           CCTK_REAL rhsval = 0;
 
           /* local linear index / global 3d index */
-          ijk = DATINDEX(pughGH,i,j,k);    
+          ijk = DATINDEX(pEx,i,j,k);    
           ig  = i + GH->cctk_lbnd[0]; 
           jg  = j + GH->cctk_lbnd[1];
           kg  = k + GH->cctk_lbnd[2];
@@ -368,7 +374,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
               for (m=-1;m<=0;m++)
                 for (n=-1;n<=0;n++)
                   if (l*m*n == 0) 
-                    if (wsp[DATINDEX(pughGH,i+l,j+m,k+n)] < 0 &&
+                    if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0 &&
                         (ig == imin || jg == jmin || kg == kmin))  
                       {
                         /* We are on an inner boundary point, eg, 
@@ -402,12 +408,12 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
                 if (l*m*n == 0) {       /* This is the 19 point ... if none are
                                          * zero, then we have no stencil here.
                                          */
-                  if (wsp[DATINDEX(pughGH,i+l,j+m,k+n)] < 0) {
+                  if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) {
                     /* This is a boundary. */
-                    rhsval += a(l,m,n) * ell_field[DATINDEX(pughGH,i+l,j+m,k+n)];
+                    rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)];
                   } 
 		  else {
-                    J = wsp[DATINDEX(pughGH,i+l,j+m,k+n)];
+                    J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
 		    CHKERRA(ierr);
                   }
@@ -577,7 +583,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
         CCTK_REAL rhsval = 0;
 
         /* these guys are easy */
-        ijk = DATINDEX(pughGH,i,j,k);
+        ijk = DATINDEX(pEx,i,j,k);
         if (wsp[ijk] >= 0) {
           ig = i + GH->cctk_lbnd[0];
           jg = j + GH->cctk_lbnd[1];