diff options
Diffstat (limited to 'src/petsc_confmetric_solver.c')
| -rw-r--r-- | src/petsc_confmetric_solver.c | 52 |
1 files changed, 34 insertions, 18 deletions
diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c index 59441d3..7ef205e 100644 --- a/src/petsc_confmetric_solver.c +++ b/src/petsc_confmetric_solver.c @@ -10,12 +10,18 @@ @@*/ +#include <stdlib.h> + #include "cctk.h" #include "cctk_Parameters.h" -#include "assert.h" -#include "mpi.h" +#include "petscversion.h" +#if PETSC_VERSION_MAJOR < 2 || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 0) +#include "sles.h" +#else #include "petscsles.h" +#endif + #include "ellpetsc.h" #include "CactusPUGH/PUGH/src/include/pugh.h" @@ -61,6 +67,11 @@ static Vec soln, b; /* approx solution, RHS */ 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); +void *GetDataPtr_NextTL(cGH *GH, const char *field); + void *GetDataPtr_NextTL(cGH *GH, const char *field) { int index; @@ -81,15 +92,13 @@ void *GetDataPtr_NextTL(cGH *GH, const char *field) { /* This passing matches that of in the convention in the elliptic registration routine see LinearElliptic.h*/ -void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, +int petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { DECLARE_CCTK_PARAMETERS /* CCTK passed parameters */ - void *GetDataPtr_NextTL(cGH *GH, const char *field); - PC pc; /* preconditioner context */ KSP ksp; /* Krylov subspace method context */ int num_A; /* Number of A-arrays as needed ny Dan's MG */ @@ -113,13 +122,13 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, int myproc; /* out processor */ /* Tolerances */ - CCTK_REAL rtol, atol, tolerance; + CCTK_REAL rtol=0, atolerance=0, tolerance; /* Values to assemble the matrix */ - CCTK_REAL two=2.0, four=4.0, tmp; - CCTK_REAL pm4, psixp,psiyp,psizp,det; + CCTK_REAL two=2.0, four=4.0, tmp, det; + CCTK_REAL pm4, psixp,psiyp,psizp; CCTK_REAL dx,dy,dz; - CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz; + CCTK_REAL uxx=0,uxy=0,uxz=0,uyy=0,uyz=0,uzz=0; CCTK_REAL Gxxx,Gxxy,Gxxz,Gxyy,Gxyz,Gxzz; /* Christoffels */ CCTK_REAL Gyxx,Gyxy,Gyxz,Gyyy,Gyyz,Gyzz; CCTK_REAL Gzxx,Gzxy,Gzxz,Gzyy,Gzyz,Gzzz; @@ -136,14 +145,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, int verbose; /* Is the solver verbose */ int debug; /* Is the solver debug-verbose */ int octant; /* Apply octant BCs inside */ - int conformal; /* Do we have conformal metric ? */ - int nabla_form; /* Which form of the nable */ + int conformal=0; /* Do we have conformal metric ? */ + int nabla_form=0; /* Which form of the nable */ int matnormalize; /* Normalize the central mat value to one? */ CCTK_REAL ac; /* Storage for a(0,0,0) for renorm */ int PetscTolStyle; /* For the upper metric form of nabla */ - CCTK_REAL *uxx3, *uxy3, *uxz3, *uyy3, *uyz3, *uzz3; + CCTK_REAL *uxx3=NULL, *uxy3=NULL, *uxz3=NULL, *uyy3=NULL, *uyz3=NULL, *uzz3=NULL; /* Pointers to the data of : petsc workspace/ the GF to solve / @@ -159,6 +168,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, int Mstorage=0, Nstorage=0; + RelTol = RelTol; + octant = CCTK_Equals(domain,"octant"); verbose = CCTK_Equals(petsc_verbose,"yes")|| CCTK_Equals(petsc_verbose,"debug"); @@ -319,7 +330,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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"); + retcode = CCTK_SyncGroup(GH,"ellpetsc::petscworkspace"); if (retcode<0) CCTK_WARN(1,"Synchronization failed\n"); @@ -942,19 +953,19 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, if (PetscTolStyle == ELLCONV_ABSOLUTE) { rtol = 1.0e-15; - atol = tolerance; + atolerance = tolerance; } else if (PetscTolStyle == ELLCONV_RELATIVE) { rtol = tolerance; - atol = 1.0e-15; + atolerance = 1.0e-15; } else if (PetscTolStyle == ELLCONV_EITHER) { rtol = tolerance; - atol = tolerance; + atolerance = tolerance; } else { printf("PETSC Solver: PetscTolStyle set incorrectly [%d]\n", PetscTolStyle); } - ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT, + ierr = KSPSetTolerances(ksp,rtol,atolerance,PETSC_DEFAULT, PETSC_DEFAULT); CHKERRQ(ierr); @@ -995,7 +1006,11 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ if (verbose) CCTK_INFO("KSPSetInitialGuess\n"); +#if PETSC_VERSION_MAJOR < 2 || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 0) ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr); +#else + ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr); +#endif /* Set runtime options. Since we don't use PETSC runtime options @@ -1050,7 +1065,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, VecRestoreArray(soln,&values); /* Sync var in CCTK speak. */ - PUGH_SyncGroup(GH,"ellpetsc::petscworkspace"); + CCTK_SyncGroup(GH,"ellpetsc::petscworkspace"); /* And finally free up the matrix memory FIXME ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr); */ @@ -1067,5 +1082,6 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, if (verbose) CCTK_INFO("LEAVING ELLPETSC"); + return (0); } |