From d2890c2f9cba7aee80a2e72992c08d5041e85af5 Mon Sep 17 00:00:00 2001 From: lanfer Date: Wed, 8 Sep 1999 13:06:51 +0000 Subject: on the way to elliptic with PETSc git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@4 1d96b42b-98df-4a6a-9d84-1b24288d4588 --- src/Startup.c | 4 +++ src/make.configuration.defn | 73 +++++++++++++++++++++++++++++++++++++------ src/petsc_confmetric_solver.c | 55 +++++++++++++++++++------------- src/petsc_wrapper.c | 22 +++++++++++-- 4 files changed, 119 insertions(+), 35 deletions(-) diff --git a/src/Startup.c b/src/Startup.c index 05c3336..62ab8a9 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -8,6 +8,10 @@ #include "CactusElliptic/LinearElliptic/src/LinearElliptic.h" +/* 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 *TolArray); diff --git a/src/make.configuration.defn b/src/make.configuration.defn index abb7803..5e11f7e 100644 --- a/src/make.configuration.defn +++ b/src/make.configuration.defn @@ -1,20 +1,73 @@ # Main make.code.defn file for thorn PETSc_Elliptic # $Header$ + +# on all systems: PETSc may have been build with X11 support, which +# needs to included here. + +# on linux system: the lapack libs (in rpms) are compiled with f2c/g77 +# and need the -lg2c or -lf2c to be linked +# keep them at the end of lapack + + +### check if PETSC_DIR/PETSC_ARCH are set, bail out if not + +ifeq ($(strip $(PETSC_DIR)), ) +$(NAME): MissingPETSC_DIR +.pseudo: MissingPETSC_DIR +MissingPETSC_DIR: + @echo "PETSc: need environment variable PETSC_DIR for compliling EllPETSc" + @echo "PETSc: set PETSC_DIR or remove EllPETSc from Thornlist" + exit 2 +endif + +ifeq ($(strip $(PETSC_ARCH)), ) +$(NAME): MissingPETSC_ARCH +.pseudo: MissingPETSC_ARCH +MissingPETSC_ARCH: + @echo "PETSc: need environment variable PETSC_ARCH for compliling EllPETSc" + @echo "PETSc: set PETSC_ARCH or remove EllPETSc from Thornlist" + exit 2 +endif + + + LIBDIRS += $(PETSC_DIR)/lib/libO/$(PETSC_ARCH) -#on linux system: the lapack libs (in rpms) are compiled with f2c/g77 -#and need the -lg2c or -lf2c to be linked -#keep them at the end of lapack +### PETSC on the T3E +ifeq ($(PETSC_ARCH),t3e) + echo "No T3E support" + exit 2 + MYLIBS := petscts petscsnes petscsles \ + petscmat petscvec petscsys X11 sci + + LIBS = $(MYLIBS) + +endif + +### PETSC on the SGI +ifeq ($(PETSC_ARCH),sgi) + echo "No SGI support" + exit 2 + MYLIBS := petscts petscsnes petscsles \ + petscmat petscvec petscsys X11 \ + complib.sgimath + + LIBS = $(MYLIBS) + +endif -MYLIBS := petsc petscdm petsccontrib \ - petscts petscsnes petscsles \ - petscmat petscvec \ - blas lapack g2c \ - mpich $(LIBS) +### PETSC on linux +ifeq ($(PETSC_ARCH),linux) + MYLIBS := petsc petscdm petsccontrib \ + petscts petscsnes petscsles \ + petscmat petscvec \ + blas lapack g2c \ + mpich $(LIBS) -LIBS = $(MYLIBS) + LIBS = $(MYLIBS) -EXTRAFLAGS = -I$(PETSC_DIR)/bmake/linux/base + EXTRAFLAGS = -I$(PETSC_DIR)/bmake/$(PETSC_ARCH)/base +endif diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c index c60312e..9a2e9d7 100644 --- a/src/petsc_confmetric_solver.c +++ b/src/petsc_confmetric_solver.c @@ -150,11 +150,12 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, CCTK_REAL *psi =NULL; CCTK_REAL *psix=NULL, *psiy=NULL, *psiz=NULL; - + int Mstorage=0, Nstorage=0; octant = CCTK_Equals(domain,"octant"); - verbose = CCTK_Equals(petsc_verbose,"yes")||CCTK_Equals(petsc_verbose,"debug"); + verbose = CCTK_Equals(petsc_verbose,"yes")|| + CCTK_Equals(petsc_verbose,"debug"); debug = CCTK_Equals(petsc_verbose,"debug"); reuse = 0; matnormalize = 0; @@ -213,9 +214,6 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, /* derive the metric data pointer from the index array. Note the ordering in the metric */ ell_field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex); - Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*MIndex); - Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex); - gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[0]); gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[1]); gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[2]); @@ -230,9 +228,20 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, psiy =(CCTK_REAL*) GetDataPtr_NextTL(GH,"einstein::psiy"); psiz =(CCTK_REAL*) GetDataPtr_NextTL(GH,"einstein::psiz"); } + /* 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) { + Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*MIndex); + Mstorage = 1; + } + if (*NIndex>=0) { + Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex); + Nstorage = 1; + } + /* Get the workspace data pointer */ wsp = GetDataPtr_NextTL(GH,"ellpetsc::wsp"); - /* initialize the linear index lookup table, after it's filled up (below), the -1 indicates a boundary */ @@ -261,11 +270,11 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, kmax=((pughGH->neighbors[pughGH->myproc][ZDP]<0) ? pughGH->lnz-1 : pughGH->lnz-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 (inpoints*sizeof(CCTK_REAL)); - uxy3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uxz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uyy3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uyz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); - uzz3 = malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uxx3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uxy3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uxz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uyy3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uyz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); + uzz3 = (CCTK_REAL*)malloc(pughGH->npoints*sizeof(CCTK_REAL)); for (i=0;inpoints;i++) { CCTK_REAL det; CCTK_REAL p12; @@ -436,8 +445,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, det = det*p12; /* Rescaling for the uppermetric solver */ - Mlin[i] = Mlin[i]*sqrt(det); - Nlin[i] = Nlin[i]*sqrt(det); + 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); @@ -736,7 +745,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } /* end nabla_form=3 */ /* M phi */ - 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 @@ -818,8 +828,8 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, } } - rhsval = -rhsval + Nlin[ijk] / ac; - /* printf("RHS(%d %d %d)+N: %f",i,j,k,rhsval);*/ + if (Nstorage) + rhsval = -rhsval + Nlin[ijk] / ac; ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); CHKERRA(ierr); @@ -1062,3 +1072,4 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize, if (verbose) CCTK_INFO("LEAVING ELLPETSC"); } + diff --git a/src/petsc_wrapper.c b/src/petsc_wrapper.c index 8cc39a4..15407c6 100644 --- a/src/petsc_wrapper.c +++ b/src/petsc_wrapper.c @@ -3,21 +3,35 @@ #include "cctk_parameters.h" #include "cctk_Flesh.h" +/* The wrapper functions for the core PETSc solver, that performs the + actual solve. + One wrapper fucntion 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 */ void petsc_metric(cGH *GH, int *MetricI, int *FieldIndex, int *MIndex, int *NIndex, int *AbsTol, int *RelTol) { + /* Arrays for the M/N size info */ + void petsc_confmetric_solver(cGH *GH, int *MetricI, int MetricISize, int *FieldIndex, int *MIndex, int *NIndex, int *AbsTol, int *RelTol); + + + 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, + petsc_confmetric_solver(GH, MetricI, 6, FieldIndex, + MIndex, NIndex, AbsTol, RelTol); } - +/* wrapper function for the class of ellitpic equations, that needs a conf. + factor */ void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldIndex, int *MIndex, int *NIndex, int *AbsTol, int *RelTol) { @@ -26,6 +40,8 @@ void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldIndex, int *FieldIndex, int *MIndex, int *NIndex, int *AbsTol, int *RelTol); + /* petsc_confmetric_solver expects an metric array, in the case of this + equation class, needs conf.factor as last entry -> size 7 */ petsc_confmetric_solver(GH, MetricPsiI, 7, FieldIndex, MIndex, NIndex, AbsTol, RelTol); -- cgit v1.2.3