From e7d7ab5dbd67565bb5b568ad648998adb080e937 Mon Sep 17 00:00:00 2001 From: allen Date: Sat, 12 May 2001 18:42:19 +0000 Subject: Working on SOR. Seems to work, but get different answers with different numbers of processors, so can't commit a testsuite. git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@71 fa3da13c-9f13-4301-a575-cf5b8c5e1907 --- src/Startup.c | 8 +++- src/sor_confmetric.c | 34 ++++----------- src/sor_flat.c | 119 ++++++++++++++++++++++++++++++--------------------- src/sor_wrapper.c | 4 ++ 4 files changed, 89 insertions(+), 76 deletions(-) diff --git a/src/Startup.c b/src/Startup.c index 80449f3..dfe7d5f 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -10,6 +10,10 @@ #include "CactusElliptic/EllBase/src/EllBase.h" #include "CactusElliptic/EllBase/src/Ell_DBstructure.h" +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_Startup_c) + void EllSOR_Register(CCTK_ARGUMENTS); void sor_confmetric(cGH *GH, int *MetricPsiI, @@ -55,8 +59,8 @@ void EllSOR_Register(CCTK_ARGUMENTS) the routines that sets them ! */ /* Register boundary SOR can handle */ - err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllSOR::Bnd::Robin"); - err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllSOR::Bnd::Const"); + err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Robin"); + err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Const"); diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c index 6fb677b..15c4a9d 100644 --- a/src/sor_confmetric.c +++ b/src/sor_confmetric.c @@ -21,6 +21,10 @@ #include "Ell_DBstructure.h" #include "EllBase.h" +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_confmetric_c) + #define SQR(a) ((a)*(a)) int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, @@ -167,18 +171,10 @@ int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, /* 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 */ - timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } + timelevel = 0; var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); - timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[0]) - 1; - if (timelevel > 0) - { - timelevel--; - } + timelevel = 0; gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[0]); gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[1]); gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[2]); @@ -188,11 +184,7 @@ int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, if (conformal) { - timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[6]) - 1; - if (timelevel > 0) - { - timelevel--; - } + timelevel = 0; psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[6]); } @@ -201,21 +193,13 @@ int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, set flag Mstorage=1, dito for N */ if (MIndex>=0) { - timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } + timelevel = 0; Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex); Mstorage = 1; } if (NIndex>=0) { - timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1; - if (timelevel > 0) - { - timelevel--; - } + timelevel = 0; Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex); Nstorage = 1; } 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;kcctk_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;kcctk_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; } diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c index 062a846..3df20c5 100644 --- a/src/sor_wrapper.c +++ b/src/sor_wrapper.c @@ -24,6 +24,10 @@ #include "CactusElliptic/EllBase/src/EllBase.h" +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_sor_wrapper_c) + int sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, -- cgit v1.2.3