From fb99316260dcc0cc232927e8be09abaa1f319a8a Mon Sep 17 00:00:00 2001 From: lanfer Date: Wed, 15 Sep 1999 09:26:19 +0000 Subject: less beta git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@10 1d96b42b-98df-4a6a-9d84-1b24288d4588 --- src/petsc_flat_solver.c | 255 ++++++++++++++++-------------------------------- 1 file changed, 84 insertions(+), 171 deletions(-) diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c index 890289a..1bbefc0 100644 --- a/src/petsc_flat_solver.c +++ b/src/petsc_flat_solver.c @@ -1,15 +1,12 @@ #include "cctk.h" #include "cctk_parameters.h" -#include "cctk_Flesh.h" -#include "assert.h" #include "mpi.h" #include "sles.h" #include "CactusPUGH/PUGH/src/include/pugh.h" - #define DEBUG /*Don't know what these actually mean! FIXME */ @@ -17,47 +14,36 @@ #define ELLCONV_RELATIVE 1 #define ELLCONV_EITHER 2 -/* A few convenient macros */ +/* A convenient macros for the stencil molecule*/ #define a(i,j,k) a[i+1 + 3*(j+1 + 3*(k+1))] -#define SQR(a) ((a)*(a)) - -/* Place the matrix in global space so we can keep it around - with the memory stripped out. Thanks, Barry! */ - static int trips=0; static Mat *A; /* linear system matrix */ static Vec soln, b; /* approx solution, RHS */ 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 */ - - void *GetDataPtr_NextTL(cGH *GH, const char *field); - int CCTK_GetArgc(void); /* PETSc wnats to know the prog. name */ - void CCTK_GetArgv(char **argv); /* and parameters, don't think they can be passed through */ - - PC pc; /* preconditioner context */ KSP ksp; /* Krylov subspace method context */ int num_A; /* Number of A-arrays as needed ny Dan's MG */ - double norm; /* norm of solution error */ int ierr; /* Check the return status of petsc */ int retcode; /* Check the return status of CCTK */ - double a[27]; /* The stencil array */ - int idx[27]; /* The column of the stencil array */ - int rank; /* Rank of the matrix/vector */ - int its; /* Number of iterations */ + int idx[27]; /* The column of the stencil array */ + int rank; /* Rank of the matrix/vector */ + int its; /* Number of iterations */ + + CCTK_REAL norm; /* norm of solution error */ + CCTK_REAL a[27]; /* The stencil array */ + CCTK_REAL ac; /* Storage for a(0,0,0) for renorm */ /* Loopers */ int i,j,k,l,m,n,ll; @@ -65,25 +51,18 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, /* loop limits put in for stencil_w !=1 Ed Evans 1/19/98 */ int imin,imax,jmin,jmax,kmin,kmax,interior; - pGH *pughGH; /* The pugh Extension handle */ - int myproc; /* out processor */ - + pGH *pughGH; /* The pugh Extension handle */ + /* Tolerances */ CCTK_REAL rtol, atol, tolerance; - + /* Values to assemble the matrix */ - CCTK_REAL two=2.0, four=4.0, tmp; - CCTK_REAL pm4, psixp,psiyp,psizp,det; + CCTK_REAL tmp, *values; + + /* Misc */ CCTK_REAL dx,dy,dz; - CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz; - 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; - CCTK_REAL dxxx,dxxy,dxxz,dxyy,dxyz,dxzz; - CCTK_REAL dyxx,dyxy,dyxz,dyyy,dyyz,dyzz; - CCTK_REAL dzxx,dzxy,dzxz,dzyy,dzyz,dzzz; - CCTK_REAL *values; - + + int myproc; int nxs,nys,nzs; /* Size of the grid with stencils off... */ int startpoint=0; /* My starting index (per proc) */ int endpoint; /* One more than my end */ @@ -96,24 +75,25 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, int nabla_form; /* Which form of the nable */ int reuse; /* Do PETSc reuse tricks? */ int matnormalize; /* Normalize the central mat value to one? */ - CCTK_REAL ac; /* Storage for a(0,0,0) for renorm */ int pctype, ksptype; /* Which PC and KSP to use. See below for monster nested if statement. */ + int I,J; /* The PETSc col and row in the matrix */ + int PetscTolStyle; - /* For the upper metric form of nabla */ - CCTK_REAL *uxx3, *uxy3, *uxz3, *uyy3, *uyz3, *uzz3; - /* Pointers to the data of : petsc workspace/ the GF to solve / - metric / psi / derivs(psi) / Mlinear / Nlinear(source) */ + Mlinear / Nlinear(source) */ CCTK_REAL *wsp =NULL, *ell_field=NULL; - CCTK_REAL *Mlin=NULL, *Nlin=NULL; + /* 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"); @@ -121,6 +101,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, reuse = 0; matnormalize = 0; + /* FIXME, the TolAbs/TolRel will be evaluated here */ PetscTolStyle=0; tolerance =AbsTol[0]; @@ -149,15 +130,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, trips++; + /* Allocate storage for system matrix, will have more of them later */ num_A = 1; A=(Mat*) malloc (sizeof(Mat)*num_A); - - /* Get the data ptr of these GFs, They all have to be - on the same timelevel */ - /* derive the metric data pointer from the index array. - Note the ordering in the metric */ - ell_field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex); + + /* Get the data pointers */ /* if we have a negative index, this GF is not needed, there for don't even look for it. when index positive, @@ -170,10 +148,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex); Nstorage = 1; } - wsp = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,CCTK_VarIndex("ellpetsc::wsp")); + ell_field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex); + wsp = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,CCTK_VarIndex("ellpetsc::wsp")); - /* initialize the linear index lookup table, after it's - filled up (below), the -1 indicates a boundary */ + + /* initialize the linear index lookup table with -1; + this will indicate a boudary later */ for (i=0;inpoints;i++) wsp[i] = -1.0; /* Get myproc from the CCTK, get the gridspacing */ @@ -238,14 +218,14 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, for (j=jmin;jPUGH_COMM_WORLD, (endpoint-startpoint), /* # of rows */ @@ -302,6 +283,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, CHKERRA(ierr); ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b); CHKERRA(ierr); + /*$}$*/ /*$else {$*/ /* ierr = MatRestoreValuesMemory(A); */ @@ -310,16 +292,27 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, }$*/ - - /* Compare the PETSc layout to Cactus, this better be a match */ ierr = VecGetOwnershipRange(soln,&pvstart,&pvend); ierr = MatGetOwnershipRange(A[0],&pstart,&pend); - 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 (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"); + + + /* Initialize the stencil array */ + for (I=0;I<27;I++) a[I] = 0.0; + a( 0, 0, 0) = -10.38; + a( 1, 0, 0) = 1.56; /* ae */ + a(-1, 0, 0) = 1.56; /* aw */ + a( 0, 1, 0) = 1.56; /* an */ + a( 0,-1, 0) = 1.56; /* as */ + a( 0, 0, 1) = 1.56; /* at */ + a( 0, 0,-1) = 1.56; /* ab */ for (k=kmin;k= 0) { /* Set up indices */ - int ijk; /* The data point for the array */ + int ijk; /* The lineax index for the array */ int ig, jg, kg; /* The position in global space */ - int I,J; /* The col and row in the matrix */ CCTK_REAL rhsval = 0; - - for (I=0;I<27;I++) a[I] = 0.0; - a( 0, 0, 0) = -10.38; - a( 1, 0, 0) = 1.56; /* ae */ - a(-1, 0, 0) = 1.56; /* aw */ - a( 0, 1, 0) = 1.56; /* an */ - a( 0,-1, 0) = 1.56; /* as */ - a( 0, 0, 1) = -1.56; /* at */ - a( 0, 0,-1) = -1.56; /* ab */ - - /* these guys are easy */ - ijk = DATINDEX(pughGH,i,j,k); /* get linear index */ + + /* local linear index / global 3d index */ + ijk = DATINDEX(pughGH,i,j,k); ig = i + GH->cctk_lbnd[0]; jg = j + GH->cctk_lbnd[1]; kg = k + GH->cctk_lbnd[2]; @@ -353,45 +336,19 @@ 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. This is - really painful due to the boundaries (here we force - dirichlet). */ + /* Great now set the values of the matrix. Here we force dirichlet. */ VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES); - /* Put in the octant boundary conditions. - - We do these by reflecting the -1 stencil if we are at a - boundary at -1. That is, imagine the 1D laplace equation - with the boundary condition a(0) = a(1). This means - the point 1 stencil (which is our first stencil in - the matrix) becomes - - a(0) + a(2) - 2 a(1) = rho(i) * deltax -> - a(2) - 2 a(1) + a(1) = rho(i) * deltax - - OK so think about this with a general stenci - - S_0 a_0 + S_1 a_1 + S_2 a_2 = rho_1 -> - S_0 a_1 + S_1 a_1 + S_2 a_2 = rho_1 -> - (S_0 + S_1) a_1 + S_2 a_2 = rho_1 - - eg, S_0 -> 0 and S_1 -> S_1 + S_0 - - Great, so implement this in 3D. It is a bit trickier, - since we don't want to zero neighbors in the wrong - direction, but not impossible. - - */ - +#ifdef DEBUG 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)); printf("\n"); } +#endif if (octant) for (l=-1;l<=0;l++) @@ -417,7 +374,8 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, ac = a(0,0,0); for (J=0;J<27;J++) a[J] = a[J] / ac; } - else ac = 1; + else + ac = 1; /* This is the new way-clever look. Note it relies @@ -442,17 +400,19 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, } } } - - if (Nstorage) - rhsval = -rhsval + Nlin[ijk] / ac; + + if (Nstorage) rhsval = -rhsval + Nlin[ijk] / ac; ierr = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES); CHKERRA(ierr); - } - } - } - } + /* reset the central value*/ + a( 0, 0, 0) = -10.38; + + } /* if (wsp>=0) */ + } /*for i */ + } /* for j */ + } /* for k */ if (verbose) CCTK_INFO ("Assembling the vectors"); ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr); @@ -462,21 +422,11 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, ierr = VecAssemblyEnd(soln); CHKERRA(ierr); ierr = VecAssemblyEnd(b); CHKERRA(ierr); ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(ierr); - - if (nabla_form == 3) { - if (verbose) - printf ("Freeing upper metric storage\n"); - free(uxx3); - free(uxy3); - free(uxz3); - free(uyy3); - free(uyz3); - free(uzz3); - } + if (trips==0) OptionsSetValue("-ksp_monitor",""); - CCTK_INFO("CREATING SLES"); + if (verbose) CCTK_INFO("CREATING SLES"); ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles); CHKERRA(ierr); @@ -489,14 +439,14 @@ 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. */ - CCTK_INFO("CREATING SLES OPERATOR"); + if (verbose) CCTK_INFO("CREATING SLES OPERATOR"); ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(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"); + if (verbose) CCTK_INFO("SLESGet KSP/PC"); ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr); ierr = SLESGetPC(sles,&pc); CHKERRA(ierr); @@ -567,48 +517,12 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, PetscTolStyle); } - ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT, - PETSC_DEFAULT); + ierr = KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, PETSC_DEFAULT); CHKERRA(ierr); - - /* We are warned in the manual that - - The default technique for orthogonalization of the Hessenberg matrix - in GMRES is the modified Gram-Schmidt method, which employs many - VecDot() operations and can thus be slow in parallel. A fast approach - is to use the unmodified Gram-Schmidt method, which can be set with - - ierr = KSPGMRESSetOrthogonalization(KSP ksp, - KSPGMRESUnmodifiedGramSchmidtOrthogonalization); - - or the options database command - -ksp_gmres_unmodifiedgramschmidt. Note that this algorithm is - numerically unstable, but may deliver much better speed - performance. One can also use unmodifed Gram-Schmidt with - iterative refinement, by setting the orthogonalization routine, - KSPGMRESIROrthog(), by using the command line option - -ksp_gmres_irorthog. - - So this my not be a smart thing to do. But lets put it here - so we can turn it on or off later. - - */ - /*$ierr = KSPGMRESSetOrthogonalization(ksp, - KSPGMRESUnmodifiedGramSchmidtOrthogonalization); - CHKERRA(ierr);$*/ - - /* We also learn that ... - - By default, KSP assumes an initial guess of zero by zeroing the - initial value for the solution vector that is given. To use a - nonzero initial guess, the user must call - - ierr = KSPSetInitialGuessNonzero(KSP ksp); - - */ - CCTK_INFO("KSPSetInitialGuess\n"); - ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr); + if (verbose) CCTK_INFO("KSPSetInitialGuess\n"); + ierr = KSPSetInitialGuessNonzero(ksp); + CHKERRA(ierr); /* Set runtime options. Since we don't use PETSC runtime options @@ -624,7 +538,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - CCTK_INFO("SLES solve..."); + if (verbose) CCTK_INFO("SLES solve..."); ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); /* Here we can form a "res = Ax - b" and find its norm to get @@ -659,7 +573,7 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, /* Now this one. "Fortran-order" the matrix. But remember we have stripped off the ghost zones. Hence ig-1 and nxs... */ - I =wsp[ijk]-startpoint; + I = wsp[ijk]-startpoint; ell_field[ijk] = values[I]; } } @@ -668,11 +582,10 @@ void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, VecRestoreArray(soln,&values); - /* Sync var in CCTK speak. */ pugh_SyncGroup(GH,"ellpetsc::petscworkspace"); - /* And finally free up the matrix memory FIXME - ierr = MatReleaseValuesMemory(A); CHKERRA(ierr); */ + /* And finally free up the matrix memory */ + /*$ierr = MatReleaseValuesMemory(A); CHKERRA(ierr);$*/ /* This code is not used anymore */ if (verbose) CCTK_INFO("Destroying Matrices"); -- cgit v1.2.3