diff options
Diffstat (limited to 'src/sor_flat.c')
| -rw-r--r-- | src/sor_flat.c | 145 |
1 files changed, 112 insertions, 33 deletions
diff --git a/src/sor_flat.c b/src/sor_flat.c index 67f0463..752277b 100644 --- a/src/sor_flat.c +++ b/src/sor_flat.c @@ -1,3 +1,14 @@ + /*@@ + @file sor_flat.c + @date Tue Aug 24 12:50:07 1999 + @author Gerd Lanfermann + @desc + SOR solver for 3D flat equation + @enddesc + @@*/ + +/*#define DEBUG_ELLIPTIC*/ + #include <stdlib.h> #include <stdio.h> #include <math.h> @@ -6,21 +17,25 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" -#include "CactusBase/Boundary/src/Boundary.h" -#include "CactusBase/CartGrid3D/src/Symmetry.h" -#include "CactusElliptic/EllBase/src/Ell_DBstructure.h" +#include "Boundary.h" +#include "Symmetry.h" +#include "Ell_DBstructure.h" +#include "EllBase.h" -void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, +int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); -void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, +int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { DECLARE_CCTK_PARAMETERS + int retval = ELL_SUCCESS; + /* The pointer to the data fields */ - CCTK_REAL *Mlin=NULL, *Nlin=NULL; + CCTK_REAL *Mlin=NULL; + CCTK_REAL *Nlin=NULL; CCTK_REAL *var =NULL; /* shortcuts for deltas,etc. */ @@ -40,6 +55,7 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, int j,js,je; int k,ks,ke,kstep; int nxyz; + int timelevel; /* stencil index */ int ijk; @@ -55,55 +71,88 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, static int firstcall = 1; CCTK_REAL dx2rec, dy2rec, dz2rec; - /* Get the reduction handle */ sum_handle = CCTK_ReductionArrayHandle("sum"); if (sum_handle<0) - CCTK_WARN(1,"Cannot get reduction handle for operation >sum<"); - + { + CCTK_WARN(1,"sor_flat_3d: Cannot get reduction handle for sum operation"); + } + /* 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")) { + 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"); } /* Only supports absolute tolerance */ tol = AbsTol[0]; if (CCTK_EQUALS(sor_accel,"const")) + { accel_const = 1; + } else if (CCTK_EQUALS(sor_accel,"cheb")) + { accel_cheb = 1; + } /* Things to do only once! */ - if (firstcall==1) { + if (firstcall==1) + { if (CCTK_Equals(elliptic_verbose, "yes")) + { + if (accel_cheb) + { + CCTK_INFO("SOR with Chebyshev acceleration with radius of 1"); + } + else if (accel_const) + { + CCTK_INFO("SOR with hardcoded omega = 1.8"); + } + else { - if (accel_cheb) - printf("SOR with Chebyshev acceleration with radius of 1\n"); - else if (accel_const) - printf("SOR with hardcoded omega = 1.8\n"); - else - printf("SOR with unaccelearted relaxation (omega = 1)\n"); + CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)"); } + } firstcall = 0; } /* Get the data ptr of these GFs, They all have to be 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 */ - var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex); - if (MIndex>=0) { - Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex); + /* get the current time level */ + timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1; + if (timelevel > 0) + { + timelevel--; + } + var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); + + if (MIndex>=0) + { + timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1; + if (timelevel > 0) + { + timelevel--; + } + Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex); Mstorage = 1; } - if (NIndex>=0) { - Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex); + if (NIndex>=0) + { + timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1; + if (timelevel > 0) + { + timelevel--; + } + Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex); Nstorage = 1; } @@ -134,26 +183,37 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, kstep = 2; /* start at 1 for historic (Fortran) reasons */ - for (sorit=1; sorit<=maxit; sorit++) { - + for (sorit=1; sorit<=maxit; sorit++) + { omega = 1.0; rjacobian = 1.0; if (accel_cheb) + { for (chebit=2;chebit<sorit;chebit++) + { omega = 1.0/(1.0 - 0.25*rjacobian*rjacobian*omega); + } + } if (accel_const) + { omega = 1.8; + } resnorm = 0.0; ks = (sorit%2)+1; if (GH->cctk_lsh[2]==3) + { ks = 2; - - for (k=ks;k<ke;k+=kstep) { - for (j=js;j<je;j++) { - for (i=is;i<ie;i++) { + } + + for (k=ks;k<ke;k+=kstep) + { + for (j=js;j<je;j++) + { + for (i=is;i<ie;i++) + { ac = ac_orig; @@ -166,7 +226,9 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, ijkm = CCTK_GFINDEX3D(GH,i ,j ,k-1); if (Mstorage) + { ac += Mlin[ijk]; + } residual = ac * var[ijk] + ae *var[ipjk] + aw*var[imjk] @@ -174,14 +236,14 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, + at *var[ijkp] + ab*var[ijkm]; if (Nstorage) + { residual +=Nlin[ijk]; + } resnorm = resnorm + fabs(residual); var[ijk] = var[ijk] - omega*residual/ac; - printf(" %d %d %d %f \n",i,j,k,var[ijk]); - } } } @@ -190,32 +252,49 @@ void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle, &resnorm, &glob_residual, CCTK_VARIABLE_REAL); if (ierr<0) + { CCTK_WARN(1,"Reduction of Norm failed"); + } glob_residual = glob_residual / (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]); + +#ifdef DEBUG_ELLIPTIC + printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol); +#endif + /* apply symmetry boundary conditions within loop */ - if (CartSymVI(GH,FieldIndex)<0) + if (CartSymVI(GH,FieldIndex)<0) + { CCTK_WARN(1,"CartSymVI failed in EllSOR loop"); + break; + } /* apply boundary conditions within loop */ if (CCTK_EQUALS(sor_bound,"robin")) + { ierr = BndRobinVI(GH, sw, finf, npow, FieldIndex); + } /* synchronization of grid variable */ CCTK_SyncGroupWithVarI(GH, FieldIndex); /* Leave iteration loop if tolerance criterium is met */ if (glob_residual<tol) + { break; + } } if (glob_residual>tol) - CCTK_WARN(2,"SOR SOLVER DID NOT CONVERGE"); + { + CCTK_WARN(1,"sor_flat: SOR Solver did not converge"); + retval = ELL_NOCONVERGENCE; + } - return; + return retval; } |