From c8ca30566098863b88a3b785d2221273f2109688 Mon Sep 17 00:00:00 2001 From: allen Date: Thu, 3 Jan 2002 21:00:00 +0000 Subject: New tidied version of SOR, Linear solver uses grid sweeping which gives consistent results for multiprocessors git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@73 fa3da13c-9f13-4301-a575-cf5b8c5e1907 --- src/ConfMetric.c | 546 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/Flat.c | 373 +++++++++++++++++++++++++++++++++++ src/Startup.c | 46 +++-- src/Wrapper.c | 159 +++++++++++++++ src/make.code.defn | 2 +- src/sor_confmetric.c | 474 -------------------------------------------- src/sor_flat.c | 321 ------------------------------ src/sor_wrapper.c | 158 --------------- 8 files changed, 1110 insertions(+), 969 deletions(-) create mode 100644 src/ConfMetric.c create mode 100644 src/Flat.c create mode 100644 src/Wrapper.c delete mode 100644 src/sor_confmetric.c delete mode 100644 src/sor_flat.c delete mode 100644 src/sor_wrapper.c diff --git a/src/ConfMetric.c b/src/ConfMetric.c new file mode 100644 index 0000000..be9b0d2 --- /dev/null +++ b/src/ConfMetric.c @@ -0,0 +1,546 @@ + /*@@ + @file ConfMetric.c + @date Tue Sep 26 11:29:18 2000 + @author Gerd Lanfermann + @desc + Provides sor solver engine routines + @enddesc + @@*/ + +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "Boundary.h" +#include "Symmetry.h" +#include "Ell_DBstructure.h" +#include "EllBase.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_ConfMetric_c) + +#define SQR(a) ((a)*(a)) + +static int firstcall = 1; + +int SORConfMetric3D(cGH *GH, int *MetricPsiI, int conformal, + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); + + /*@@ + @routine SORConfMetric3D + @date Tue Sep 26 11:28:08 2000 + @author Joan Masso, Paul Walker, Gerd Lanfermann + @desc + This is a standalone sor solver engine, + called by the wrapper functions + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int SORConfMetric3D(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; + int origin_sign; + + /* The pointer to the metric fields */ + CCTK_REAL *gxx =NULL, *gxy =NULL; + CCTK_REAL *gxz =NULL, *gyy =NULL; + CCTK_REAL *gyz =NULL, *gzz =NULL; + CCTK_REAL *psi =NULL; + CCTK_REAL *Mlin=NULL, *Nlin=NULL; + CCTK_REAL *var =NULL; + + /* The inverse metric, allocated here */ + CCTK_REAL *uxx=NULL, *uyy=NULL, + *uzz=NULL, *uxz=NULL, + *uxy=NULL, *uyz=NULL; + + /* shortcuts for metric, psi, deltas, etc. */ + CCTK_REAL pm4, p12, det; + + CCTK_REAL dx,dy,dz; + CCTK_REAL dxdx, dydy, dzdz, + dxdy, dxdz, dydz; + + /* Some physical variables */ + CCTK_REAL omega, resnorm=0.0, residual=0.0; + CCTK_REAL glob_residual=0.0, rjacobian=0.0; + CCTK_REAL finf; + CCTK_INT npow; + CCTK_REAL tol; + + /* Iteration / stepping variables */ + int sorit, maxit; + int i,j,k; + int nxyz; + + /* 19 point stencil index */ + int ijk; + int ipjk, ijpk, ijkp, imjk, ijmk, ijkm; + int ipjpk, ipjmk, imjpk, imjmk; + int ipjkp, ipjkm, imjkp, imjkm; + int ijpkp, ijpkm, ijmkp, ijmkm; + + /* 19 point stencil coefficients */ + CCTK_REAL ac; + CCTK_REAL ae,aw,an,as,at,ab; + CCTK_REAL ane, anw, ase, asw, ate, atw, abe, abw; + CCTK_REAL atn, ats, abn, absol; + + /* Miscellaneous */ + int sum_handle=-1; + int sw[3], ierr; + int Mstorage=0, Nstorage=0; + size_t varsize; + CCTK_REAL detrec_pm4, sqrtdet; + char *robin=NULL; + char *sor_accel=NULL; + + + /* Get the reduction handle */ + sum_handle = CCTK_ReductionArrayHandle("sum"); + if (sum_handle<0) + { + CCTK_WARN(1,"SORConfMetric3D: " + "Cannot get handle for sum reduction"); + retval = ELL_FAILURE; + } + + /* 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 (Ell_GetStrKey(&robin, "EllLinConfMetric::Bnd::Robin")< 0) + { + CCTK_WARN(1,"SORConfMetric3D: Robin keys not set"); + } + + if (CCTK_EQUALS(robin,"yes")) + { + 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"); + } + + /* Get the maximum number of iterations allowed. */ + if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET) + { + CCTK_WARN(1,"SORConfMetric3D: Ell::SORmaxit not set. " + "Maximum allowed iterations being set to 100."); + maxit = 100; + } + + /* Only supports absolute tolerance */ + tol = AbsTol[0]; + + /* Get the type of acceleration to use. */ + if (Ell_GetStrKey(&sor_accel, "Ell::SORaccel") == ELLGET_NOTSET) + { + const char tmpstr[6] = "const"; + CCTK_WARN(1,"SORConfMetric3D: Ell::SORaccel not set. " + "Omega being set to a constant value of 1.8."); + sor_accel = strdup(tmpstr); + } + + /* Things to do only once! */ + if (firstcall==1) + { + if (CCTK_Equals(elliptic_verbose, "yes")) + { + if (CCTK_Equals(sor_accel, "cheb")) + { + CCTK_INFO("SOR with Chebyshev acceleration"); + } + else if (CCTK_Equals(sor_accel, "const")) + { + CCTK_INFO("SOR with hardcoded omega = 1.8"); + } + else if (CCTK_Equals(sor_accel, "none")) + { + CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)"); + } + else + { + CCTK_WARN(0, "SORConfMetric3D: sor_accel not set"); + } + } + firstcall = 0; + } + + /* + 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 = 0; + var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); + + 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]); + 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) + { + timelevel = 0; + 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) + { + timelevel = 0; + Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex); + if (Mlin) + { + Mstorage = 1; + } + else + { + CCTK_WARN(0, "SORConfMetric3D: Unable to get pointer to M."); + } + } + if (NIndex>=0) + { + timelevel = 0; + Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex); + if (Nlin) + { + Nstorage = 1; + } + else + { + CCTK_WARN(0, "SORConfMetric3D: Unable to get pointer to N."); + } + } + + /* 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]; + + /* Allocate the inverse metric */ + varsize = (size_t)CCTK_VarTypeSize(CCTK_VarTypeI(FieldIndex))*nxyz; + uxx = (CCTK_REAL*) malloc(varsize); + uxy = (CCTK_REAL*) malloc(varsize); + uxz = (CCTK_REAL*) malloc(varsize); + uyy = (CCTK_REAL*) malloc(varsize); + uyz = (CCTK_REAL*) malloc(varsize); + uzz = (CCTK_REAL*) malloc(varsize); + + if (!uxx || !uxy || !uxz || !uyy || !uyz || !uzz) + { + CCTK_WARN(0,"SORConfMetric3D: Memory allocation failed for upper metric"); + } + + /* calculate the differential coefficient */ + dxdx = 0.5/(dx*dx); + dydy = 0.5/(dy*dy); + dzdz = 0.5/(dz*dz); + dxdy = 0.25/(dx*dy); + dxdz = 0.25/(dx*dz); + dydz = 0.25/(dy*dz); + + /* Calculate the inverse metric */ + for (ijk=0; ijkcctk_lbnd[0] + GH->cctk_lbnd[1] + GH->cctk_lbnd[2])%2; + + /* start at 1 for historic (Fortran) reasons */ + for (sorit=1; sorit<=maxit; sorit++) + { + resnorm = 0.0; + + for (k=1; kcctk_lsh[2]-1; k++) + { + for (j=1; jcctk_lsh[1]-1; j++) + { + for (i=1; icctk_lsh[0]-1; i++) + { + if ((origin_sign+i+j+k)%2 == sorit%2) + { + 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; + + } + } + } + } + + /* reduction operation on processor-local residual values */ + if (CCTK_ReduceLocScalar(GH, -1, sum_handle, + &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) + { + CCTK_WARN(1,"SORConfMetric3D: 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,"SORConfMetric3D: CartSymVI failed in EllSOR loop"); + } + + /* apply Robin boundary conditions within loop */ + if (CCTK_EQUALS(robin,"yes")) + { + if (BndRobinVI(GH, sw, finf, npow, FieldIndex)<0) + { + CCTK_WARN(1, "SORConfMetric3D: BndRobinVI failed in EllSOR loop"); + break; + } + } + + /* synchronization of grid variable */ + ierr = CCTK_SyncGroupWithVarI(GH, FieldIndex); + if (ierr < 0) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ConfMetric3D: Synchronization of %s failed (Error: %d)", + CCTK_VarName(FieldIndex),ierr); + } + + /* Leave iteration loop if tolerance criterium is met */ + if (glob_residualtol) + { + CCTK_WARN(1, "SORConfMetric3D: SOR Solver did not converge"); + retval = ELL_NOCONVERGENCE; + } + + if (uxx) + { + free(uxx); + } + if (uyy) + { + free(uyy); + } + if (uzz) + { + free(uzz); + } + if (uxy) + { + free(uxy); + } + if (uxz) + { + free(uxz); + } + if (uyz) + { + free(uyz); + } + if (robin) + { + free(robin); + } + if (sor_accel) + { + free(sor_accel); + } + + return retval; +} diff --git a/src/Flat.c b/src/Flat.c new file mode 100644 index 0000000..6fcbb0e --- /dev/null +++ b/src/Flat.c @@ -0,0 +1,373 @@ + /*@@ + @file 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 +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "Boundary.h" +#include "Symmetry.h" +#include "Ell_DBstructure.h" +#include "EllBase.h" + +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); + +int SORFlat3D(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; + CCTK_REAL *Nlin=NULL; + CCTK_REAL *var =NULL; + + /* shortcuts for deltas,etc. */ + CCTK_REAL dx,dy,dz; + + /* Some physical variables */ + CCTK_REAL omega, resnorm, residual, glob_residual, rjacobian; + CCTK_REAL finf; + CCTK_INT npow; + CCTK_REAL tol; + + /* Iteration and stepping variables */ + int sorit, maxit; + int i, j, k; + int origin_sign; + + /* stencil index */ + int ijk; + int ipjk, ijpk, ijkp, imjk, ijmk, ijkm; + + /* Coefficients for the stencil... */ + CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab; + + /* Miscellaneous */ + int sum_handle; + int sw[3]; + int Mstorage=0, Nstorage=0; + static int firstcall = 1; + CCTK_REAL dx2rec, dy2rec, dz2rec; + char *robin = NULL; + char *sor_accel = NULL; + + /* Boundary conditions */ + int use_robin = 0; + + /* Avoid compiler warnings */ + RelTol = RelTol; + + /* Get the reduction handle */ + sum_handle = CCTK_ReductionArrayHandle("sum"); + if (sum_handle<0) + { + CCTK_WARN(1,"SORFlat3D: 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 (Ell_GetStrKey(&robin,"EllLinFlat::Bnd::Robin") < 0) + { + CCTK_WARN(2,"SORFlat3D: Robin keys not set"); + } + else + { + if (CCTK_Equals(robin,"yes")) + { + use_robin = 1; + + sw[0]=1; + sw[1]=1; + sw[2]=1; + + if (Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf") == + ELLGET_NOTSET) + { + 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) + { + CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::falloff not set, " + "setting to 1"); + npow = 1; + } + } + } + + /* Get the maximum number of iterations allowed. */ + if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET) + { + CCTK_WARN(1,"SORFlat3D: Ell::SORmaxit not set. " + "Maximum allowed iterations being set to 100."); + maxit = 100; + } + + /* Only supports absolute tolerance */ + tol = AbsTol[0]; + + /* Get the type of acceleration to use. */ + if (Ell_GetStrKey(&sor_accel, "Ell::SORaccel") == ELLGET_NOTSET) + { + const char tmpstr[6] = "const"; + CCTK_WARN(1, "SORFlat3D: Ell::SORaccel not set. " + "Omega being set to a constant value of 1.8."); + sor_accel = strdup(tmpstr); + } + + /* Things to do only once! */ + /* TODO: Need to handle this differently, since it may be called + for different reasons by the same code. */ + if (firstcall==1) + { + if (CCTK_Equals(elliptic_verbose, "yes")) + { + if (CCTK_Equals(sor_accel, "cheb")) + { + CCTK_INFO("SOR with Chebyshev acceleration"); + } + else if (CCTK_Equals(sor_accel, "const")) + { + CCTK_INFO("SOR with hardcoded omega = 1.8"); + } + else if (CCTK_Equals(sor_accel, "none")) + { + CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)"); + } + else + { + CCTK_WARN(0, "SORFlat3D: sor_accel not set"); + } + } + 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 (var == NULL) + { + CCTK_WARN(0, "SORFlat3D: Unable to get pointer to GF variable"); + } + + if (MIndex>=0) + { + Mlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, MIndex); + if (Mlin) + { + Mstorage = 1; + } + else + { + CCTK_WARN(0, "SORFlat3D: Unable to get pointer to M"); + } + } + + if (NIndex>=0) + { + Nlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, NIndex); + if (Nlin) + { + Nstorage = 1; + } + else + { + CCTK_WARN(0, "SORFlat3D: Unable to get pointer to N"); + } + } + + /* Shortcuts */ + dx = GH->cctk_delta_space[0]; + dy = GH->cctk_delta_space[1]; + dz = GH->cctk_delta_space[2]; + + dx2rec = 1.0/(dx*dx); + dy2rec = 1.0/(dy*dy); + dz2rec = 1.0/(dz*dz); + + ae = dx2rec; + aw = dx2rec; + an = dy2rec; + as = dy2rec; + at = dz2rec; + ab = dz2rec; + + ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec; + + /* Initialize omega. */ + /* TODO: Make it so that the value of the constant omega can be set. */ + if (CCTK_Equals(sor_accel, "const")) + { + omega = 1.8; + } + else + { + omega = 1.; + } + + /* Set the spectral radius of the Jacobi iteration. */ + /* TODO: I think this can be easily computed for flat metrics? */ + rjacobian = 0.999; + + /* Compute whether the origin of this processor's sub grid is "red" or + "black". */ + origin_sign = (GH->cctk_lbnd[0] + GH->cctk_lbnd[1] + GH->cctk_lbnd[2])%2; + + /* start at 1 for historic (Fortran) reasons */ + for (sorit=1; sorit<=maxit; sorit++) + { + resnorm = 0.; + + for (k=1; kcctk_lsh[2]-1; k++) + { + for (j=1; jcctk_lsh[1]-1; j++) + { + for (i=1; icctk_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; + } + } + } + } + + /* reduction operation on processor-local residual values */ + + if (CCTK_ReduceLocScalar(GH, -1, sum_handle, + &resnorm, &glob_residual, CCTK_VARIABLE_REAL)<0) + { + CCTK_WARN(1,"SORFlat3D: Reduction of Norm failed"); + } + +#ifdef DEBUG_ELLIPTIC + printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol); +#endif + + 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) + { + CCTK_WARN(1,"SORFlat3D: CartSymVI failed in EllSOR loop"); + break; + } + + /* apply Robin boundary conditions within loop */ + if (use_robin) + { + if (BndRobinVI(GH, sw, finf, npow, FieldIndex)<0) + { + CCTK_WARN(1,"SORFlat3D: BndRobinVI failed in EllSOR loop"); + break; + } + } + + /* synchronization of grid variable */ + CCTK_SyncGroupWithVarI(GH, FieldIndex); + + /* Leave iteration loop if tolerance criterium is met */ + if (glob_residualtol) + { + CCTK_WARN(1,"SORFlat3D: SOR Solver did not converge"); + retval = ELL_NOCONVERGENCE; + } + + if (use_robin) + { + free(robin); + } + if (sor_accel) + { + free(sor_accel); + } + + return retval; +} + + diff --git a/src/Startup.c b/src/Startup.c index dfe7d5f..19a75bd 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -15,21 +15,21 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusElliptic_EllSOR_Startup_c) void EllSOR_Register(CCTK_ARGUMENTS); -void sor_confmetric(cGH *GH, +void SORConfMetric(cGH *GH, int *MetricPsiI, int FieldI, int MI, int NI, CCTK_REAL *AbsTol, CCTK_REAL *RelTol); -void sor_flat(cGH *GH, - int FieldI, - int MI, - int NI, - CCTK_REAL *AbsTol, - CCTK_REAL *RelTol); - -/* routine registers the SOR solver "sor_confmetric" under the name "sor" +void SORFlat(cGH *GH, + int FieldI, + int MI, + int NI, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol); + +/* routine registers the SOR solver "SORConfMetric" under the name "sor" with the Elliptic Class "LinConfMetric". */ @@ -39,20 +39,20 @@ void EllSOR_Register(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - int err; /* Register the solver with the elliptic classes */ - err = Ell_RegisterSolver(sor_confmetric,"sor","Ell_LinConfMetric"); + err = Ell_RegisterSolver(SORConfMetric,"sor","Ell_LinConfMetric"); if (!err==ELL_SUCCESS) { - CCTK_WARN(0,"Failed to register sor for Ell_LinConfMetric"); + CCTK_WARN(0,"EllSOR_Register: Failed to register sor for Ell_LinConfMetric"); } - err = Ell_RegisterSolver(sor_flat,"sor","Ell_LinFlat"); + + err = Ell_RegisterSolver(SORFlat,"sor","Ell_LinFlat"); if (!err==ELL_SUCCESS) { - CCTK_WARN(0,"Failed to register sor for Ell_LinFlat"); + CCTK_WARN(0,"EllSOR_Register: Failed to register sor for Ell_LinFlat"); } /* These "keys" have to be same in other elliptic solver and in @@ -62,6 +62,22 @@ void EllSOR_Register(CCTK_ARGUMENTS) err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Robin"); err = Ell_CreateKey(CCTK_VARIABLE_STRING,"EllLinFlat::Bnd::Const"); - + /* Create a key to be associated with the maximum number of iterations + allowed. */ + err = Ell_CreateKey(CCTK_VARIABLE_INT, "Ell::SORmaxit"); + if (err != ELL_SUCCESS) + { + CCTK_WARN(0, "EllSOR_Register: " + "Failed to create integer key Ell::SORmaxit"); + } + + /* Create a key to be associated with the type of acceleration to be + used. */ + err = Ell_CreateKey(CCTK_VARIABLE_STRING, "Ell::SORaccel"); + if (err != ELL_SUCCESS) + { + CCTK_WARN(0, "EllSOR_Register: " + "Failed to create integer key Ell::SORaccel"); + } } diff --git a/src/Wrapper.c b/src/Wrapper.c new file mode 100644 index 0000000..b3ca1c5 --- /dev/null +++ b/src/Wrapper.c @@ -0,0 +1,159 @@ + /*@@ + @file Wrapper.c + @date Tue Aug 24 12:50:07 1999 + @author Gerd Lanfermann + @desc + The C wrapper, which calles the core Fortran routine, which + performs the actual solve. + We cannot derive the pointers to the GF data from the indices in + Fortran. So we do this here in C and then pass the everything + over to the Fortran routine. + + This wrapper is registers with the Elliptic solver registry + (not the Fortran file) , as coded up in ./CactusElliptic/EllBase + @enddesc + @@*/ + +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_FortranString.h" + +#include "EllBase.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusElliptic_EllSOR_Wrapper_c) + +int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); +int SORConfMetric3D(cGH *GH, int *MetricPsiI, int conformal, + int FieldIndex, int MIndex, int NIndex, + CCTK_REAL *AbsTol, CCTK_REAL *RelTol); +int SORConfMetric(cGH *GH, + int *MetricPsiI, + int FieldIndex, + int MIndex, + int NIndex, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol); +int SORFlat(cGH *GH, + int FieldIndex, + int MIndex, + int NIndex, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol); + +/*@@ + @routine SORConfMetric + @date Tue Sep 26 11:31:42 2000 + @author Gerd Lanfermann + @desc + elliptic solver wrapper which provides a interface to the + different n-dimensional SOR solvers for the conformal metric, + of which only 3D is coded currently. + + 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. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int SORConfMetric(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,"SORConfMetric: No 1D SOR solver implemented"); + break; + case 2: + CCTK_WARN(0,"SORConfMetric: No 2D SOR solver implemented"); + break; + case 3: + retval = SORConfMetric3D(GH, MetricPsiI, 1, + FieldIndex, MIndex, NIndex, + AbsTol, RelTol); + break; + default: + CCTK_WARN(1,"SORConfMetric: Solver only implemented for 3D"); + break; + } + + return retval; +} + +/*@@ + @routine SORFlat + @date Tue Sep 26 11:31:42 2000 + @author Gerd Lanfermann + @desc + Elliptic solver wrapper which provides a interface to the + different n-dimensional SOR solvers (flat case), + of which only 3D is coded currently. + + This wrapper is registered and if it is being called with + a n-dimensional grid function, it then picks the correct solver. + + 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 equations. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int SORFlat(cGH *GH, + int FieldIndex, + int MIndex, + int NIndex, + CCTK_REAL *AbsTol, + CCTK_REAL *RelTol) +{ + int retval; + + retval = ELL_NOSOLVER; + + switch (CCTK_GroupDimFromVarI(FieldIndex)) + { + case 1: + CCTK_WARN(1,"SORFlat: No 1D SOR solver implemented"); + break; + case 2: + CCTK_WARN(1,"SORFlat: No 2D SOR solver implemented"); + break; + case 3: + retval = SORFlat3D(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol); + break; + default: + CCTK_WARN(1,"SORFlat: Solver only implemented for 3D"); + break; + } + + return retval; + +} + diff --git a/src/make.code.defn b/src/make.code.defn index 4e86518..117e3e5 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = Startup.c sor_wrapper.c sor_confmetric.c sor_flat.c +SRCS = Startup.c Wrapper.c ConfMetric.c Flat.c # Subdirectories containing source files SUBDIRS = diff --git a/src/sor_confmetric.c b/src/sor_confmetric.c deleted file mode 100644 index 15c4a9d..0000000 --- a/src/sor_confmetric.c +++ /dev/null @@ -1,474 +0,0 @@ - /*@@ - @file sor_confmetric.c - @date Tue Sep 26 11:29:18 2000 - @author Gerd Lanfermann - @desc - Provides sor solver engine routines - @enddesc - @@*/ - - -#include -#include -#include - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -#include "Boundary.h" -#include "Symmetry.h" -#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, - int FieldIndex, int MIndex, int NIndex, - CCTK_REAL *AbsTol, CCTK_REAL *RelTol); - - /*@@ - @routine sor_confmetric - @date Tue Sep 26 11:28:08 2000 - @author Joan Masso, Paul Walker, Gerd Lanfermann - @desc - This is a standalone sor solver engine, - called by the wrapper functions - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - -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; - CCTK_REAL *gyz =NULL, *gzz =NULL; - CCTK_REAL *psi =NULL; - CCTK_REAL *Mlin=NULL, *Nlin=NULL; - CCTK_REAL *var =NULL; - - /* The inverse metric, allocated here */ - CCTK_REAL *uxx=NULL, *uyy=NULL, - *uzz=NULL, *uxz=NULL, - *uxy=NULL, *uyz=NULL; - - /* shortcuts for metric, psi, deltas, etc. */ - CCTK_REAL pm4, p12, det; - - CCTK_REAL dx,dy,dz; - CCTK_REAL dxdx, dydy, dzdz, - dxdy, dxdz, dydz; - - /* Some physical variables */ - int accel_cheb=0, accel_const=0; - int chebit; - CCTK_REAL omega, resnorm=0.0, residual=0.0; - CCTK_REAL glob_residual=0.0, rjacobian=0.0; - CCTK_REAL finf; - CCTK_INT npow; - CCTK_REAL tol; - - - /* Iteration / stepping variables */ - int sorit; - int i,is,ie; - int j,js,je; - int k,ks,ke,kstep; - int nxyz; - - /* 19 point stencil index */ - int ijk; - int ipjk, ijpk, ijkp, imjk, ijmk, ijkm; - int ipjpk, ipjmk, imjpk, imjmk; - int ipjkp, ipjkm, imjkp, imjkm; - int ijpkp, ijpkm, ijmkp, ijmkm; - - /* 19 point stencil coefficients */ - CCTK_REAL ac; - CCTK_REAL ae,aw,an,as,at,ab; - CCTK_REAL ane, anw, ase, asw, ate, atw, abe, abw; - CCTK_REAL atn, ats, abn, absol; - - /* Miscellaneous */ - int sum_handle=-1; - int sw[3], ierr; - int Mstorage=0, Nstorage=0; - size_t varsize; - static int firstcall = 1; - CCTK_REAL detrec_pm4, sqrtdet; - - - /* 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")) - { - 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 (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 - { - printf("SOR with unaccelearted relaxation (omega = 1)\n"); - } - } - firstcall = 0; - } - - - /* 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 = 0; - var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex); - - 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]); - 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) - { - timelevel = 0; - 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) - { - timelevel = 0; - Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex); - Mstorage = 1; - } - if (NIndex>=0) - { - timelevel = 0; - Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex); - Nstorage = 1; - } - - /* 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]; - - - /* Allocate the inverse metric */ - varsize = (size_t)CCTK_VarTypeSize(CCTK_VarTypeI(FieldIndex))*nxyz; - if (!(uxx = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - if (!(uyy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - if (!(uzz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - if (!(uxy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - if (!(uxz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - if (!(uyz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed"); - - /* PreCalc the differential coeff. */ - dxdx = 1.0/(2.0*dx*dx); - dydy = 1.0/(2.0*dy*dy); - dzdz = 1.0/(2.0*dz*dz); - dxdy = 1.0/(4.0*dx*dy); - dxdz = 1.0/(4.0*dx*dz); - dydz = 1.0/(4.0*dy*dz); - - /* Calculate the inverse metric */ - - for (ijk=0;ijkcctk_lsh[2]-1;k++) { - for (j=1;jcctk_lsh[1]-1;j++) - { - for (i=1;icctk_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]); - } - } - }$*/ - - - /* iteration boundaries (ks set inside loop)*/ - 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++) { - - omega = 1.0; - rjacobian = 1.0; - - if (accel_cheb) - { - for (chebit=2;chebitcctk_lsh[2]==3) - ks = 1; - - for (k=ks;kcctk_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_residualtol) - { - 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); - if (uzz) free(uzz); - if (uxy) free(uxy); - if (uxz) free(uxz); - if (uyz) free(uyz); - - return retval; -} - - diff --git a/src/sor_flat.c b/src/sor_flat.c deleted file mode 100644 index 1b30f1d..0000000 --- a/src/sor_flat.c +++ /dev/null @@ -1,321 +0,0 @@ - /*@@ - @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 -#include -#include - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -#include "Boundary.h" -#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); - -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; - CCTK_REAL *Nlin=NULL; - CCTK_REAL *var =NULL; - - /* shortcuts for deltas,etc. */ - CCTK_REAL dx,dy,dz; - - /* Some physical variables */ - int accel_cheb=0, accel_const=0; - int chebit; - CCTK_REAL omega, resnorm, residual, glob_residual, rjacobian; - 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,step; - int timelevel; - - /* stencil index */ - int ijk; - int ipjk, ijpk, ijkp, imjk, ijmk, ijkm; - - /* Coeeficients for the stencil... */ - CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab; - - /* Miscellaneous */ - int sum_handle=-1; - int sw[3], ierr; - 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) - { - 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) */ - - 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")) - { - accel_const = 1; - } - else if (CCTK_EQUALS(sor_accel,"cheb")) - { - accel_cheb = 1; - } - - /* Things to do only once! */ - 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 - { - 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 */ - /* get the current time level */ - - var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, FieldIndex); - - if (MIndex>=0) - { - Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex); - Mstorage = 1; - } - else - { - CCTK_WARN(0,"Stop"); - } - if (NIndex>=0) - { - 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]; - - dx2rec = 1.0/(dx*dx); - dy2rec = 1.0/(dy*dy); - dz2rec = 1.0/(dz*dz); - - ae = dx2rec; - aw = dx2rec; - an = dy2rec; - as = dy2rec; - at = dz2rec; - ab = dz2rec; - - ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec; - - - ie = GH->cctk_lsh[0]-1; - je = GH->cctk_lsh[1]-1; - ke = GH->cctk_lsh[2]-1; - - /* start at 1 for historic (Fortran) reasons */ - for (sorit=1; sorit<=maxit; sorit++) - { - omega = 1.0; - rjacobian = 1.0; - - if (accel_cheb) - { - for (chebit=2;chebitcctk_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]); - - -#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,"sor_flat: CartSymVI failed in EllSOR loop"); - break; - } - - BndScalarVI(GH,sw,zero,FieldIndex); - /* apply boundary conditions within loop */ - /* 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); - - /* Leave iteration loop if tolerance criterium is met */ - if (glob_residualtol) - { - CCTK_WARN(1,"sor_flat: SOR Solver did not converge"); - retval = ELL_NOCONVERGENCE; - } - - if (ans) free(ans); - return retval; -} - - diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c deleted file mode 100644 index 3df20c5..0000000 --- a/src/sor_wrapper.c +++ /dev/null @@ -1,158 +0,0 @@ - /*@@ - @file jacobi_wrapper.c - @date Tue Aug 24 12:50:07 1999 - @author Gerd Lanfermann - @desc - The C wrapper, which calles the core Fortran routine, which - performs the actual solve. - We cannot derive the pointers to the GF data from the indeces in - Fortran. So we do this here in C and then pass the everything - over to the Fortran routine. - - This wrapper is registers with the Elliptic solver registry - (not the Fortran file) , as coded up in ./CactusElliptic/EllBase - @enddesc - @@*/ - -#include -#include -#include - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_FortranString.h" - -#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, - 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, - CCTK_REAL *AbsTol, - CCTK_REAL *RelTol); - -/*@@ - @routine sor_confmetric - @date Tue Sep 26 11:31:42 2000 - @author Gerd Lanfermann - @desc - elliptic solver wrapper which provides a interface to the - different n-dimensional SOR solvers for the conformal metric, - of which only 3D is coded currently. - - 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. - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - -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,"sor_confmetric: No 1D SOR solver implemented"); - break; - case 2: - CCTK_WARN(0,"sor_confmetric: No 2D SOR solver implemented"); - break; - case 3: - retval = sor_confmetric_3d(GH, MetricPsiI, 1, - FieldIndex, MIndex, NIndex, - AbsTol, RelTol); - break; - default: - CCTK_WARN(1,"sor_confmetric: Requested SOR solver for dimension equal " - "zero or gt. three not implemented!"); - break; - } - - return retval; -} - -/*@@ - @routine sor_flat - @date Tue Sep 26 11:31:42 2000 - @author Gerd Lanfermann - @desc - elliptic solver wrapper which provides a interface to the - different n-dimensional SOR solvers (flat case), - of which only 3D is coded currently. - - 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 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 - @history - - @endhistory - -@@*/ - -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,"sor_flat: No 1D SOR solver implemented"); - break; - case 2: - CCTK_WARN(1,"sor_flat: No 2D SOR solver implemented"); - break; - case 3: - retval = sor_flat_3d(GH, FieldIndex, MIndex, NIndex, AbsTol, RelTol); - break; - default: - CCTK_WARN(1,"sor_flat: Requested SOR solver for dimension equal" - "zero or gt. three not implemented!"); - break; - } - - return retval; - -} - -- cgit v1.2.3