From ff1ddc9e1ca1aab1b700bfa8469cf57e333d9ef3 Mon Sep 17 00:00:00 2001 From: schnetter Date: Wed, 5 May 2004 15:00:00 +0000 Subject: Remove unused file. git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@86 1d96b42b-98df-4a6a-9d84-1b24288d4588 --- src/petsc_confmetric.c | 942 ------------------------------------------------- 1 file changed, 942 deletions(-) delete mode 100644 src/petsc_confmetric.c 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;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=((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 (irnx[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;kdata[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= 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= 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"); - -} -- cgit v1.2.3