diff options
author | allen <allen@fa3da13c-9f13-4301-a575-cf5b8c5e1907> | 2001-04-13 14:30:23 +0000 |
---|---|---|
committer | allen <allen@fa3da13c-9f13-4301-a575-cf5b8c5e1907> | 2001-04-13 14:30:23 +0000 |
commit | d8f42b55c08565f64087d7c1fe54cff04207ddde (patch) | |
tree | 6654a2cf1d6153bc86e5654f7981002a0e7bbb5c | |
parent | 371749064f2bd37f05898ee66128c72e5c6bc2d4 (diff) |
Fixed a number of problems, notably that timelevels weren't being used properly, and
error codes not returned.
Still some testing to do, then I'll add a testsuite for SOR.
git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@69 fa3da13c-9f13-4301-a575-cf5b8c5e1907
-rw-r--r-- | src/sor_confmetric.c | 310 | ||||
-rw-r--r-- | src/sor_flat.c | 145 | ||||
-rw-r--r-- | src/sor_wrapper.c | 76 |
3 files changed, 346 insertions, 185 deletions
diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c index 4f5a2e4..6fb677b 100644 --- a/src/sor_confmetric.c +++ b/src/sor_confmetric.c @@ -16,13 +16,14 @@ #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" #define SQR(a) ((a)*(a)) -void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, +int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); @@ -42,14 +43,15 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, @@*/ - - -void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, +int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { DECLARE_CCTK_PARAMETERS + int retval = ELL_SUCCESS; + int timelevel; + /* The pointer to the metric fields */ CCTK_REAL *gxx =NULL, *gxy =NULL; CCTK_REAL *gxz =NULL, *gyy =NULL; @@ -112,12 +114,15 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, /* Get the reduction handle */ sum_handle = CCTK_ReductionArrayHandle("sum"); if (sum_handle<0) + { CCTK_WARN(1,"Cannot get reduction handle for operation >sum<"); - + } + /* 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; @@ -129,21 +134,32 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, /* 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) + { + printf("SOR with Chebyshev acceleration with radius of 1\n"); + } + else if (accel_const) + { + printf("SOR with hardcoded omega = 1.8\n"); + } + 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"); + printf("SOR with unaccelearted relaxation (omega = 1)\n"); } + } firstcall = 0; } @@ -151,27 +167,56 @@ void 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 */ - var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex); - 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]); - gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[3]); - gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[4]); - gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[5]); - + timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1; + if (timelevel > 0) + { + timelevel--; + } + var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); + + timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[0]) - 1; + if (timelevel > 0) + { + timelevel--; + } + 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]); + gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[3]); + gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[4]); + gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[5]); if (conformal) - psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[6]); + { + timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[6]) - 1; + if (timelevel > 0) + { + timelevel--; + } + psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[6]); + } /* if we have a negative index for M/N, 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); + 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; } @@ -201,7 +246,8 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, /* Calculate the inverse metric */ - for (ijk=0;ijk<nxyz;ijk++) { + for (ijk=0;ijk<nxyz;ijk++) + { /* determinant */ det = -(SQR(gxz[ijk])*gyy[ijk]) + @@ -210,10 +256,13 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, SQR(gxy[ijk])*gzz[ijk] + gxx[ijk]*gyy[ijk]*gzz[ijk]; - if (conformal) { + if (conformal) + { pm4 = 1.0/pow(psi[ijk],4.0); p12 = pow(psi[ijk],12.0); - } else { + } + else + { pm4 = 1.0; p12 = 1.0; } @@ -247,8 +296,10 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, } /*$for (k=1;k<GH->cctk_lsh[2]-1;k++) { - for (j=1;j<GH->cctk_lsh[1]-1;j++) { - for (i=1;i<GH->cctk_lsh[0]-1;i++) { + for (j=1;j<GH->cctk_lsh[1]-1;j++) + { + for (i=1;i<GH->cctk_lsh[0]-1;i++) + { ijk = CCTK_GFINDEX3D(GH,i,j,k); printf("U G : %d %d %d %7.8g %7.8g %7.8g \n", i,j,k,uxx[ijk],uyy[ijk],gxx[ijk]); @@ -272,99 +323,111 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, 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 = 1; - for (k=ks;k<ke;k+=kstep) { - for (j=js;j<je;j++) { - for (i=is;i<ie;i++) { - - ijk = CCTK_GFINDEX3D(GH,i ,j ,k ); - - ipjk = CCTK_GFINDEX3D(GH,i+1,j ,k ); - imjk = CCTK_GFINDEX3D(GH,i-1,j ,k ); - ijpk = CCTK_GFINDEX3D(GH,i ,j+1,k ); - ijmk = CCTK_GFINDEX3D(GH,i ,j-1,k ); - ijkp = CCTK_GFINDEX3D(GH,i ,j ,k+1); - ijkm = CCTK_GFINDEX3D(GH,i ,j ,k-1); - - ipjpk = CCTK_GFINDEX3D(GH,i+1,j+1,k ); - ipjmk = CCTK_GFINDEX3D(GH,i+1,j-1,k ); - imjpk = CCTK_GFINDEX3D(GH,i-1,j+1,k ); - imjmk = CCTK_GFINDEX3D(GH,i-1,j-1,k ); - - ijpkp = CCTK_GFINDEX3D(GH,i ,j+1,k+1); - ijpkm = CCTK_GFINDEX3D(GH,i ,j+1,k-1); - ijmkp = CCTK_GFINDEX3D(GH,i ,j-1,k+1); - ijmkm = CCTK_GFINDEX3D(GH,i ,j-1,k-1); - - ipjkp = CCTK_GFINDEX3D(GH,i+1,j ,k+1); - ipjkm = CCTK_GFINDEX3D(GH,i+1,j ,k-1); - imjkp = CCTK_GFINDEX3D(GH,i-1,j ,k+1); - imjkm = CCTK_GFINDEX3D(GH,i-1,j ,k-1); - - - ac = -1.0*uxx[ipjk] - 2.0*uxx[ijk] - 1.0*uxx[imjk] - -1.0*uyy[ijpk] - 2.0*uyy[ijk] - 1.0*uyy[ijmk] - -1.0*uzz[ijkp] - 2.0*uzz[ijk] - 1.0*uzz[ijkm]; - - - if (Mstorage) - ac += Mlin[ijk]; - - - ae = uxx[ipjk]+uxx[ijk]; - aw = uxx[imjk]+uxx[ijk]; - an = uyy[ijpk]+uyy[ijk]; - as = uyy[ijmk]+uyy[ijk]; - at = uzz[ijkp]+uzz[ijk]; - ab = uzz[ijkm]+uzz[ijk]; - - ane = uxy[ijpk]+uxy[ipjk]; - anw =-uxy[imjk]-uxy[ijpk]; - ase =-uxy[ipjk]-uxy[ijmk]; - asw = uxy[imjk]+uxy[ijmk]; - - ate = uxz[ijkp]+uxz[ipjk]; - atw =-uxz[imjk]-uxz[ijkp]; - abe =-uxz[ipjk]-uxz[ijkm]; - abw = uxz[imjk]+uxz[ijkm]; - - atn = uyz[ijpk]+uyz[ijkp]; - ats =-uyz[ijkp]-uyz[ijmk]; - abn =-uyz[ijkm]-uyz[ijpk]; - absol = uyz[ijkm]+uyz[ijmk]; - - residual = ac * var[ijk] - + ae *var[ipjk] + aw*var[imjk] - + an *var[ijpk] + as*var[ijmk] - + at *var[ijkp] + ab*var[ijkm]; - - residual = residual - + ane*var[ipjpk] + anw*var[imjpk]; - - residual = residual - + ase*var[ipjmk] + asw*var[imjmk]; - - residual = residual - + ate*var[ipjkp] + atw*var[imjkp] - + abe*var[ipjkm] + abw*var[imjkm] - + atn*var[ijpkp] + ats*var[ijmkp] - + abn*var[ijpkm] + absol*var[ijmkm]; - - if (Nstorage) - residual +=Nlin[ijk]; + for (k=ks;k<ke;k+=kstep) + { + for (j=js;j<je;j++) + { + for (i=is;i<ie;i++) + { + + ijk = CCTK_GFINDEX3D(GH,i ,j ,k ); + + ipjk = CCTK_GFINDEX3D(GH,i+1,j ,k ); + imjk = CCTK_GFINDEX3D(GH,i-1,j ,k ); + ijpk = CCTK_GFINDEX3D(GH,i ,j+1,k ); + ijmk = CCTK_GFINDEX3D(GH,i ,j-1,k ); + ijkp = CCTK_GFINDEX3D(GH,i ,j ,k+1); + ijkm = CCTK_GFINDEX3D(GH,i ,j ,k-1); + + ipjpk = CCTK_GFINDEX3D(GH,i+1,j+1,k ); + ipjmk = CCTK_GFINDEX3D(GH,i+1,j-1,k ); + imjpk = CCTK_GFINDEX3D(GH,i-1,j+1,k ); + imjmk = CCTK_GFINDEX3D(GH,i-1,j-1,k ); + + ijpkp = CCTK_GFINDEX3D(GH,i ,j+1,k+1); + ijpkm = CCTK_GFINDEX3D(GH,i ,j+1,k-1); + ijmkp = CCTK_GFINDEX3D(GH,i ,j-1,k+1); + ijmkm = CCTK_GFINDEX3D(GH,i ,j-1,k-1); + + ipjkp = CCTK_GFINDEX3D(GH,i+1,j ,k+1); + ipjkm = CCTK_GFINDEX3D(GH,i+1,j ,k-1); + imjkp = CCTK_GFINDEX3D(GH,i-1,j ,k+1); + imjkm = CCTK_GFINDEX3D(GH,i-1,j ,k-1); + + + ac = -1.0*uxx[ipjk] - 2.0*uxx[ijk] - 1.0*uxx[imjk] + -1.0*uyy[ijpk] - 2.0*uyy[ijk] - 1.0*uyy[ijmk] + -1.0*uzz[ijkp] - 2.0*uzz[ijk] - 1.0*uzz[ijkm]; + + + if (Mstorage) + { + ac += Mlin[ijk]; + } + + ae = uxx[ipjk]+uxx[ijk]; + aw = uxx[imjk]+uxx[ijk]; + an = uyy[ijpk]+uyy[ijk]; + as = uyy[ijmk]+uyy[ijk]; + at = uzz[ijkp]+uzz[ijk]; + ab = uzz[ijkm]+uzz[ijk]; + + ane = uxy[ijpk]+uxy[ipjk]; + anw =-uxy[imjk]-uxy[ijpk]; + ase =-uxy[ipjk]-uxy[ijmk]; + asw = uxy[imjk]+uxy[ijmk]; + + ate = uxz[ijkp]+uxz[ipjk]; + atw =-uxz[imjk]-uxz[ijkp]; + abe =-uxz[ipjk]-uxz[ijkm]; + abw = uxz[imjk]+uxz[ijkm]; + + atn = uyz[ijpk]+uyz[ijkp]; + ats =-uyz[ijkp]-uyz[ijmk]; + abn =-uyz[ijkm]-uyz[ijpk]; + absol = uyz[ijkm]+uyz[ijmk]; + + residual = ac * var[ijk] + + ae *var[ipjk] + aw*var[imjk] + + an *var[ijpk] + as*var[ijmk] + + at *var[ijkp] + ab*var[ijkm]; + + residual = residual + + ane*var[ipjpk] + anw*var[imjpk]; + + residual = residual + + ase*var[ipjmk] + asw*var[imjmk]; + + residual = residual + + ate*var[ipjkp] + atw*var[imjkp] + + abe*var[ipjkm] + abw*var[imjkm] + + atn*var[ijpkp] + ats*var[ijmkp] + + abn*var[ijpkm] + absol*var[ijmkm]; + + if (Nstorage) + { + residual +=Nlin[ijk]; + } - resnorm = resnorm + fabs(residual); - var[ijk] = var[ijk] - omega*residual/ac; + resnorm = resnorm + fabs(residual); + var[ijk] = var[ijk] - omega*residual/ac; } } } @@ -373,35 +436,46 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle, &resnorm, &glob_residual, CCTK_VARIABLE_REAL); if (ierr<0) + { CCTK_WARN(1,"Reduction of residual failed"); + } glob_residual = glob_residual / (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]); /* apply symmetry boundary conditions within loop */ if (CartSymVI(GH,FieldIndex)<0) + { CCTK_WARN(1,"CartSymVI failed in EllSOR loop"); - + } + /* 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; + } } /* Information for the user if the solve did not converge within the given constraints of max.iteration and tolerance */ if (residual>tol) + { + retval = ELL_NOCONVERGENCE; CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, "SOR SOLVER DID NOT CONVERGE within given " "tolerance/max.number of iterations.\n " "max. iteration %d residual (tolerance): %g (%g)\n", maxit, glob_residual, tol ); + } if (uxx) free(uxx); if (uyy) free(uyy); @@ -410,7 +484,7 @@ void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, if (uxz) free(uxz); if (uyz) free(uyz); - return; + return retval; } 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; } diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c index 4f2a8b0..062a846 100644 --- a/src/sor_wrapper.c +++ b/src/sor_wrapper.c @@ -22,19 +22,21 @@ #include "cctk_Parameters.h" #include "cctk_FortranString.h" -void sor_flat_3d(cGH *GH, int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); -void sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal, - int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); -void sor_confmetric(cGH *GH, - int *MetricPsiI, - int FieldIndex, - int MIndex, - int NIndex, - CCTK_REAL *AbsTol, - CCTK_REAL *RelTol); -void sor_flat(cGH *GH, +#include "CactusElliptic/EllBase/src/EllBase.h" + +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, + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); +int sor_confmetric(cGH *GH, + int *MetricPsiI, + int FieldIndex, + int MIndex, + int NIndex, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol); +int sor_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, @@ -65,32 +67,36 @@ void sor_flat(cGH *GH, @@*/ -void sor_confmetric(cGH *GH, - int *MetricPsiI, - int FieldIndex, - int MIndex, - int NIndex, - CCTK_REAL *AbsTol, - CCTK_REAL *RelTol) +int sor_confmetric(cGH *GH, + int *MetricPsiI, + int FieldIndex, + int MIndex, + int NIndex, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol) { + int retval = ELL_NOSOLVER; switch (CCTK_GroupDimFromVarI(FieldIndex)) { case 1: - CCTK_WARN(0,"No 1D SOR solver implemented"); + CCTK_WARN(0,"sor_confmetric: No 1D SOR solver implemented"); break; case 2: - CCTK_WARN(0,"No 2D SOR solver implemented"); + CCTK_WARN(0,"sor_confmetric: No 2D SOR solver implemented"); break; case 3: - sor_confmetric_3d(GH, MetricPsiI, 1, FieldIndex, MIndex, NIndex, - AbsTol, RelTol); + retval = sor_confmetric_3d(GH, MetricPsiI, 1, + FieldIndex, MIndex, NIndex, + AbsTol, RelTol); break; default: - CCTK_WARN(1,"Requested SOR solver for dimension equal " + CCTK_WARN(1,"sor_confmetric: Requested SOR solver for dimension equal " "zero or gt. three not implemented!"); break; } + + return retval; } /*@@ @@ -105,9 +111,9 @@ void sor_confmetric(cGH *GH, This wrapper is registered and if it is being called with a n-dim. grid function, it goes of and picks the correct solver. - We pass in the arguments that are neccessary for this class of elliptic eq. - this solver is intended to solve. See ./CactusElliptic/EllBase/src/ for the - classes of elliptic eq. + We pass in the arguments that are necessary for this class of + elliptic equations this solver is intended to solve. + See ./CactusElliptic/EllBase/src/ for the classes of elliptic eq. @enddesc @calls @calledby @@ -117,30 +123,32 @@ void sor_confmetric(cGH *GH, @@*/ -void sor_flat(cGH *GH, +int sor_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { + int retval = ELL_NOSOLVER; switch (CCTK_GroupDimFromVarI(FieldIndex)) { case 1: - CCTK_WARN(1,"No 1D SOR solver implemented"); + CCTK_WARN(1,"sor_flat: No 1D SOR solver implemented"); break; case 2: - CCTK_WARN(1,"No 2D SOR solver implemented"); + CCTK_WARN(1,"sor_flat: No 2D SOR solver implemented"); break; case 3: - printf(" wrapper: %d %d %d \n",FieldIndex, MIndex, NIndex); - sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol); + retval = sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol); break; default: - CCTK_WARN(1,"Requested SOR solver for dimension equal" + CCTK_WARN(1,"sor_flat: Requested SOR solver for dimension equal" "zero or gt. three not implemented!"); break; } + return retval; + } |