aboutsummaryrefslogtreecommitdiff
path: root/src/sor_confmetric.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/sor_confmetric.c')
-rw-r--r--src/sor_confmetric.c474
1 files changed, 0 insertions, 474 deletions
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 <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_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;ijk<nxyz;ijk++)
- {
-
- /* determinant */
- det = -(SQR(gxz[ijk])*gyy[ijk]) +
- 2*gxy[ijk]*gxz[ijk]*gyz[ijk] -
- gxx[ijk]*SQR(gyz[ijk]) -
- SQR(gxy[ijk])*gzz[ijk] +
- gxx[ijk]*gyy[ijk]*gzz[ijk];
-
- if (conformal)
- {
- pm4 = 1.0/pow(psi[ijk],4.0);
- p12 = pow(psi[ijk],12.0);
- }
- else
- {
- pm4 = 1.0;
- p12 = 1.0;
- }
-
-
- /* try to avoid constly division */
- detrec_pm4 = 1.0/det*pm4;
- sqrtdet = sqrt(det);
-
- uxx[ijk]=(-SQR(gyz[ijk]) + gyy[ijk]*gzz[ijk])*detrec_pm4;
- uxy[ijk]=( gxz[ijk]*gyz[ijk] - gxy[ijk]*gzz[ijk])*detrec_pm4;
- uxz[ijk]=(-gxz[ijk]*gyy[ijk] + gxy[ijk]*gyz[ijk])*detrec_pm4;
- uyy[ijk]=(-SQR(gxz[ijk]) + gxx[ijk]*gzz[ijk])*detrec_pm4;
- uyz[ijk]=( gxy[ijk]*gxz[ijk] - gxx[ijk]*gyz[ijk])*detrec_pm4;
- uzz[ijk]=(-SQR(gxy[ijk]) + gxx[ijk]*gyy[ijk])*detrec_pm4;
-
- det = det*p12;
- sqrtdet= sqrt(det);
-
- /* Rescaling for the uppermetric solver */
- if (Mstorage) Mlin[ijk] = Mlin[ijk]*sqrt(det);
- if (Nstorage) Nlin[ijk] = Nlin[ijk]*sqrt(det);
-
- uxx[ijk]=uxx[ijk]*dxdx*sqrtdet;
- uyy[ijk]=uyy[ijk]*dydy*sqrtdet;
- uzz[ijk]=uzz[ijk]*dzdz*sqrtdet;
- uxy[ijk]=uxy[ijk]*dxdy*sqrtdet;
- uxz[ijk]=uxz[ijk]*dxdz*sqrtdet;
- uyz[ijk]=uyz[ijk]*dydz*sqrtdet;
-
- }
-
- /*$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++)
- {
- 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;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];
- }
-
- resnorm = resnorm + fabs(residual);
- var[ijk] = var[ijk] - omega*residual/ac;
- }
- }
- }
-
- /* reduction operation on processor-local residual values */
- 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);
- if (uzz) free(uzz);
- if (uxy) free(uxy);
- if (uxz) free(uxz);
- if (uyz) free(uyz);
-
- return retval;
-}
-
-
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 02:04:36 +0000