on the way to elliptic with PETSc

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@4 1d96b42b-98df-4a6a-9d84-1b24288d4588

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 (i<myproc), 
-     so we count them here. (that's minus the ghostszones (we use PUGH's "rn" 
-     variables (number of gp on proc[#] in direction [i]) and subtract the
-     ghostzones if nec.) */
+  /* 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
+     (i<myproc), so we count them here (that's minus the ghostszones
+     (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++) {
     
@@ -402,12 +411,12 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
       if (verbose) 
 	CCTK_INFO("Forming nabla with upper g and finite difference of dg \n");
       nabla_form = 3;
-      uxx3 = malloc(pughGH->npoints*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;i<pughGH->npoints;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);