diff options
Diffstat (limited to 'src/petsc_confmetric_solver.c')
| -rw-r--r-- | src/petsc_confmetric_solver.c | 55 |
1 files changed, 33 insertions, 22 deletions
diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c index c60312e..9a2e9d7 100644 --- a/src/petsc_confmetric_solver.c +++ b/src/petsc_confmetric_solver.c @@ -150,11 +150,12 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_REAL *psi =NULL; CCTK_REAL *psix=NULL, *psiy=NULL, *psiz=NULL; - + int Mstorage=0, Nstorage=0; octant = CCTK_Equals(domain,"octant"); - verbose = CCTK_Equals(petsc_verbose,"yes")||CCTK_Equals(petsc_verbose,"debug"); + verbose = CCTK_Equals(petsc_verbose,"yes")|| + CCTK_Equals(petsc_verbose,"debug"); debug = CCTK_Equals(petsc_verbose,"debug"); reuse = 0; matnormalize = 0; @@ -213,9 +214,6 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, /* 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]); @@ -230,9 +228,20 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 */ @@ -261,11 +270,11 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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.) */ + /* 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++) { @@ -402,12 +411,12 @@ 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 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uxy3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uxz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uyy3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uyz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uzz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); + 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++) { CCTK_REAL det; CCTK_REAL p12; @@ -436,8 +445,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, det = det*p12; /* Rescaling for the uppermetric solver */ - Mlin[i] = Mlin[i]*sqrt(det); - Nlin[i] = Nlin[i]*sqrt(det); + 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); @@ -736,7 +745,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } /* end nabla_form=3 */ /* M phi */ - a(0,0,0) = a(0,0,0) - Mlin[ijk]; + 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 @@ -818,8 +828,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } } - rhsval = -rhsval + Nlin[ijk] / ac; - /* printf("RHS(%d %d %d)+N: %f",i,j,k,rhsval);*/ + if (Nstorage) + rhsval = -rhsval + Nlin[ijk] / ac; ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); CHKERRA(ierr); @@ -1062,3 +1072,4 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, if (verbose) CCTK_INFO("LEAVING ELLPETSC"); } + |