diff options
Diffstat (limited to 'src/Flat.c')
| -rw-r--r-- | src/Flat.c | 94 |
1 files changed, 46 insertions, 48 deletions
@@ -27,10 +27,10 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllSOR_Flat_c) int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol) + CCTK_REAL *AbsTol, CCTK_REAL *RelTol) { DECLARE_CCTK_PARAMETERS @@ -104,18 +104,18 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, sw[2]=1; if (Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf") == - ELLGET_NOTSET) + ELLGET_NOTSET) { - CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::inf not set, " - "setting to 1"); - finf = 1; + CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::inf not set, " + "setting to 1"); + finf = 1; } if (Ell_GetIntKey(&npow, "EllLinFlat::Bnd::Robin::falloff") - == ELLGET_NOTSET) + == ELLGET_NOTSET) { - CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::falloff not set, " - "setting to 1"); - npow = 1; + CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::falloff not set, " + "setting to 1"); + npow = 1; } } } @@ -124,7 +124,7 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET) { CCTK_WARN(1,"SORFlat3D: Ell::SORmaxit not set. " - "Maximum allowed iterations being set to 100."); + "Maximum allowed iterations being set to 100."); maxit = 100; } @@ -136,7 +136,7 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, { const char tmpstr[6] = "const"; CCTK_WARN(3, "SORFlat3D: Ell::SORaccel not set. " - "Omega being set to a constant value of 1.8."); + "Omega being set to a constant value of 1.8."); sor_accel = strdup(tmpstr); } @@ -149,15 +149,15 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, { if (CCTK_Equals(sor_accel, "cheb")) { - CCTK_INFO("SOR with Chebyshev acceleration"); + CCTK_INFO("SOR with Chebyshev acceleration"); } else if (CCTK_Equals(sor_accel, "const")) { - CCTK_INFO("SOR with hardcoded omega = 1.8"); + CCTK_INFO("SOR with hardcoded omega = 1.8"); } else if (CCTK_Equals(sor_accel, "none")) { - CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)"); + CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)"); } else { @@ -249,47 +249,47 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, { for (j=1; j<GH->cctk_lsh[1]-1; j++) { - for (i=1; i<GH->cctk_lsh[0]-1; i++) + for (i=1; i<GH->cctk_lsh[0]-1; i++) { if ((origin_sign + i + j + k)%2 == sorit%2) { - ac = ac_orig; - - 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); - - if (Mstorage) - { - ac += Mlin[ijk]; - } - - residual = ac * var[ijk] - + ae*var[ipjk] + aw*var[imjk] - + an*var[ijpk] + as*var[ijmk] - + at*var[ijkp] + ab*var[ijkm]; - - if (Nstorage) - { - residual += Nlin[ijk]; - } - - resnorm += fabs(residual); - - var[ijk] -= omega*residual/ac; + ac = ac_orig; + + 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); + + if (Mstorage) + { + ac += Mlin[ijk]; + } + + residual = ac * var[ijk] + + ae*var[ipjk] + aw*var[imjk] + + an*var[ijpk] + as*var[ijmk] + + at*var[ijkp] + ab*var[ijkm]; + + if (Nstorage) + { + residual += Nlin[ijk]; + } + + resnorm += fabs(residual); + + var[ijk] -= omega*residual/ac; } - } + } } } /* reduction operation on processor-local residual values */ if (CCTK_ReduceLocScalar(GH, -1, sum_handle, - &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) + &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) { CCTK_WARN(1,"SORFlat3D: Reduction of Norm failed"); } @@ -370,5 +370,3 @@ int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, return retval; } - - |