diff options
Diffstat (limited to 'src/sor_flat.c')
| -rw-r--r-- | src/sor_flat.c | 119 |
1 files changed, 70 insertions, 49 deletions
diff --git a/src/sor_flat.c b/src/sor_flat.c index 752277b..1b30f1d 100644 --- a/src/sor_flat.c +++ b/src/sor_flat.c @@ -21,6 +21,10 @@ #include "Symmetry.h" #include "Ell_DBstructure.h" #include "EllBase.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_flat_c) int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); @@ -48,13 +52,13 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL finf; CCTK_INT npow; CCTK_REAL tol; + CCTK_REAL zero=0.0; /* Iteration / stepping variables */ int sorit; int i,is,ie; int j,js,je; - int k,ks,ke,kstep; - int nxyz; + int k,ks,ke,step; int timelevel; /* stencil index */ @@ -70,7 +74,8 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, int Mstorage=0, Nstorage=0; static int firstcall = 1; CCTK_REAL dx2rec, dy2rec, dz2rec; - + char *ans; + /* Get the reduction handle */ sum_handle = CCTK_ReductionArrayHandle("sum"); if (sum_handle<0) @@ -81,17 +86,20 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, /* IF Robin BCs are set, prepare for a boundary call: setup stencil width and get Robin constants (set by the routine which is calling the solver interface) */ - if (CCTK_EQUALS(sor_bound,"robin")) - { - sw[0]=1; - sw[1]=1; - sw[2]=1; - - ierr = Ell_GetRealKey(&finf, "EllLinConfMetric::Bnd::Robin::inf"); - ierr = Ell_GetIntKey (&npow, "EllLinConfMetric::Bnd::Robin::falloff"); - } + sw[0]=1; + sw[1]=1; + sw[2]=1; + + Ell_GetStrKey(&ans,"EllLinFlat::Bnd::Robin"); + if (CCTK_Equals(ans,"yes")) + { + CCTK_WARN(1,"Robin boundary conditions are not working properly yet"); + ierr = Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf"); + ierr = Ell_GetIntKey (&npow, "EllLinFlat::Bnd::Robin::falloff"); + } + /* Only supports absolute tolerance */ tol = AbsTol[0]; if (CCTK_EQUALS(sor_accel,"const")) @@ -128,39 +136,32 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, on the same timelevel; if we have a negative index for M/N, this GF is not set, there for don't even look for it and flag it */ /* get the current time level */ - timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } - var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); + + var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex); if (MIndex>=0) { - timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } - Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex); + Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex); Mstorage = 1; } + else + { + CCTK_WARN(0,"Stop"); + } if (NIndex>=0) { - timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } - Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex); + Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex); Nstorage = 1; } + else + { + CCTK_WARN(0,"Stop"); + } /* Shortcuts */ dx = GH->cctk_delta_space[0]; dy = GH->cctk_delta_space[1]; dz = GH->cctk_delta_space[2]; - nxyz = GH->cctk_lsh[0]*GH->cctk_lsh[1]*GH->cctk_lsh[2]; dx2rec = 1.0/(dx*dx); dy2rec = 1.0/(dy*dy); @@ -175,12 +176,10 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec; - is = 1; - js = 1; + ie = GH->cctk_lsh[0]-1; je = GH->cctk_lsh[1]-1; ke = GH->cctk_lsh[2]-1; - kstep = 2; /* start at 1 for historic (Fortran) reasons */ for (sorit=1; sorit<=maxit; sorit++) @@ -200,21 +199,26 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, omega = 1.8; } - resnorm = 0.0; + resnorm = 0.0; - ks = (sorit%2)+1; - if (GH->cctk_lsh[2]==3) - { - ks = 2; - } + is = 1; + js = 1; + ks = 1; - for (k=ks;k<ke;k+=kstep) + /*is = ((GH->cctk_lbnd[0])%2)+(sorit%2)+1; + js = ((GH->cctk_lbnd[1])%2)+(sorit%2)+1; + ks = ((GH->cctk_lbnd[2])%2)+(sorit%2)+1;*/ + + /*step = 2;*/ + step = 1; + + for (k=ks;k<ke;k=k+step) { - for (j=js;j<je;j++) + for (j=js;j<je;j=j+step) { - for (i=is;i<ie;i++) + for (i=is;i<ie;i=i+step) { - + ac = ac_orig; ijk = CCTK_GFINDEX3D(GH,i ,j ,k ); @@ -248,7 +252,6 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, } } - /* reduction operation on processor-local residual values */ ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle, &resnorm, &glob_residual, CCTK_VARIABLE_REAL); if (ierr<0) @@ -256,6 +259,10 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_WARN(1,"Reduction of Norm failed"); } +#ifdef DEBUG_ELLIPTIC + printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol); +#endif + glob_residual = glob_residual / (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]); @@ -263,19 +270,25 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, #ifdef DEBUG_ELLIPTIC printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol); #endif + /* CCTK_WARN(0,"stop");*/ /* apply symmetry boundary conditions within loop */ if (CartSymVI(GH,FieldIndex)<0) { - CCTK_WARN(1,"CartSymVI failed in EllSOR loop"); + CCTK_WARN(1,"sor_flat: CartSymVI failed in EllSOR loop"); break; } + BndScalarVI(GH,sw,zero,FieldIndex); /* apply boundary conditions within loop */ - if (CCTK_EQUALS(sor_bound,"robin")) - { - ierr = BndRobinVI(GH, sw, finf, npow, FieldIndex); + /* FIXME */ + /* + if (BndRobinVI(GH, sw, finf, npow, FieldIndex)<0) + { + CCTK_WARN(1,"sor_flat: BndRobinVI failed in EllSOR loop"); + break; } + */ /* synchronization of grid variable */ CCTK_SyncGroupWithVarI(GH, FieldIndex); @@ -287,6 +300,13 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, } } + + if (elliptic_verbose) + { + CCTK_VInfo("EllSOR","Required SOR tolerence was set at %f",tol); + CCTK_VInfo("EllSOR","Achieved SOR residual was %f",glob_residual); + CCTK_VInfo("EllSOR","Number of iterations was %d (max: %d)",sorit,maxit); + } if (glob_residual>tol) { @@ -294,6 +314,7 @@ int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, retval = ELL_NOCONVERGENCE; } + if (ans) free(ans); return retval; } |