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.c310
1 files changed, 192 insertions, 118 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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 02:37:08 +0000