aboutsummaryrefslogtreecommitdiff
path: root/src/petsc_confmetric_solver.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/petsc_confmetric_solver.c')
-rw-r--r--src/petsc_confmetric_solver.c254
1 files changed, 127 insertions, 127 deletions
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;i<pEx->npoints;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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 05:43:50 +0000