Remove unused file.

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

1 files changed, 0 insertions, 942 deletions

diff --git a/src/petsc_confmetric.c b/src/petsc_confmetric.c
deleted file mode 100644
index 34e1ced..0000000
--- a/src/petsc_confmetric.c
+++ /dev/null
@@ -1,942 +0,0 @@
-/*@@
-   @file      PETScEll_conformal.c
-   @date      Wed Apr  9 09:31:24 1997
-   @author    Paul Walker
-   @desc
-   A nabla phi - M phi = N elliptic solver based around PETSc.
-   For information, see the documentation for @seeroutine petscEll_conformal
-   @enddesc
-   @version $Header$
-@@*/
-
-#include "cctk.h"
-
-#include "cctk_Parameters.h"
-
-#include "assert.h"
-#include "mpi.h"
-#include "petscsles.h"
-#include "ellpetsc.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_confmetric_c)
-
-#define DEBUG
-
-/*Don't know what these actually mean! FIXME */
-#define ELLCONV_ABSOLUTE 0
-#define ELLCONV_RELATIVE 1
-#define ELLCONV_EITHER   2
-
-/* A few convenient macros */
-#define a(i,j,k) a[i+1 + 3*(j+1 + 3*(k+1))]
-#define SQR(a) ((a)*(a))
-
-  
-/* 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))]
-  
-
-/* 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 *GetDataPtr_NextTL(cGH *GH, const char *field) {
-  int index, tlevel;
-  char *err;
-
-  index = CCTK_VarIndex(field);
-  if (index<0) {
-    err = (char*) malloc( 256*sizeof(char));
-    sprintf(err,"ERROR: Index for >%s< not found",field);
-    CCTK_WARN(1,err);
-    if (err) free(err);
-  }
-  return((CCTK_REAL*)CCTK_VarDataPtrI(GH,0,index));
-}
-
-
-/* This passing matches that of in the convention in the
-   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) {
-
-  DECLARE_CCTK_PARAMETERS    /* CCTK passed parameters */
-
-
-  void *GetDataPtr_NextTL(cGH *GH, const char *field);
-
-  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 */
-  
-  /* Loopers */
-  int i,j,k,l,m,n;
-
-  /* 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 */
-
-  /* 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 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 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 */
-  int pstart, pend;     /* A check for PETSc layout */
-  int pvstart, pvend;   /* A check for PETSc layout */
-  int verbose;          /* Is the solver verbose */
-  int debug;            /* Is the solver debug-verbose */
-  int octant;           /* Apply octant BCs inside */
-  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 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) */
-  CCTK_REAL *wsp =NULL, *ell_field=NULL; 
-  CCTK_REAL *gxx =NULL, *gxy =NULL; 
-  CCTK_REAL *gxz =NULL, *gyy =NULL; 
-  CCTK_REAL *gyz =NULL, *gzz =NULL; 
-  CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
-  CCTK_REAL *psi =NULL;
-  CCTK_REAL *psix=NULL, *psiy=NULL, *psiz=NULL; 
-
-  /* Get some parameters */
-  octant       = CCTK_Equals(domain,"octant");
-  verbose      = CCTK_Equals(petsc_verbose,"yes")||CCTK_Equals(petsc_verbose,"debug");
-  debug        = CCTK_Equals(petsc_verbose,"debug");
-  reuse        = 0; /*$CCTK_Equals(petsc_reuse,"yes");$*/
-  matnormalize = 0; /*$CCTK_Equals(petsc_coeff_to_one,"yes");$*/ 
-
-  /* FIXME, the TolAbs/TolRel will be evaluated here */
-  PetscTolStyle=0;
-  tolerance    =AbsTol[0];
-
-  printf("\n**** PETSC  !!!!****\n\n");
-
-
-  /* Get the link to pugh Extension */
-  pughGH = (pGH*)GH->extensions[CCTK_GHExtensionHandle("PUGH")];
-  if (!pughGH) CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle\n");
-
-  /* Things to do on first iteration */
-  if (trips==0) {
-    int argc;
-    char **argv; 
-
-#ifdef DEBUG
-    if (debug)
-      printf("PETSc: initial trip: %d \n",trips);
-#endif
-
-    /* Get the commandline arguments */
-    argc = CCTK_CommandLine(&argv);
-
-    /* Set the PETSc communicator to set of PUGH and
-       initialzie PETSc */
-    ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRA(ierr);
-    PetscInitialize(&argc,&argv,NULL,NULL); 
-
-    CCTK_INFO("PETSc initialized");
-  }
-
-  trips++;
-
-  /* Create a array of matrices A */
-  /* num_A = MultiGridCount(); */
-  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);
-  Mlin      = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*MIndex);
-  Nlin      = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex);
-
-  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]);
-
-  /* FIXME: evolved derivatives should go! -> Einstein specific*/  
-  if (conformal) {
-    psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[6]);
-    psix =(CCTK_REAL*) GetDataPtr_NextTL(GH,"einstein::psix");
-    psiy =(CCTK_REAL*) GetDataPtr_NextTL(GH,"einstein::psiy");
-    psiz =(CCTK_REAL*) GetDataPtr_NextTL(GH,"einstein::psiz");
-  }
-
-  wsp  = GetDataPtr_NextTL(GH,"ellpetsc::wsp");
-
-  
-  /* 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;
-
-  /* Get myproc from the CCTK, get the gridspacing */
-  myproc = PUGH_MyProc(GH);
-  dx     = GH->cctk_delta_space[0];
-  dy     = GH->cctk_delta_space[1];
-  dz     = GH->cctk_delta_space[2];
-
-  /* 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 : 
-	 pughGH->nghostzones);
-  imax=((pughGH->neighbors[pughGH->myproc][XDP]<0) ? pughGH->lnx-1 : 
-	 pughGH->lnx-GH->cctk_nghostzones[0]);
-
-  jmin=((pughGH->neighbors[pughGH->myproc][YDM]<0) && !(octant) ? 1 : 
-	GH->cctk_nghostzones[1]);
-  jmax=((pughGH->neighbors[pughGH->myproc][YDP]<0) ? pughGH->lny-1 : 
-	 pughGH->lny-GH->cctk_nghostzones[1]);
-  
-  kmin=((pughGH->neighbors[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 
-	 GH->cctk_nghostzones[2]);
-  kmax=((pughGH->neighbors[pughGH->myproc][ZDP]<0) ? pughGH->lnz-1 : 
-	 pughGH->lnz-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 (i<myproc), 
-     so we count them here. (that's minus the ghostszones (we use PUGH's "rn" 
-     variables (number of gp on proc[#] in direction [i]) and subtract the
-     ghostzones if nec.) */
-
-  for (i=0;i<myproc;i++) {
-    
-    nxs=pughGH->rnx[i];
-    nxs=((pughGH->neighbors[i][XDM]<0) && !(octant) ? nxs-1 : 
-         nxs-GH->cctk_nghostzones[0]); 
-    nxs=((pughGH->neighbors[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]);
-    
-    nys=pughGH->rny[i];
-    nys=((pughGH->neighbors[i][YDM]<0) && !(octant) ? nys-1 : 
-         nys-GH->cctk_nghostzones[1]);
-    nys=((pughGH->neighbors[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]);
-    
-    nzs=pughGH->rnz[i];
-    nzs=((pughGH->neighbors[i][ZDM]<0) && !(octant) ? nzs-1 : 
-         nzs-GH->cctk_nghostzones[2]);
-    nzs=((pughGH->neighbors[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;
-  }
-
-#ifdef DEBUG 
-  printf("STARTPOINT: %d \n",startpoint);
-#endif 
-
-  /* ...we save that as a start ...*/
-  endpoint = startpoint;
-    
-  /* ...and continue running over our own region to get our endpoint */
-  for (k=imin;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
-
-  /* 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 */
-  retcode = PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
-  
-  
-  /* So woohoo. Now for each point in our ijk space, we have
-     information about our row in the matrix (workspace->data[ijk])
-     as well as the column for the stencil of the neighbors
-     (workspace->data[DATINDEX(GH,i+1,j,k)] etc...). So onwards */
-
-  /* FIXME I don't think this works for stencil=2 ! */
-#ifdef DEBUG
-  if (debug)
-    CCTK_DumpGHInfo(GH);
-#endif
-
-  nxs = pughGH->nx-2;
-  nys = pughGH->ny-2;
-  nzs = pughGH->nz-2;
-
-#ifdef DEBUG
-  if (debug) {
-    printf("nxyzs %d %d %d \n",nxs,nys,nzs);
-    printf("lnxyz: %d %d %d \n",pughGH->lnx,pughGH->lny,pughGH->lnz);
-  }
-#endif
-
-  /* Suck off boundaries... */
-  rank = nxs*nys*nzs;
-  
-  
-  /* Create the petsc matrix A and vectors x and b efficiently.
-     
-     Now for efficiency we have to be really careful here. By
-     setting up the index array, we are abandoning the 19-way
-     tri-diagonal system so we can lay out our matrices on the
-     processors with the data. (One one proc we will still
-     get 19 way tridiagonal). So we want each processor
-     to own only part of the matrix/vector but we know exactly
-     what it is, eg, from startpoint to endpoint-1 is on our
-     local processor (that is, row-wise ; all columns are local).
-
-     So we can use MatCreateMPIAIJ for exactly this, and
-     VecCreateMPI for the vectors.
-   */
-
-  /* FIXME: perhaps only set up only on first iteration, this was a hack
-     by PAUL, using some inofficical petsc code, check if this is 
-     "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 */
-			 (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);
-  ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
-  CHKERRQ(ierr);
-  /*$}$*/ 
-  /*$else {$*/
-  /*  ierr = MatRestoreValuesMemory(A);  */
-  /*$CCTK_WARN(0,"no reuse");
-    CHKERRQ(ierr);
-  }$*/
-
-
-  /* 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 (pstart != startpoint && pend != endpoint) 
-    CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why??)\n");
-
-  nabla_form = 2;
-  if (verbose)
-    CCTK_INFO("Forming nabla with lower g and finite difference of g");
-
-  if (verbose) 
-    CCTK_INFO("Creating the coefficient matrix");
-
-  if (verbose && matnormalize) 
-    CCTK_INFO ("With diagonal renormalized to one");
-
-  for (k=kmin;k<kmax;k++) {
-    for (j=jmin;j<jmax;j++) {
-      for (i=imin;i<imax;i++) {
-
-        if (wsp[DATINDEX(pughGH,i,j,k)] >= 0) { /* FIXME */
-	  
-          CCTK_REAL tdxgxx, tdxgxy, tdxgxz, tdxgyy, tdxgyz, tdxgzz;
-          CCTK_REAL tdygxx, tdygxy, tdygxz, tdygyy, tdygyz, tdygzz;
-          CCTK_REAL tdzgxx, tdzgxy, tdzgxz, tdzgyy, tdzgyz, tdzgzz;
-
-          /* Set up indices */
-          int ijk;              /* The data point 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;
-          
-          /* these guys are easy */
-          ijk = DATINDEX(pughGH,i,j,k);               /* get linear index */
-          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 (nabla_form == 2) {
-            /* Use finite differences of g for the d's */
-            int ijkp, ijkm;
-            
-            /* X derivatives */
-            ijkp = DATINDEX(pughGH,i+1,j,k);
-            ijkm = DATINDEX(pughGH,i-1,j,k);
-            
-            tdxgxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dx);
-            tdxgxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dx);
-            tdxgxz = (gxz[ijkp] - gxz[ijkm])/(2.0*dx);
-            tdxgyy = (gyy[ijkp] - gyy[ijkm])/(2.0*dx);
-            tdxgyz = (gyz[ijkp] - gyz[ijkm])/(2.0*dx);
-            tdxgzz = (gzz[ijkp] - gzz[ijkm])/(2.0*dx);
-              
-              
-            /* Y derivatives */
-            ijkp = DATINDEX(pughGH,i,j+1,k);
-            ijkm = DATINDEX(pughGH,i,j-1,k);
-            
-            tdygxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dy);
-            tdygxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dy);
-            tdygxz = (gxz[ijkp] - gxz[ijkm])/(2.0*dy);
-            tdygyy = (gyy[ijkp] - gyy[ijkm])/(2.0*dy);
-            tdygyz = (gyz[ijkp] - gyz[ijkm])/(2.0*dy);
-            tdygzz = (gzz[ijkp] - gzz[ijkm])/(2.0*dy);
-            
-            /* X derivatives */
-            ijkp = DATINDEX(pughGH,i,j,k+1);
-            ijkm = DATINDEX(pughGH,i,j,k-1);
-            
-            tdzgxx = (gxx[ijkp] - gxx[ijkm])/(2.0*dz);
-            tdzgxy = (gxy[ijkp] - gxy[ijkm])/(2.0*dz);
-            tdzgxz = (gxz[ijkp] - gxz[ijkm])/(2.0*dz);
-            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.
-               If this looks a lot like sor_conformal.F it is... */
-            
-            /* Form upper metric - compute 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];
-            
-            /*invert metric. This is the conformal upper metric. */
-            uxx=(-SQR(gyz[ijk]) + gyy[ijk]*gzz[ijk])/det;
-            uxy=(gxz[ijk]*gyz[ijk] - gxy[ijk]*gzz[ijk])/det;
-            uyy=(-SQR(gxz[ijk]) + gxx[ijk]*gzz[ijk])/det;
-            uxz=(-gxz[ijk]*gyy[ijk] + gxy[ijk]*gyz[ijk])/det;
-            uyz=(gxy[ijk]*gxz[ijk] - gxx[ijk]*gyz[ijk])/det;
-            uzz=(-SQR(gxy[ijk]) + gxx[ijk]*gyy[ijk])/det;
-            
-            /* Coeff. Contributions from second derivative */
-            
-            /* X derivative */
-            a(-1,0,0) =         pm4 * uxx / (dx*dx);
-            a(0,0,0)  =  -two * pm4 * uxx / (dx*dx);
-            a(1,0,0)  =         pm4 * uxx / (dx*dx);
-            /* Y derivative */
-            a(0,-1,0) =         pm4 * uyy / (dy*dy);
-            a(0,0,0)  = a(0,0,0)  -  two * pm4 * uyy / (dy*dy);
-            a(0,1,0)  =         pm4 * uyy / (dy*dy);
-            /* Z derivative */
-            a(0,0,-1) =         pm4 * uzz / (dz*dz);
-            a(0,0,0)  = a(0,0,0)  -  two * pm4 * uzz / (dz*dz);
-            a(0,0,1)  =         pm4 * uzz / (dz*dz);
-            /* Mixed XY */
-            a(1,1,0)  =   two * pm4 * uxy / (dx*dy);
-            a(-1,-1,0)=   two * pm4 * uxy / (dx*dy);
-            a(1,-1,0) =  -two * pm4 * uxy / (dx*dy);
-            a(-1,1,0) =  -two * pm4 * uxy / (dx*dy);
-            /* Mixed XZ */
-            a(1,0,1)  =   two * pm4 * uxz / (dx*dz);
-            a(-1,0,-1)=   two * pm4 * uxz / (dx*dz);
-            a(1,0,-1) =  -two * pm4 * uxz / (dx*dz);
-            a(-1,0,1) =  -two * pm4 * uxz / (dx*dz);
-            /* Mixed YZ */
-            a(0,1,1)  =   two * pm4 * uyz / (dz*dy);
-            a(0,-1,-1)=   two * pm4 * uyz / (dz*dy);
-            a(0,-1,1) =  -two * pm4 * uyz / (dz*dy);
-            a(0,1,-1) =  -two * pm4 * uyz / (dz*dy);
-            
-            /*     Great so now form christoffels. Remember that
-                   
-                   G_kij = psi^4 (2 * (psi_j/psi g_ik + psi_i/psi g_jk - 
-                   psi_k/psi g_ij)
-                   + D_jik + D_ijk - D_kij)
-                   
-                   Since these are the DOWN christoffels, store them in d...
-                   
-                   NOTE however that since we will up these, we can drop the 
-                   psi^4 here and psi^-4 from the up metric. */
-            
-            
-            /*     These three have lots of cancelations */
-            dxxx = (two * psixp * gxx[ijk] + tdxgxx);
-            dxxy = (two * psiyp * gxx[ijk] + tdygxx);
-            dxxz = (two * psizp * gxx[ijk] + tdzgxx);
-            /* This one has a reduction of two identical terms */
-            dxyy = (four * psiyp * gxy[ijk] -
-                    two  * psixp * gyy[ijk] +
-                    two  * tdygxy - tdxgyy);
-            /* As does this one */
-            dxzz = (four * psizp * gxz[ijk] -
-                    two  * psixp * gzz[ijk] +
-                    two  * tdzgxz - tdxgzz);
-            /* And this one is completely general */
-            dxyz = (two * psiyp * gxz[ijk] +
-                    two * psizp * gxy[ijk] -
-                    two * psixp * gyz[ijk] +
-                    tdzgxy + tdygxz - tdxgyz);
-            
-            /* Now do it twice more without the explanations */
-            dyyy = (two * psiyp * gyy[ijk] + tdygyy);
-            dyxy = (two * psixp * gyy[ijk] + tdxgyy);
-            dyyz = (two * psizp * gyy[ijk] + tdzgyy);
-            dyxx = (four * psixp * gxy[ijk] -
-                    two  * psiyp * gxx[ijk] +
-                    two  * tdxgxy - tdygxx);
-            dyzz = (four * psizp * gyz[ijk] -
-                    two  * psiyp * gzz[ijk] +
-                    two  * tdzgyz - tdygzz);
-            dyxz = (two * psizp * gxy[ijk] +
-                    two * psixp * gyz[ijk] -
-                    two * psiyp * gxz[ijk] +
-                    tdzgxy + tdxgyz - tdygxz);
-            
-            dzzz = (two * psizp * gzz[ijk] + tdzgzz);
-            dzxz = (two * psixp * gzz[ijk] + tdxgzz);
-            dzyz = (two * psiyp * gzz[ijk] + tdygzz);
-            dzxx = (four * psixp * gxz[ijk] -
-                    two  * psizp * gxx[ijk] +
-                    two  * tdxgxz - tdzgxx);
-            dzyy = (four * psiyp * gyz[ijk] -
-                    two  * psizp * gyy[ijk] +
-                    two  * tdygyz - tdzgyy);
-            dzxy = (two * psiyp * gxz[ijk] +
-                    two * psixp * gyz[ijk] -
-                    two * psizp * gxy[ijk] +
-                    tdxgyz + tdygxz - tdzgxy);
-            
-            /* And now raise the first index */
-            Gxxx = uxx*dxxx + uxy*dyxx + uxz*dzxx;
-            Gxxy = uxx*dxxy + uxy*dyxy + uxz*dzxy;
-            Gxxz = uxx*dxxz + uxy*dyxz + uxz*dzxz;
-            Gxyy = uxx*dxyy + uxy*dyyy + uxz*dzyy;
-            Gxyz = uxx*dxyz + uxy*dyyz + uxz*dzyz;
-            Gxzz = uxx*dxzz + uxy*dyzz + uxz*dzzz;
-            
-            Gyxx = uxy*dxxx + uyy*dyxx + uyz*dzxx;
-            Gyxy = uxy*dxxy + uyy*dyxy + uyz*dzxy;
-            Gyxz = uxy*dxxz + uyy*dyxz + uyz*dzxz;
-            Gyyy = uxy*dxyy + uyy*dyyy + uyz*dzyy;
-            Gyyz = uxy*dxyz + uyy*dyyz + uyz*dzyz;
-            Gyzz = uxy*dxzz + uyy*dyzz + uyz*dzzz;
-            
-            Gzxx = uxz*dxxx + uyz*dyxx + uzz*dzxx;
-            Gzxy = uxz*dxxy + uyz*dyxy + uzz*dzxy;
-            Gzxz = uxz*dxxz + uyz*dyxz + uzz*dzxz;
-            Gzyy = uxz*dxyy + uyz*dyyy + uzz*dzyy;
-            Gzyz = uxz*dxyz + uyz*dyyz + uzz*dzyz;
-            Gzzz = uxz*dxzz + uyz*dyzz + uzz*dzzz;
-               
-
-            /*     Great. So now start adding the summed contributions
-                   from the first derivative. Note that these all have a
-                   sign change since the term comes in - ... */
-            
-            /*  g^ij G ^x_ij */
-            tmp = uxx * Gxxx + uyy * Gxyy + uzz * Gxzz +
-              two * uxy * Gxxy +
-              two * uxz * Gxxz + 
-              two * uyz * Gxyz;
-            a(1,0,0)  = a(1,0,0)  - pm4 * tmp / (two*dx);
-            a(-1,0,0) = a(-1,0,0) + pm4 * tmp / (two*dx);
-            
-            /* g^ij G^y_ij */
-            tmp = uxx * Gyxx + uyy * Gyyy + uzz * Gyzz +
-              two * uxy * Gyxy +
-              two * uxz * Gyxz + 
-              two * uyz * Gyyz;
-            a(0,1,0)  = a(0,1,0)  - pm4 * tmp / (two*dy);
-            a(0,-1,0) = a(0,-1,0) + pm4 * tmp / (two*dy);
-            
-            /* g^ij G^z_ij */
-            tmp = uxx * Gzxx + uyy * Gzyy + uzz * Gzzz +
-              two * uxy * Gzxy +
-              two * uxz * Gzxz + 
-              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);
-	      
-          } 
-
-          /* M phi */
-          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);
-          
-            /* 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.
-             
-            */
-
-          if (octant)
-            for (l=-1;l<=0;l++)
-              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 &&
-                        (ig == imin || jg == jmin || kg == kmin))  
-                      {
-                        /* We are on an inner boundary point, eg, 
-                           at an inner face */
-                        int ll,mm,nn;
-                        /* Only zero the guys at the boundaries */
-                        ll = (ig == imin ? 0 : l);
-                        mm = (jg == jmin ? 0 : m);
-                        nn = (kg == kmin ? 0 : n);
-                        a(ll,mm,nn) += a(l,m,n);
-                        a(l,m,n)     = 0.0;
-                      }
-      
-          /* renormalize */
-          if (matnormalize) {
-            ac = a(0,0,0);
-            for (J=0;J<27;J++) a[J] = a[J] / ac;
-          } 
-          else ac = 1;
-          
-
-          /* This is the new way-clever look. Note it relies
-             heavily on the index array in the workspace and
-             on multiple processors it will *NOT* make a
-             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.
-                                         */
-                  if (wsp[DATINDEX(pughGH,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)];
-                  } 
-		  else {
-                    J = wsp[DATINDEX(pughGH,i+l,j+m,k+n)];
-                    ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRQ(ierr);
-                  }
-                }
-              }
-
-          rhsval = -rhsval - Nlin[ijk] / ac;
-	  /* printf("RHS(%d %d %d)+N: %f",i,j,k,rhsval);*/
-
-          ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRQ(ierr);
-          
-        }
-      }
-    }
-  }
-  
-  if (verbose) CCTK_INFO ("Assembling the vectors");
-  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
-  ierr = VecAssemblyBegin(soln);                    CHKERRQ(ierr);
-  ierr = VecAssemblyBegin(b);                       CHKERRQ(ierr)
-  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);
-  ierr = VecAssemblyEnd(soln);                      CHKERRQ(ierr);
-  ierr = VecAssemblyEnd(b);                         CHKERRQ(ierr);
-  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr);
-
-  OptionsSetValue("-ksp_monitor","");
-  
-  CCTK_INFO("CREATING SLES");
-  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRQ(ierr);
-  
-  
-  /* 
-     Set operators. Here the matrix that defines the linear system
-     also serves as the preconditioning matrix. At a later date
-     we can probably optimize this using SAME_NONZERO_PATTERN
-     and a static trip-though flag (eg, on the second trip
-     through do something different). Also we need to think about
-     using A to precondition itself. Since I don't know what this
-     means, we'll leave it for now.  */
-
-  CCTK_INFO("CREATING SLES OPERATOR");
-  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRQ(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");
-  ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr);
-  ierr = SLESGetPC(sles,&pc);   CHKERRQ(ierr);
-  
- /* Get the PC Type */
-  if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) 
-    ierr = PCSetType(pc,PCJACOBI);
-  else if (CCTK_Equals(petsc_PC_type,"PCBJACOBI"))
-    ierr = PCSetType(pc,PCBJACOBI);
-  else if (CCTK_Equals(petsc_PC_type,"PCICC"))
-    ierr = PCSetType(pc,PCICC);
-  else if (CCTK_Equals(petsc_PC_type,"PCILU"))
-    ierr = PCSetType(pc,PCILU);
-  else if (CCTK_Equals(petsc_PC_type,"PCASM"))
-    ierr = PCSetType(pc,PCASM);
-  else if (CCTK_Equals(petsc_PC_type,"PCLU"))
-    ierr = PCSetType(pc,PCLU);
-  else if (CCTK_Equals(petsc_PC_type,"PCNONE"))
-    ierr = PCSetType(pc,PCNONE);
-  else {
-    CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n");
-    ierr = PCSetType(pc,PCNONE);
-  }
-  CHKERRQ(ierr);
-
-  
-  /* Now the same thing for the KSP Type */
-  if (CCTK_Equals(petsc_KSP_type,"KSPBCGS")) 
-    ierr = KSPSetType(ksp,KSPBCGS); 
-  else if (CCTK_Equals(petsc_KSP_type,"KSPRICHARDSON"))
-    ierr = KSPSetType(ksp,KSPRICHARDSON);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPCHEBYCHEV"))
-    ierr = KSPSetType(ksp,KSPCHEBYCHEV);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPCG"))
-    ierr = KSPSetType(ksp,KSPCG);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPGMRES"))
-    ierr = KSPSetType(ksp,KSPGMRES);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPCGS"))
-    ierr = KSPSetType(ksp,KSPCGS);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPTFQMR"))
-    ierr = KSPSetType(ksp,KSPTFQMR);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPTCQMR"))
-    ierr = KSPSetType(ksp,KSPTCQMR);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPCR"))
-    ierr = KSPSetType(ksp,KSPCR);
-  else if (CCTK_Equals(petsc_KSP_type,"KSPLSQR"))
-    ierr = KSPSetType(ksp,KSPLSQR);
-  else {
-    CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n");
-    ierr = KSPSetType(ksp,KSPBCGS);
-  }
-  CHKERRQ(ierr);
-
-
-
- /* Set up tolerances */
-  
-  if (PetscTolStyle == ELLCONV_ABSOLUTE) {
-    rtol = 1.0e-15; 
-    atol = tolerance;
-  } else if (PetscTolStyle == ELLCONV_RELATIVE) {
-    rtol = tolerance; 
-    atol = 1.0e-15;
-  } else if (PetscTolStyle == ELLCONV_EITHER) {
-    rtol = tolerance; 
-    atol = tolerance;
-  } else {
-    printf("PETSC Solver: PetscTolStyle set incorrectly [%d]\n",
-            PetscTolStyle);
-  }
-  
-  ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT,
-                          PETSC_DEFAULT);
-  CHKERRQ(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); 
-  CHKERRQ(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); CHKERRQ(ierr);
-  
-  /* 
-    Set runtime options. Since we don't use PETSC runtime options
-    but rather use the CACTUS parameter parser, we do this before
-    we parse things out. But that may not be such a good idea.
-    So lets put it here.
-  */
-  ierr = SLESSetFromOptions(sles);   CHKERRQ(ierr);
-    
-
-  /* OK so finally we are able to ... */
-
-  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
-                      Solve the linear system
-     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-  CCTK_INFO("SLES solve...");
-  ierr = SLESSolve(sles,b,soln,&its); CHKERRQ(ierr); 
-  
-  /* Here we can form a "res = Ax  - b" and find its norm to get
-     an idea of how well we've converged. Put this in later
-     but we can do it using MatMultAdd() after we flip the sign
-     on b with VecScale() (into res, a new vector, then find the
-     norm of res with VecNorm()
-  */
-  
-  /* Now since we now have local layout, we can just get the values
-     from the soln vector now matter how many processors we are
-     on (eg, the VecCreateMPI has made the solution local) */
-
-  VecGetArray(soln,&values); 
-  for (k=kmin;k<kmax;k++) {
-    for (j=jmin;j<jmax;j++) {
-      for (i=imin;i<imax;i++) {
-        /* Set up indices */
-        int ijk;                /* The data point 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;
-
-        /* these guys are easy */
-        ijk = DATINDEX(pughGH,i,j,k);
-        if (wsp[ijk] >= 0) {
-          ig = i + GH->cctk_lbnd[0];
-          jg = j + GH->cctk_lbnd[1];
-          kg = k + GH->cctk_lbnd[2];
-          
-          /* 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;
-          ell_field[ijk] = values[I];
-        }
-      }
-    }
-  }
-  
-  VecRestoreArray(soln,&values);   
-
-  /* Sync var in CCTK speak. */ 
-  PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
-  
-  /* And finally free up the matrix memory   FIXME 
-  ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr); */
-
-  /* This code is not used anymore */
-  if (verbose) CCTK_INFO("Destroying Matrices");
-  ierr = SLESDestroy(sles); CHKERRQ(ierr);
-  ierr = VecDestroy(soln);  CHKERRQ(ierr);
-  ierr = VecDestroy(b);     CHKERRQ(ierr);  
-  ierr = MatDestroy(A[0]);  CHKERRQ(ierr);
-  free(A);
-    
-  /*$PetscFinalize();$*/
-
-  if (verbose) CCTK_INFO("LEAVING ELLPETSC");
-  
-}