From 8a173b605446372ef742d879610f43793d58a891 Mon Sep 17 00:00:00 2001 From: allen Date: Thu, 3 Jan 2002 20:44:52 +0000 Subject: Changes for PETSc 2.1.0 git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@61 1d96b42b-98df-4a6a-9d84-1b24288d4588 --- src/Startup.c | 31 ++---- src/ellpetsc.h | 21 ++++ src/make.configuration.defn | 2 +- src/petsc_confmetric.c | 62 +++++------ src/petsc_confmetric_solver.c | 74 +++++++------ src/petsc_flat_solver.c | 251 +++++++++++++++++++++++++++--------------- src/petsc_wrapper.c | 14 ++- 7 files changed, 274 insertions(+), 181 deletions(-) create mode 100644 src/ellpetsc.h diff --git a/src/Startup.c b/src/Startup.c index 84c8b9c..24306df 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -9,38 +9,25 @@ 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); + +void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI, + int *NI, int *AbsTol, int *RelTol); + +void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, + int *AbsTol, int *RelTol); + /* Registration of the petsc solvers with the Elliptic solver registry. This routine registers petsc_confmetric under the name "petsc" for the class of elliptic equations "LinConfMetric" */ int EllPETSc_Register(cGH *GH) { - void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI, - int *NI, int *AbsTol, int *RelTol); - - void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI, - int *NI, int *AbsTol, int *RelTol); - - void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, - int *AbsTol, int *RelTol); DECLARE_CCTK_PARAMETERS - if (CCTK_Equals(elliptic_verbose,"yes")) - { - printf("PETSc: Registering petsc for Ell_LinConfMetric...\n"); - } Ell_RegisterSolver(petsc_confmetric,"petsc","Ell_LinConfMetric"); - - if (CCTK_Equals(elliptic_verbose,"yes")) - { - printf("PETSc: Registering petsc for Ell_LinMetric...\n"); - } Ell_RegisterSolver(petsc_metric,"petsc","Ell_LinMetric"); - - if (CCTK_Equals(elliptic_verbose,"yes")) - { - printf("PETSc: Registering petsc for Ell_LinFlat...\n"); - } Ell_RegisterSolver(petsc_flat,"petsc","Ell_LinFlat"); return 0; diff --git a/src/ellpetsc.h b/src/ellpetsc.h new file mode 100644 index 0000000..a1dceb2 --- /dev/null +++ b/src/ellpetsc.h @@ -0,0 +1,21 @@ +/*@@ + @header ellpetsc.h + @author Gabrielle Allen + @date 17th September 2001 + @desc + Include file for thorn EllPETSc + @enddesc + @version $Header$ +@@*/ + +#ifndef _ELLPETSC_H_ +#define _ELLPETSC_H_ 1 + +#define XDM 0 +#define XDP 1 +#define YDM 2 +#define YDP 3 +#define ZDM 4 +#define ZDP 5 + +#endif /* _ELLPETSC_H_ */ diff --git a/src/make.configuration.defn b/src/make.configuration.defn index 4fc63b8..7ea39c3 100644 --- a/src/make.configuration.defn +++ b/src/make.configuration.defn @@ -51,7 +51,7 @@ endif ### PETSC on linux ifeq ($(PETSC_ARCH),linux) - PLATFORM_LIBS = lapack blas g2c mpich + PLATFORM_LIBS = flapack fblas g2c mpich endif ### Otherwise diff --git a/src/petsc_confmetric.c b/src/petsc_confmetric.c index 96f1dbb..d85307a 100644 --- a/src/petsc_confmetric.c +++ b/src/petsc_confmetric.c @@ -15,8 +15,8 @@ #include "assert.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$"; @@ -354,16 +354,16 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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); /*$}$*/ /*$else {$*/ /* ierr = MatRestoreValuesMemory(A); */ /*$CCTK_WARN(0,"no reuse"); - CHKERRA(ierr); + CHKERRQ(ierr); }$*/ @@ -702,7 +702,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, else { J = wsp[DATINDEX(pughGH,i+l,j+m,k+n)]; ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); - CHKERRA(ierr); + CHKERRQ(ierr); } } } @@ -711,7 +711,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, /* printf("RHS(%d %d %d)+N: %f",i,j,k,rhsval);*/ ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); - CHKERRA(ierr); + CHKERRQ(ierr); } } @@ -719,18 +719,18 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } 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); + 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); OptionsSetValue("-ksp_monitor",""); CCTK_INFO("CREATING SLES"); - ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles); CHKERRA(ierr); + ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles); CHKERRQ(ierr); /* @@ -743,15 +743,15 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, means, we'll leave it for now. */ CCTK_INFO("CREATING SLES OPERATOR"); - ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr); + 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. */ CCTK_INFO("SLESGet KSP/PC"); - ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr); - ierr = SLESGetPC(sles,&pc); CHKERRA(ierr); + ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); + ierr = SLESGetPC(sles,&pc); CHKERRQ(ierr); /* Get the PC Type */ if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) @@ -772,7 +772,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n"); ierr = PCSetType(pc,PCNONE); } - CHKERRA(ierr); + CHKERRQ(ierr); /* Now the same thing for the KSP Type */ @@ -800,7 +800,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n"); ierr = KSPSetType(ksp,KSPBCGS); } - CHKERRA(ierr); + CHKERRQ(ierr); @@ -822,7 +822,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT, PETSC_DEFAULT); - CHKERRA(ierr); + CHKERRQ(ierr); /* We are warned in the manual that @@ -849,7 +849,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ /*$ierr = KSPGMRESSetOrthogonalization(ksp, KSPGMRESUnmodifiedGramSchmidtOrthogonalization); - CHKERRA(ierr);$*/ + CHKERRQ(ierr);$*/ /* We also learn that ... @@ -861,7 +861,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ CCTK_INFO("KSPSetInitialGuess\n"); - ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr); + ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr); /* Set runtime options. Since we don't use PETSC runtime options @@ -869,7 +869,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 ... */ @@ -878,7 +878,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, Solve the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ CCTK_INFO("SLES solve..."); - ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); + 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 @@ -925,14 +925,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, PUGH_SyncGroup(GH,"ellpetsc::petscworkspace"); /* And finally free up the matrix memory FIXME - 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); + ierr = SLESDestroy(sles); CHKERRQ(ierr); + ierr = VecDestroy(soln); CHKERRQ(ierr); + ierr = VecDestroy(b); CHKERRQ(ierr); + ierr = MatDestroy(A[0]); CHKERRQ(ierr); free(A); /*$PetscFinalize();$*/ diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c index 0219ed3..59441d3 100644 --- a/src/petsc_confmetric_solver.c +++ b/src/petsc_confmetric_solver.c @@ -15,8 +15,8 @@ #include "assert.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$"; @@ -200,7 +200,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, /* Set the PETSc communicator to set of PUGH and initialzie PETSc */ - 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"); @@ -371,11 +371,11 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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); @@ -819,7 +819,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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); } } } @@ -828,7 +828,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, rhsval = -rhsval - Nlin[ijk] / ac; ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); - CHKERRA(ierr); + CHKERRQ(ierr); } } @@ -836,13 +836,13 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } 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); + 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 (nabla_form == 3) { if (verbose) @@ -855,10 +855,16 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, free(uzz3); } 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); /* @@ -871,15 +877,15 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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); + 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); + ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); + ierr = SLESGetPC(sles,&pc); CHKERRQ(ierr); /* Get the PC Type */ if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) @@ -900,7 +906,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n"); ierr = PCSetType(pc,PCNONE); } - CHKERRA(ierr); + CHKERRQ(ierr); /* Now the same thing for the KSP Type */ @@ -928,7 +934,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n"); ierr = KSPSetType(ksp,KSPBCGS); } - CHKERRA(ierr); + CHKERRQ(ierr); @@ -950,7 +956,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT, PETSC_DEFAULT); - CHKERRA(ierr); + CHKERRQ(ierr); /* We are warned in the manual that @@ -977,7 +983,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ /*$ierr = KSPGMRESSetOrthogonalization(ksp, KSPGMRESUnmodifiedGramSchmidtOrthogonalization); - CHKERRA(ierr);$*/ + CHKERRQ(ierr);$*/ /* We also learn that ... @@ -989,7 +995,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, */ if (verbose) CCTK_INFO("KSPSetInitialGuess\n"); - ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr); + ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr); /* Set runtime options. Since we don't use PETSC runtime options @@ -997,7 +1003,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, 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 ... */ @@ -1006,7 +1012,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, Solve the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (verbose) CCTK_INFO("SLES solve..."); - ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); + 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 @@ -1047,14 +1053,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, PUGH_SyncGroup(GH,"ellpetsc::petscworkspace"); /* And finally free up the matrix memory FIXME - 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); + ierr = SLESDestroy(sles); CHKERRQ(ierr); + ierr = VecDestroy(soln); CHKERRQ(ierr); + ierr = VecDestroy(b); CHKERRQ(ierr); + ierr = MatDestroy(A[0]); CHKERRQ(ierr); free(A); /*$PetscFinalize();$*/ 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;inpoints;i++) wsp[i] = -1.0; + for (i=0;inpoints;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;irnsize[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= 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= 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"); + } } diff --git a/src/petsc_wrapper.c b/src/petsc_wrapper.c index 1f25b3e..74d3512 100644 --- a/src/petsc_wrapper.c +++ b/src/petsc_wrapper.c @@ -7,17 +7,17 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_wrapper_c) -/* The wrapper functions for the core PETSc solver, that performs the actual solve. - One wrapper function for each class, because these functions are registered - with a fixed set up arguments. Still, they all call petsc_confmetric_solver - as the core solver. */ +/* The wrapper functions for the core PETSc solver, that performs the + actual solve. One wrapper function for each class, because these + functions are registered with a fixed set up arguments. Still, they + all call petsc_confmetric_solver as the core solver. */ -/* Wrapper function for the class of ellitpic equations that needs a metric */ +/* 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) { - /* Arrays for the M/N size info */ +/* Arrays for the M/N size info */ void petsc_confmetric_solver(cGH *GH, int *MetricI, int MetricISize, int FieldIndex, int MIndex, int NIndex, @@ -25,7 +25,9 @@ void petsc_metric(cGH *GH, int *MetricI, int FieldIndex, if (GH->cctk_dim>3) + { CCTK_WARN(0,"This elliptic solver implementation does not do dimension>3!"); + } petsc_confmetric_solver(GH, MetricI, 6, FieldIndex, MIndex, NIndex, -- cgit v1.2.3