aboutsummaryrefslogtreecommitdiff
path: root/src/sor_flat.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/sor_flat.c')
-rw-r--r--src/sor_flat.c321
1 files changed, 0 insertions, 321 deletions
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 <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-
-#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;chebit<sorit;chebit++)
- {
- omega = 1.0/(1.0 - 0.25*rjacobian*rjacobian*omega);
- }
- }
- if (accel_const)
- {
- omega = 1.8;
- }
-
- resnorm = 0.0;
-
- is = 1;
- js = 1;
- ks = 1;
-
- /*is = ((GH->cctk_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;k<ke;k=k+step)
- {
- for (j=js;j<je;j=j+step)
- {
- for (i=is;i<ie;i=i+step)
- {
-
- 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 = resnorm + fabs(residual);
-
- var[ijk] = var[ijk] - omega*residual/ac;
-
- }
- }
- }
-
- ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
- &resnorm, &glob_residual, CCTK_VARIABLE_REAL);
- if (ierr<0)
- {
- CCTK_WARN(1,"Reduction of Norm failed");
- }
-
-#ifdef DEBUG_ELLIPTIC
- printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
-#endif
-
- 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
- /* 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_residual<tol)
- {
- break;
- }
-
- }
-
- 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)
- {
- CCTK_WARN(1,"sor_flat: SOR Solver did not converge");
- retval = ELL_NOCONVERGENCE;
- }
-
- if (ans) free(ans);
- return retval;
-}
-
-
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 02:36:22 +0000