aboutsummaryrefslogtreecommitdiff
path: root/src/petsc_flat_solver.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/petsc_flat_solver.c')
-rw-r--r--src/petsc_flat_solver.c251
1 files changed, 164 insertions, 87 deletions
diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c
index 74e0216..0a8ec9e 100644
--- a/src/petsc_flat_solver.c
+++ b/src/petsc_flat_solver.c
@@ -15,8 +15,8 @@
#include "cctk_Parameters.h"
#include "mpi.h"
-#include "sles.h"
-
+#include "petscsles.h"
+#include "ellpetsc.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"
static const char *rcsid = "$Header$";
@@ -48,7 +48,8 @@ static SLES sles; /* linear solver context */
void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
- int *AbsTol, int *RelTol) {
+ int *AbsTol, int *RelTol)
+{
DECLARE_CCTK_PARAMETERS /* CCTK passed parameters */
@@ -105,10 +106,6 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
/* flags to signal if storage is ON or OFF (used/not-used) for M/N*/
int Mstorage=0, Nstorage=0;
-
-
-
-
octant = CCTK_Equals(domain,"octant");
verbose = CCTK_Equals(petsc_verbose,"yes")||
CCTK_Equals(petsc_verbose,"debug");
@@ -122,14 +119,18 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
/* 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");
+ if (!pughGH)
+ {
+ CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle");
+ }
/* Get the extras extension for 3D grid functions */
pEx = pughGH->GFExtras[2];
pCon = pughGH->Connectivity[2];
/* Things to do on first iteration */
- if (trips==0) {
+ if (trips==0)
+ {
int argc;
char **argv;
@@ -139,7 +140,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
#endif
argc = CCTK_CommandLine(&argv);
- ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRA(ierr);
+ ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRQ(ierr);
PetscInitialize(&argc,&argv,NULL,NULL);
CCTK_INFO("PETSc initialized");
@@ -157,11 +158,13 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
/* 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) {
+ if (MIndex>=0)
+ {
Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
Mstorage = 1;
}
- if (NIndex>=0) {
+ if (NIndex>=0)
+ {
Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex);
Nstorage = 1;
}
@@ -171,7 +174,10 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
/* initialize the linear index lookup table with -1;
this will indicate a boudary later */
- for (i=0;i<pEx->npoints;i++) wsp[i] = -1.0;
+ for (i=0;i<pEx->npoints;i++)
+ {
+ wsp[i] = -1.0;
+ }
/* Get myproc from the CCTK, get the gridspacing */
myproc = PUGH_MyProc(GH);
@@ -199,7 +205,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
(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++) {
+ for (i=0;i<myproc;i++)
+ {
nxs=pEx->rnsize[i][0];
nxs=((pCon->neighbours[i][XDM]<0) && !(octant) ? nxs-1 :
@@ -241,7 +248,10 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
do a sync that means each processor knows the indiced
in the ghost zones (eg, on its neighbors) in a FD stencil of 1 */
retcode = PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
- if (retcode<0) CCTK_WARN(1,"Synchronization failed\n");
+ if (retcode<0)
+ {
+ CCTK_WARN(1,"Synchronization failed");
+ }
/* So woohoo. Now for each point in our ijk space, we have
information about our row in the matrix (workspace->data[ijk])
@@ -253,7 +263,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
nzs = pEx->nsize[2]-2;
#ifdef DEBUG
- if (debug) {
+ if (debug)
+ {
printf("nxyzs %d %d %d \n",nxs,nys,nzs);
printf("lnxyz: %d %d %d \n",pEx->lnsize[0],pEx->lnsize[1],pEx->lnsize[2]);
}
@@ -291,23 +302,24 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
19,PETSC_NULL, /* Diagonal */
19,PETSC_NULL, /* Off diagonal */
&A[0]); /* The output */
- CHKERRA(ierr);
+ CHKERRQ(ierr);
ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
- CHKERRA(ierr);
+ CHKERRQ(ierr);
ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
- CHKERRA(ierr);
+ CHKERRQ(ierr);
/* Compare the PETSc layout to Cactus, this better be a match */
ierr = VecGetOwnershipRange(soln,&pvstart,&pvend);
ierr = MatGetOwnershipRange(A[0],&pstart,&pend);
- if (debug) {
+ if (debug)
+ {
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");
+ CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why??)");
/* Initialize the stencil array */
@@ -321,11 +333,15 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
a( 0, 0,-1) = 1.0; /* ab */
- for (k=kmin;k<kmax;k++) {
- for (j=jmin;j<jmax;j++) {
- for (i=imin;i<imax;i++) {
+ for (k=kmin;k<kmax;k++)
+ {
+ for (j=jmin;j<jmax;j++)
+ {
+ for (i=imin;i<imax;i++)
+ {
- if (wsp[DATINDEX(pEx,i,j,k)] >= 0) {
+ if (wsp[DATINDEX(pEx,i,j,k)] >= 0)
+ {
/* Set up indices */
int ijk; /* The lineax index for the array */
@@ -342,16 +358,27 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
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];
+ }
- /* Great now set the values of the matrix. Here we force dirichlet. */
+ /* Now set the values of the matrix. Here we force dirichlet. */
VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
#ifdef DEBUG
- if (ijk==1918) {
+ 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));
+ {
+ for (n=-1;n<=1;n++)
+ {
+ printf(" (%d %d %d): %f \n",l,m,n,a(l,m,n));
+ }
+ }
+ }
printf("\n");
}
#endif
@@ -376,7 +403,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
}
/* renormalize */
- if (matnormalize) {
+ if (matnormalize)
+ {
ac = a(0,0,0);
for (J=0;J<27;J++) a[J] = a[J] / ac;
}
@@ -390,27 +418,37 @@ void 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.
- */
- if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) {
+ {
+ 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(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);
- CHKERRA(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);
- CHKERRA(ierr);
+ CHKERRQ(ierr);
/* reset the central value*/
a( 0, 0, 0) = -6.0;
@@ -420,20 +458,29 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
} /* for j */
} /* for k */
- if (verbose) CCTK_INFO ("Assembling the vectors");
- ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
- ierr = VecAssemblyBegin(soln); CHKERRA(ierr);
- ierr = VecAssemblyBegin(b); CHKERRA(ierr)
- ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
- ierr = VecAssemblyEnd(soln); CHKERRA(ierr);
- ierr = VecAssemblyEnd(b); CHKERRA(ierr);
- ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(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);
if (trips==0)
- OptionsSetValue("-ksp_monitor","");
-
- if (verbose) CCTK_INFO("CREATING SLES");
- ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles); CHKERRA(ierr);
+ {
+ PetscOptionsSetValue("-ksp_monitor","");
+ }
+
+ if (verbose)
+ {
+ CCTK_INFO("CREATING SLES");
+ }
+
+ ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles); CHKERRQ(ierr);
/*
@@ -445,16 +492,22 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
using A to precondition itself. Since I don't know what this
means, we'll leave it for now. */
- if (verbose) CCTK_INFO("CREATING SLES OPERATOR");
- ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr);
+ if (verbose)
+ {
+ 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. */
- if (verbose) CCTK_INFO("SLESGet KSP/PC");
- ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr);
- ierr = SLESGetPC(sles,&pc); CHKERRA(ierr);
+ if (verbose)
+ {
+ 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"))
@@ -471,11 +524,12 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
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");
+ else
+ {
+ CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE");
ierr = PCSetType(pc,PCNONE);
}
- CHKERRA(ierr);
+ CHKERRQ(ierr);
/* Now the same thing for the KSP Type */
@@ -499,36 +553,47 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
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");
+ else
+ {
+ CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB");
ierr = KSPSetType(ksp,KSPBCGS);
}
- CHKERRA(ierr);
+ CHKERRQ(ierr);
/* Set up tolerances */
- if (PetscTolStyle == ELLCONV_ABSOLUTE) {
+ if (PetscTolStyle == ELLCONV_ABSOLUTE)
+ {
rtol = 1.0e-15;
atol = tolerance;
- } else if (PetscTolStyle == ELLCONV_RELATIVE) {
+ }
+ else if (PetscTolStyle == ELLCONV_RELATIVE)
+ {
rtol = tolerance;
atol = 1.0e-15;
- } else if (PetscTolStyle == ELLCONV_EITHER) {
+ }
+ else if (PetscTolStyle == ELLCONV_EITHER)
+ {
rtol = tolerance;
atol = tolerance;
- } else {
+ }
+ else
+ {
printf("PETSC Solver: PetscTolStyle set incorrectly [%d]\n",
PetscTolStyle);
}
ierr = KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, PETSC_DEFAULT);
- CHKERRA(ierr);
+ CHKERRQ(ierr);
- if (verbose) CCTK_INFO("KSPSetInitialGuess\n");
+ if (verbose)
+ {
+ CCTK_INFO("KSPSetInitialGuess");
+ }
ierr = KSPSetInitialGuessNonzero(ksp);
- CHKERRA(ierr);
+ CHKERRQ(ierr);
/*
Set runtime options. Since we don't use PETSC runtime options
@@ -536,7 +601,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
we parse things out. But that may not be such a good idea.
So lets put it here.
*/
- ierr = SLESSetFromOptions(sles); CHKERRA(ierr);
+ ierr = SLESSetFromOptions(sles); CHKERRQ(ierr);
/* OK so finally we are able to ... */
@@ -544,8 +609,12 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
- if (verbose) CCTK_INFO("SLES solve...");
- ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr);
+ if (verbose)
+ {
+ 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
@@ -559,20 +628,22 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
on (eg, the VecCreateMPI has made the solution local) */
VecGetArray(soln,&values);
- for (k=kmin;k<kmax;k++) {
- for (j=jmin;j<jmax;j++) {
- for (i=imin;i<imax;i++) {
+ for (k=kmin;k<kmax;k++)
+ {
+ for (j=jmin;j<jmax;j++)
+ {
+ for (i=imin;i<imax;i++)
+ {
/* Set up indices */
int ijk; /* The data point for the array */
int I; /* The col in the matrix */
/* these guys are easy */
ijk = DATINDEX(pEx,i,j,k);
- if (wsp[ijk] >= 0) {
-
+ if (wsp[ijk] >= 0)
+ {
/* Now this one. "Fortran-order" the matrix. But remember
- we have stripped off the ghost zones. Hence ig-1 and nxs...
- */
+ we have stripped off the ghost zones. Hence ig-1 and nxs */
I = wsp[ijk]-startpoint;
ell_field[ijk] = values[I];
}
@@ -585,19 +656,25 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
/* And finally free up the matrix memory */
- /*$ierr = MatReleaseValuesMemory(A); CHKERRA(ierr);$*/
+ /*$ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr);$*/
/* This code is not used anymore */
- if (verbose) CCTK_INFO("Destroying Matrices");
- ierr = SLESDestroy(sles); CHKERRA(ierr);
- ierr = VecDestroy(soln); CHKERRA(ierr);
- ierr = VecDestroy(b); CHKERRA(ierr);
- ierr = MatDestroy(A[0]); CHKERRA(ierr);
+ 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");
+ if (verbose)
+ {
+ CCTK_INFO("LEAVING ELLPETSC");
+ }
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 03:48:11 +0000