diff options
Diffstat (limited to 'src/petsc_flat_solver.c')
| -rw-r--r-- | src/petsc_flat_solver.c | 116 |
1 files changed, 57 insertions, 59 deletions
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); } - - |