From be2349567ed422ab31f244f7113628cada2d2c6c Mon Sep 17 00:00:00 2001 From: goodale Date: Mon, 17 May 2004 12:28:57 +0000 Subject: Untabified. Please try to adhere to the Cactus coding guidelines when modifying or adding files. git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@87 1d96b42b-98df-4a6a-9d84-1b24288d4588 --- src/Startup.c | 6 +- src/petsc_confmetric_solver.c | 254 +++++++++++++++++++++--------------------- src/petsc_flat_solver.c | 116 ++++++++++--------- src/petsc_wrapper.c | 26 ++--- 4 files changed, 200 insertions(+), 202 deletions(-) diff --git a/src/Startup.c b/src/Startup.c index 05212b1..bdadc2c 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -12,13 +12,13 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllPETSc_Startup_c) void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI, - int *NI, int *AbsTol, int *RelTol); + int *NI, int *AbsTol, int *RelTol); void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI, - int *NI, int *AbsTol, int *RelTol); + int *NI, int *AbsTol, int *RelTol); void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, - int *AbsTol, int *RelTol); + int *AbsTol, int *RelTol); /* Registration of the petsc solvers with the Elliptic solver registry. This routine registers petsc_confmetric under the name "petsc" for diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c index c4bfb5c..7600fd8 100644 --- a/src/petsc_confmetric_solver.c +++ b/src/petsc_confmetric_solver.c @@ -80,8 +80,8 @@ 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); + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); void *GetDataPtr_NextTL(cGH *GH, const char *field); @@ -105,8 +105,8 @@ void *GetDataPtr_NextTL(cGH *GH, const char *field) { 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) { + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { DECLARE_CCTK_PARAMETERS /* CCTK passed parameters */ @@ -284,19 +284,19 @@ int 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=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[0]); + GH->cctk_nghostzones[0]); imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : - pEx->lnsize[0]-GH->cctk_nghostzones[0]); + pEx->lnsize[0]-GH->cctk_nghostzones[0]); jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[1]); + GH->cctk_nghostzones[1]); jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : - pEx->lnsize[1]-GH->cctk_nghostzones[1]); + pEx->lnsize[1]-GH->cctk_nghostzones[1]); kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[2]); + GH->cctk_nghostzones[2]); kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : - pEx->lnsize[2]-GH->cctk_nghostzones[2]); + 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 @@ -388,12 +388,12 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 */ + (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); @@ -421,7 +421,7 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, else if (CCTK_EQUALS(petsc_nablaform,"up")) { if (verbose) - CCTK_INFO("Forming nabla with upper g and finite difference of dg \n"); + 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)); @@ -430,42 +430,42 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 det; - 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; - } + CCTK_REAL det; + 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; + /*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); + /* 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"); @@ -498,22 +498,22 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 (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; @@ -551,9 +551,9 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 */ + + /* 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]) - @@ -707,64 +707,64 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 */ + 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]; + 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); + VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES); - /* Put in the octant boundary conditions. + /* 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 @@ -789,12 +789,12 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ #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"); - } + 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) @@ -839,19 +839,19 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, /* This is a boundary. */ rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)]; } - else { + 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); + CHKERRQ(ierr); } } } - if (Nstorage) - rhsval = -rhsval - Nlin[ijk] / ac; + if (Nstorage) + rhsval = -rhsval - Nlin[ijk] / ac; ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); - CHKERRQ(ierr); + CHKERRQ(ierr); } } @@ -869,7 +869,7 @@ int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, if (nabla_form == 3) { if (verbose) - printf ("Freeing upper metric storage\n"); + printf ("Freeing upper metric storage\n"); free(uxx3); free(uxy3); free(uxz3); diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c index dc09aa3..645502a 100644 --- a/src/petsc_flat_solver.c +++ b/src/petsc_flat_solver.c @@ -64,10 +64,10 @@ static Vec soln, b; /* approx solution, RHS */ static SLES sles; /* linear solver context */ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, - int *AbsTol, int *RelTol); + int *AbsTol, int *RelTol); int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, - int *AbsTol, int *RelTol) + int *AbsTol, int *RelTol) { DECLARE_CCTK_PARAMETERS /* CCTK passed parameters */ @@ -113,7 +113,7 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, int octant; /* Apply octant BCs inside */ int matnormalize; /* Normalize the central mat value to one? */ int I,J; /* The PETSc col and row in the matrix */ - + int PetscTolStyle; @@ -206,19 +206,19 @@ int 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=((pCon->neighbours[pughGH->myproc][XDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[0]); + GH->cctk_nghostzones[0]); imax=((pCon->neighbours[pughGH->myproc][XDP]<0) ? pEx->lnsize[0]-1 : - pEx->lnsize[0]-GH->cctk_nghostzones[0]); + pEx->lnsize[0]-GH->cctk_nghostzones[0]); jmin=((pCon->neighbours[pughGH->myproc][YDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[1]); + GH->cctk_nghostzones[1]); jmax=((pCon->neighbours[pughGH->myproc][YDP]<0) ? pEx->lnsize[1]-1 : - pEx->lnsize[1]-GH->cctk_nghostzones[1]); + pEx->lnsize[1]-GH->cctk_nghostzones[1]); kmin=((pCon->neighbours[pughGH->myproc][ZDM]<0) && !(octant) ? 1 : - GH->cctk_nghostzones[2]); + GH->cctk_nghostzones[2]); kmax=((pCon->neighbours[pughGH->myproc][ZDP]<0) ? pEx->lnsize[2]-1 : - pEx->lnsize[2]-GH->cctk_nghostzones[2]); + 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 @@ -317,12 +317,12 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 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 */ + (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); @@ -374,34 +374,34 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 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]; /* M phi */ - if (Mstorage) - { - a(0,0,0) = a(0,0,0) + Mlin[ijk]; - } + if (Mstorage) + { + a(0,0,0) = a(0,0,0) + Mlin[ijk]; + } /* Now set the values of the matrix. Here we force dirichlet. */ - VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES); + VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES); #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"); - } + 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) @@ -425,12 +425,12 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, /* renormalize */ if (matnormalize) - { + { ac = a(0,0,0); for (J=0;J<27;J++) a[J] = a[J] / ac; } else - ac = 1; + ac = 1; /* This is the new way-clever look. Note it relies @@ -439,40 +439,40 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 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. - */ + { /* 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 - { + 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); + CHKERRQ(ierr); } } } - } - } - if (Nstorage) - { - rhsval = -rhsval - Nlin[ijk] / ac; - } + } + } + if (Nstorage) + { + rhsval = -rhsval - Nlin[ijk] / ac; + } ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); - CHKERRQ(ierr); + CHKERRQ(ierr); - /* reset the central value*/ - a( 0, 0, 0) = -6.0; + /* reset the central value*/ + a( 0, 0, 0) = -6.0; } /* if (wsp>=0) */ } /*for i */ @@ -667,7 +667,7 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, int ijk; /* The data point for the array */ int I; /* The col in the matrix */ - /* these guys are easy */ + /* these guys are easy */ ijk = DATINDEX(pEx,i,j,k); if (wsp[ijk] >= 0) { @@ -707,5 +707,3 @@ int petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, return (0); } - - diff --git a/src/petsc_wrapper.c b/src/petsc_wrapper.c index 8d4fbcf..0f6fb65 100644 --- a/src/petsc_wrapper.c +++ b/src/petsc_wrapper.c @@ -8,8 +8,8 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_wrapper_c) void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex, - int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); + int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); /* The wrapper functions for the core PETSc solver, that performs the actual solve. One wrapper function for each class, because these @@ -18,14 +18,14 @@ void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex, /* Wrapper function for the class of elliptic equations that needs a metric */ void petsc_metric(cGH *GH, int *MetricI, int FieldIndex, - int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { + int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { /* Arrays for the M/N size info */ int petsc_confmetric_solver(cGH *GH, int *MetricI, int MetricISize, - int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); if (GH->cctk_dim>3) @@ -34,22 +34,22 @@ void petsc_metric(cGH *GH, int *MetricI, int FieldIndex, } petsc_confmetric_solver(GH, MetricI, 6, FieldIndex, - MIndex, NIndex, - AbsTol, RelTol); + MIndex, NIndex, + AbsTol, RelTol); } /* wrapper function for the class of elliptic equations, that needs a conf. factor */ void petsc_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex, - int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { + int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { int petsc_confmetric_solver(cGH *, int *, int, int, int, int, - CCTK_REAL *, CCTK_REAL *); + CCTK_REAL *, CCTK_REAL *); /* petsc_confmetric_solver expects an metric array, in the case of this equation class, needs conf.factor as last entry -> size 7 */ petsc_confmetric_solver(GH, MetricPsiI, 7, - FieldIndex, MIndex, NIndex, - AbsTol, RelTol); + FieldIndex, MIndex, NIndex, + AbsTol, RelTol); } -- cgit v1.2.3