/*@@ @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 #include "cctk.h" #include "cctk_Parameters.h" #include #include /* Include some definitions for backwards compatibility */ #include "petsc_defines.h" #include "ellpetsc.h" #include "CactusPUGH/PUGH/src/include/pugh.h" static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_confmetric_solver_c) /*#define DEBUG*/ /*Don't know what these actually mean! FIXME */ #define ELLCONV_ABSOLUTE 0 #define ELLCONV_RELATIVE 1 #define ELLCONV_EITHER 2 /* direction macros */ #define XDP 1 #define XDM 0 #define YDP 3 #define YDM 2 #define ZDP 5 #define ZDM 4 /* 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(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 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 */ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); void *GetDataPtr_NextTL(cGH *GH, const char *field); void *GetDataPtr_NextTL(cGH *GH, const char *field) { int index; 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 */ PC pc; /* preconditioner context */ KSP ksp; /* Krylov subspace method context */ int num_A; /* Number of A-arrays as needed ny Dan's MG */ int ierr; /* Check the return status of petsc */ int retcode; /* Check the return status of CCTK */ double a[27]; /* 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; pGH *pughGH; /* The pugh Extension handle */ pGExtras *pEx; pConnectivity *pCon; int myproc; /* out processor */ /* Tolerances */ CCTK_REAL rtol=0, atolerance=0, tolerance; /* Values to assemble the matrix */ CCTK_REAL two=2.0, four=4.0, tmp, det; CCTK_REAL pm4, psixp,psiyp,psizp; CCTK_REAL dx,dy,dz; CCTK_REAL uxx=0,uxy=0,uxz=0,uyy=0,uyz=0,uzz=0; 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 conformal=0; /* Do we have conformal metric ? */ int nabla_form=0; /* Which form of the nable */ int matnormalize; /* Normalize the central mat value to one? */ CCTK_REAL ac; /* Storage for a(0,0,0) for renorm */ int PetscTolStyle; /* For the upper metric form of nabla */ CCTK_REAL *uxx3=NULL, *uxy3=NULL, *uxz3=NULL, *uyy3=NULL, *uyz3=NULL, *uzz3=NULL; /* 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; int Mstorage=0, Nstorage=0; RelTol = RelTol; octant = CCTK_Equals(domain,"octant"); verbose = CCTK_Equals(petsc_verbose,"yes")|| CCTK_Equals(petsc_verbose,"debug"); debug = CCTK_Equals(petsc_verbose,"debug"); matnormalize = 0; if (MetricPsiISize==7) conformal=1; else if (MetricPsiISize==7) conformal=0; else CCTK_WARN(0,"Size of the Metric must be either 7 (metric+conformal) or 6 (metric)"); /* FIXME, the TolAbs/TolRel will be evaluated here */ PetscTolStyle=0; tolerance =AbsTol[0]; /* Get the link to pugh Extension */ pughGH = PUGH_pGH(GH); 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; char **argv; if (debug) { #ifdef 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); CHKERRQ(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); 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"); } /* if we have a negative index, this GF is not needed, there for don't even look for it. when index positive, set flag Mstorage=1, dito for N */ if (MIndex>=0) { Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex); Mstorage = 1; } if (NIndex>=0) { Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex); Nstorage = 1; } /* Get the workspace data pointer */ 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;inpoints;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=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : GH->cctk_nghostzones[0]); imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : pEx->lnsize[0]-GH->cctk_nghostzones[0]); jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : GH->cctk_nghostzones[1]); jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : pEx->lnsize[1]-GH->cctk_nghostzones[1]); kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : 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 (irnsize[i][0]; nxs=((pCon->neighbours[i][XDM]<0) && !(octant) ? nxs-1 : nxs-GH->cctk_nghostzones[0]); nxs=((pCon->neighbours[i][XDP]<0) ? nxs-1 : nxs-GH->cctk_nghostzones[0]); nys=pEx->rnsize[i][1]; nys=((pCon->neighbours[i][YDM]<0) && !(octant) ? nys-1 : nys-GH->cctk_nghostzones[1]); nys=((pCon->neighbours[i][YDP]<0) ? nys-1 : nys-GH->cctk_nghostzones[1]); nzs=pEx->rnsize[i][2]; nzs=((pCon->neighbours[i][ZDM]<0) && !(octant) ? 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; } #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=kmin;kdata[ijk]) as well as the column for the stencil of the neighbors (workspace->data[DATINDEX(GH,i+1,j,k)] etc...). So onwards */ 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",pEx->lnsize[0],pEx->lnsize[1],pEx->lnsize[2]); } #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 (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); /* Compare the PETSc layout to Cactus, this better be a match */ ierr = VecGetOwnershipRange(soln,&pvstart,&pvend); ierr = MatGetOwnershipRange(A[0],&pstart,&pend); #ifdef DEBUG printf("CAC M-Layout: %d %d \n",startpoint, endpoint); printf("PET M-Layout: %d %d \n",pstart, pend); printf("PET V-Layout: %d %d \n",pvstart,pvend); #endif if (pstart != startpoint && pend != endpoint) CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why?)\n"); /* Decide on the nabla form in PETSc: */ if (CCTK_EQUALS(petsc_nablaform,"down")) { if (verbose) CCTK_INFO("Forming nabla with lower g and finite difference of dg \n"); nabla_form = 2; } else if (CCTK_EQUALS(petsc_nablaform,"up")) { if (verbose) 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)); 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;inpoints;i++) { 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; /* 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"); if (verbose) CCTK_INFO("Creating the coefficient matrix..."); if (verbose && matnormalize) CCTK_INFO ("...with diagonal renormalized to one"); for (k=kmin;k= 0) { 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(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]; /* 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(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); 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(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); 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(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); 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. 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); } /* 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 */ /* M phi */ 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); /* Put in the octant boundary conditions. We do these by reflecting the -1 stencil if we are at a boundary at -1. That is, imagine the 1D laplace equation with the boundary condition a(0) = a(1). This means the point 1 stencil (which is our first stencil in the matrix) becomes a(0) + a(2) - 2 a(1) = rho(i) * deltax -> a(2) - 2 a(1) + a(1) = rho(i) * deltax OK so think about this with a general stenci S_0 a_0 + S_1 a_1 + S_2 a_2 = rho_1 -> S_0 a_1 + S_1 a_1 + S_2 a_2 = rho_1 -> (S_0 + S_1) a_1 + S_2 a_2 = rho_1 eg, S_0 -> 0 and S_1 -> S_1 + S_0 Great, so implement this in 3D. It is a bit trickier, since we don't want to zero neighbors in the wrong direction, but not impossible. */ #ifdef DEBUG if (ijk==1918) { for (l=-1;l<=1;l++) for (m=-1;m<=1;m++) for (n=-1;n<=1;n++) printf(" (%d %d %d): %f \n",l,m,n,a(l,m,n)); printf("\n"); } #endif if (octant) for (l=-1;l<=0;l++) for (m=-1;m<=0;m++) for (n=-1;n<=0;n++) if (l*m*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, 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(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 { 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); } } } if (Nstorage) rhsval = -rhsval - Nlin[ijk] / ac; 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); #if PETSC_VERSION_MAJOR < 3 ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr); #else ierr = MatSetOption(A[0],MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); CHKERRQ(ierr); #endif if (nabla_form == 3) { if (verbose) printf ("Freeing upper metric storage\n"); free(uxx3); free(uxy3); free(uxz3); free(uyy3); free(uyz3); free(uzz3); } if (trips==0) { #if PETSC_VERSION_MAJOR < 2 || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 0) OptionsSetValue("-ksp_monitor",""); #else PetscOptionsSetValue("-ksp_monitor",""); #endif } if (verbose) { 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. */ if (verbose) 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. */ if (verbose) 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; atolerance = tolerance; } else if (PetscTolStyle == ELLCONV_RELATIVE) { rtol = tolerance; atolerance = 1.0e-15; } else if (PetscTolStyle == ELLCONV_EITHER) { rtol = tolerance; atolerance = tolerance; } else { printf("PETSC Solver: PetscTolStyle set incorrectly [%d]\n", PetscTolStyle); } ierr = KSPSetTolerances(ksp,rtol,atolerance,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); */ if (verbose) CCTK_INFO("KSPSetInitialGuess\n"); #if PETSC_VERSION_MAJOR < 2 || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR <= 0 && PETSC_VERSION_SUBMINOR < 5) ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr); #else ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr); #endif /* 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (verbose) 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= 0) { /* 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. */ CCTK_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"); return (0); }