diff options
-rw-r--r-- | src/petsc_confmetric.c | 3 | ||||
-rw-r--r-- | src/petsc_confmetric_solver.c | 115 | ||||
-rw-r--r-- | src/petsc_flat_solver.c | 70 |
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]; |