diff options
| author | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-02-13 15:27:22 -0600 |
|---|---|---|
| committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-02-13 15:27:22 -0600 |
| commit | 9eea7b4930f6f85e477139187eb93e4cef2bef3d (patch) | |
| tree | 5a1283291ceed40abbeece2695535b707814e0ce /ML_ADM | |
| parent | 76d87a0dd9b19a671d422fcefc17c7aabd223c08 (diff) | |
Update auto-generated code
Diffstat (limited to 'ML_ADM')
| -rw-r--r-- | ML_ADM/src/Differencing.h | 24 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_Minkowski.c | 122 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_RHS.c | 566 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_boundary.c | 122 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_constraints.c | 450 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_constraints_boundary.c | 94 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_convertFromADMBase.c | 186 | ||||
| -rw-r--r-- | ML_ADM/src/ML_ADM_convertToADMBase.c | 202 |
8 files changed, 876 insertions, 890 deletions
diff --git a/ML_ADM/src/Differencing.h b/ML_ADM/src/Differencing.h index ae119d2..aa21060 100644 --- a/ML_ADM/src/Differencing.h +++ b/ML_ADM/src/Differencing.h @@ -1,12 +1,12 @@ -#define PDstandardNth1(u,i,j,k) (p1o12dx*(-8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(k))] + 8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(k))])) -#define PDstandardNth2(u,i,j,k) (p1o12dy*(-8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(k))] + 8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(k))])) -#define PDstandardNth3(u,i,j,k) (p1o12dz*(-8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(-1 + k))] + 8*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(-2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(2 + k))])) -#define PDstandardNth11(u,i,j,k) (pm1o12dx2*(-16*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(k))] + 30*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(k))])) -#define PDstandardNth22(u,i,j,k) (pm1o12dy2*(-16*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(k))] + 30*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(k))])) -#define PDstandardNth33(u,i,j,k) (pm1o12dz2*(-16*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(1 + k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(2 + k))] + 30*(u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(j),(int)(k))])) -#define PDstandardNth12(u,i,j,k) (p1o144dxdy*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(-1 + j),(int)(k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(-1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(1 + j),(int)(k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(-1 + j),(int)(k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(-1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(1 + j),(int)(k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(-2 + j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(2 + j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(2 + j),(int)(k))])) -#define PDstandardNth13(u,i,j,k) (p1o144dxdz*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(-1 + k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(1 + k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(-1 + k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(1 + k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(-2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(2 + k))])) -#define PDstandardNth21(u,i,j,k) (p1o144dxdy*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(-1 + j),(int)(k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(-1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(1 + j),(int)(k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(-1 + j),(int)(k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(-1 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(1 + j),(int)(k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(-2 + j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(2 + j),(int)(k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(-2 + j),(int)(k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(2 + j),(int)(k))])) -#define PDstandardNth23(u,i,j,k) (p1o144dydz*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(-1 + k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(1 + k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(-1 + k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(1 + k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(-2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(2 + k))])) -#define PDstandardNth31(u,i,j,k) (p1o144dxdz*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(-1 + k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(1 + k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(-1 + k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(-1 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(1 + i),(int)(j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(1 + k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(-2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(-2 + i),(int)(j),(int)(2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(2 + i),(int)(j),(int)(2 + k))])) -#define PDstandardNth32(u,i,j,k) (p1o144dydz*(-64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(-1 + k))]) + 64*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(1 + k))]) + 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(-1 + k))]) - 8*((u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-1 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(1 + j),(int)(2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(-1 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(1 + k))]) + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(-2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(-2 + j),(int)(2 + k))] - (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(-2 + k))] + (u)[CCTK_GFINDEX3D(cctkGH,(int)(i),(int)(2 + j),(int)(2 + k))])) +#define PDstandardNth1(u,i,j,k) (p1o12dx*(-8*(u)[index+di*(-1)+dj*(0)+dk*(0)] + 8*(u)[index+di*(1)+dj*(0)+dk*(0)] + (u)[index+di*(-2)+dj*(0)+dk*(0)] - (u)[index+di*(2)+dj*(0)+dk*(0)])) +#define PDstandardNth2(u,i,j,k) (p1o12dy*(-8*(u)[index+di*(0)+dj*(-1)+dk*(0)] + 8*(u)[index+di*(0)+dj*(1)+dk*(0)] + (u)[index+di*(0)+dj*(-2)+dk*(0)] - (u)[index+di*(0)+dj*(2)+dk*(0)])) +#define PDstandardNth3(u,i,j,k) (p1o12dz*(-8*(u)[index+di*(0)+dj*(0)+dk*(-1)] + 8*(u)[index+di*(0)+dj*(0)+dk*(1)] + (u)[index+di*(0)+dj*(0)+dk*(-2)] - (u)[index+di*(0)+dj*(0)+dk*(2)])) +#define PDstandardNth11(u,i,j,k) (pm1o12dx2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(-1)+dj*(0)+dk*(0)] + (u)[index+di*(1)+dj*(0)+dk*(0)]) + (u)[index+di*(-2)+dj*(0)+dk*(0)] + (u)[index+di*(2)+dj*(0)+dk*(0)])) +#define PDstandardNth22(u,i,j,k) (pm1o12dy2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(0)+dj*(-1)+dk*(0)] + (u)[index+di*(0)+dj*(1)+dk*(0)]) + (u)[index+di*(0)+dj*(-2)+dk*(0)] + (u)[index+di*(0)+dj*(2)+dk*(0)])) +#define PDstandardNth33(u,i,j,k) (pm1o12dz2*(30*(u)[index+di*(0)+dj*(0)+dk*(0)] - 16*((u)[index+di*(0)+dj*(0)+dk*(-1)] + (u)[index+di*(0)+dj*(0)+dk*(1)]) + (u)[index+di*(0)+dj*(0)+dk*(-2)] + (u)[index+di*(0)+dj*(0)+dk*(2)])) +#define PDstandardNth12(u,i,j,k) (p1o144dxdy*(-64*((u)[index+di*(-1)+dj*(1)+dk*(0)] + (u)[index+di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[index+di*(-1)+dj*(-1)+dk*(0)] + (u)[index+di*(1)+dj*(1)+dk*(0)]) + 8*((u)[index+di*(-1)+dj*(2)+dk*(0)] + (u)[index+di*(1)+dj*(-2)+dk*(0)] + (u)[index+di*(-2)+dj*(1)+dk*(0)] + (u)[index+di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[index+di*(-1)+dj*(-2)+dk*(0)] + (u)[index+di*(1)+dj*(2)+dk*(0)] + (u)[index+di*(-2)+dj*(-1)+dk*(0)] + (u)[index+di*(2)+dj*(1)+dk*(0)]) + (u)[index+di*(-2)+dj*(-2)+dk*(0)] - (u)[index+di*(-2)+dj*(2)+dk*(0)] - (u)[index+di*(2)+dj*(-2)+dk*(0)] + (u)[index+di*(2)+dj*(2)+dk*(0)])) +#define PDstandardNth13(u,i,j,k) (p1o144dxdz*(-64*((u)[index+di*(-1)+dj*(0)+dk*(1)] + (u)[index+di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[index+di*(-1)+dj*(0)+dk*(-1)] + (u)[index+di*(1)+dj*(0)+dk*(1)]) + 8*((u)[index+di*(-1)+dj*(0)+dk*(2)] + (u)[index+di*(1)+dj*(0)+dk*(-2)] + (u)[index+di*(-2)+dj*(0)+dk*(1)] + (u)[index+di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[index+di*(-1)+dj*(0)+dk*(-2)] + (u)[index+di*(1)+dj*(0)+dk*(2)] + (u)[index+di*(-2)+dj*(0)+dk*(-1)] + (u)[index+di*(2)+dj*(0)+dk*(1)]) + (u)[index+di*(-2)+dj*(0)+dk*(-2)] - (u)[index+di*(-2)+dj*(0)+dk*(2)] - (u)[index+di*(2)+dj*(0)+dk*(-2)] + (u)[index+di*(2)+dj*(0)+dk*(2)])) +#define PDstandardNth21(u,i,j,k) (p1o144dxdy*(-64*((u)[index+di*(-1)+dj*(1)+dk*(0)] + (u)[index+di*(1)+dj*(-1)+dk*(0)]) + 64*((u)[index+di*(-1)+dj*(-1)+dk*(0)] + (u)[index+di*(1)+dj*(1)+dk*(0)]) + 8*((u)[index+di*(-1)+dj*(2)+dk*(0)] + (u)[index+di*(1)+dj*(-2)+dk*(0)] + (u)[index+di*(-2)+dj*(1)+dk*(0)] + (u)[index+di*(2)+dj*(-1)+dk*(0)]) - 8*((u)[index+di*(-1)+dj*(-2)+dk*(0)] + (u)[index+di*(1)+dj*(2)+dk*(0)] + (u)[index+di*(-2)+dj*(-1)+dk*(0)] + (u)[index+di*(2)+dj*(1)+dk*(0)]) + (u)[index+di*(-2)+dj*(-2)+dk*(0)] - (u)[index+di*(-2)+dj*(2)+dk*(0)] - (u)[index+di*(2)+dj*(-2)+dk*(0)] + (u)[index+di*(2)+dj*(2)+dk*(0)])) +#define PDstandardNth23(u,i,j,k) (p1o144dydz*(-64*((u)[index+di*(0)+dj*(-1)+dk*(1)] + (u)[index+di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[index+di*(0)+dj*(-1)+dk*(-1)] + (u)[index+di*(0)+dj*(1)+dk*(1)]) + 8*((u)[index+di*(0)+dj*(-1)+dk*(2)] + (u)[index+di*(0)+dj*(1)+dk*(-2)] + (u)[index+di*(0)+dj*(-2)+dk*(1)] + (u)[index+di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[index+di*(0)+dj*(-1)+dk*(-2)] + (u)[index+di*(0)+dj*(1)+dk*(2)] + (u)[index+di*(0)+dj*(-2)+dk*(-1)] + (u)[index+di*(0)+dj*(2)+dk*(1)]) + (u)[index+di*(0)+dj*(-2)+dk*(-2)] - (u)[index+di*(0)+dj*(-2)+dk*(2)] - (u)[index+di*(0)+dj*(2)+dk*(-2)] + (u)[index+di*(0)+dj*(2)+dk*(2)])) +#define PDstandardNth31(u,i,j,k) (p1o144dxdz*(-64*((u)[index+di*(-1)+dj*(0)+dk*(1)] + (u)[index+di*(1)+dj*(0)+dk*(-1)]) + 64*((u)[index+di*(-1)+dj*(0)+dk*(-1)] + (u)[index+di*(1)+dj*(0)+dk*(1)]) + 8*((u)[index+di*(-1)+dj*(0)+dk*(2)] + (u)[index+di*(1)+dj*(0)+dk*(-2)] + (u)[index+di*(-2)+dj*(0)+dk*(1)] + (u)[index+di*(2)+dj*(0)+dk*(-1)]) - 8*((u)[index+di*(-1)+dj*(0)+dk*(-2)] + (u)[index+di*(1)+dj*(0)+dk*(2)] + (u)[index+di*(-2)+dj*(0)+dk*(-1)] + (u)[index+di*(2)+dj*(0)+dk*(1)]) + (u)[index+di*(-2)+dj*(0)+dk*(-2)] - (u)[index+di*(-2)+dj*(0)+dk*(2)] - (u)[index+di*(2)+dj*(0)+dk*(-2)] + (u)[index+di*(2)+dj*(0)+dk*(2)])) +#define PDstandardNth32(u,i,j,k) (p1o144dydz*(-64*((u)[index+di*(0)+dj*(-1)+dk*(1)] + (u)[index+di*(0)+dj*(1)+dk*(-1)]) + 64*((u)[index+di*(0)+dj*(-1)+dk*(-1)] + (u)[index+di*(0)+dj*(1)+dk*(1)]) + 8*((u)[index+di*(0)+dj*(-1)+dk*(2)] + (u)[index+di*(0)+dj*(1)+dk*(-2)] + (u)[index+di*(0)+dj*(-2)+dk*(1)] + (u)[index+di*(0)+dj*(2)+dk*(-1)]) - 8*((u)[index+di*(0)+dj*(-1)+dk*(-2)] + (u)[index+di*(0)+dj*(1)+dk*(2)] + (u)[index+di*(0)+dj*(-2)+dk*(-1)] + (u)[index+di*(0)+dj*(2)+dk*(1)]) + (u)[index+di*(0)+dj*(-2)+dk*(-2)] - (u)[index+di*(0)+dj*(-2)+dk*(2)] - (u)[index+di*(0)+dj*(2)+dk*(-2)] + (u)[index+di*(0)+dj*(2)+dk*(2)])) diff --git a/ML_ADM/src/ML_ADM_Minkowski.c b/ML_ADM/src/ML_ADM_Minkowski.c index 4420395..ace3b38 100644 --- a/ML_ADM/src/ML_ADM_Minkowski.c +++ b/ML_ADM/src/ML_ADM_Minkowski.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_Minkowski_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_Minkowski_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_Minkowski_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_IN /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,18 +87,18 @@ void ML_ADM_Minkowski_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_IN i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ /* Declare local copies of grid functions */ - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; - CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + // CCTK_REAL alphaL = INITVALUE; + // CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ @@ -116,37 +114,37 @@ void ML_ADM_Minkowski_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_IN /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - g11L = 1; + CCTK_REAL const g11L = 1; - g12L = 0; + CCTK_REAL const g12L = 0; - g13L = 0; + CCTK_REAL const g13L = 0; - g22L = 1; + CCTK_REAL const g22L = 1; - g23L = 0; + CCTK_REAL const g23L = 0; - g33L = 1; + CCTK_REAL const g33L = 1; - K11L = 0; + CCTK_REAL const K11L = 0; - K12L = 0; + CCTK_REAL const K12L = 0; - K13L = 0; + CCTK_REAL const K13L = 0; - K22L = 0; + CCTK_REAL const K22L = 0; - K23L = 0; + CCTK_REAL const K23L = 0; - K33L = 0; + CCTK_REAL const K33L = 0; - alphaL = 1; + CCTK_REAL const alphaL = 1; - beta1L = 0; + CCTK_REAL const beta1L = 0; - beta2L = 0; + CCTK_REAL const beta2L = 0; - beta3L = 0; + CCTK_REAL const beta3L = 0; /* Copy local copies back to grid functions */ diff --git a/ML_ADM/src/ML_ADM_RHS.c b/ML_ADM/src/ML_ADM_RHS.c index 8c25677..4170199 100644 --- a/ML_ADM/src/ML_ADM_RHS.c +++ b/ML_ADM/src/ML_ADM_RHS.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_RHS_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT cons /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,404 +87,404 @@ void ML_ADM_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT cons i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ - CCTK_REAL detg = INITVALUE; - CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; - CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; - CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; - CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; - CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; - CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; - CCTK_REAL trK = INITVALUE; + // CCTK_REAL detg = INITVALUE; + // CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; + // CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; + // CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; + // CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; + // CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; + // CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; + // CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; + // CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; + // CCTK_REAL trK = INITVALUE; /* Declare local copies of grid functions */ - CCTK_REAL alphaL = INITVALUE, alpharhsL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g11rhsL = INITVALUE, g12L = INITVALUE, g12rhsL = INITVALUE, g13L = INITVALUE, g13rhsL = INITVALUE; - CCTK_REAL g22L = INITVALUE, g22rhsL = INITVALUE, g23L = INITVALUE, g23rhsL = INITVALUE, g33L = INITVALUE, g33rhsL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K11rhsL = INITVALUE, K12L = INITVALUE, K12rhsL = INITVALUE, K13L = INITVALUE, K13rhsL = INITVALUE; - CCTK_REAL K22L = INITVALUE, K22rhsL = INITVALUE, K23L = INITVALUE, K23rhsL = INITVALUE, K33L = INITVALUE, K33rhsL = INITVALUE; + // CCTK_REAL alphaL = INITVALUE, alpharhsL = INITVALUE; + // CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g11rhsL = INITVALUE, g12L = INITVALUE, g12rhsL = INITVALUE, g13L = INITVALUE, g13rhsL = INITVALUE; + // CCTK_REAL g22L = INITVALUE, g22rhsL = INITVALUE, g23L = INITVALUE, g23rhsL = INITVALUE, g33L = INITVALUE, g33rhsL = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K11rhsL = INITVALUE, K12L = INITVALUE, K12rhsL = INITVALUE, K13L = INITVALUE, K13rhsL = INITVALUE; + // CCTK_REAL K22L = INITVALUE, K22rhsL = INITVALUE, K23L = INITVALUE, K23rhsL = INITVALUE, K33L = INITVALUE, K33rhsL = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ - CCTK_REAL PDstandardNth1alpha = INITVALUE; - CCTK_REAL PDstandardNth2alpha = INITVALUE; - CCTK_REAL PDstandardNth3alpha = INITVALUE; - CCTK_REAL PDstandardNth11alpha = INITVALUE; - CCTK_REAL PDstandardNth22alpha = INITVALUE; - CCTK_REAL PDstandardNth33alpha = INITVALUE; - CCTK_REAL PDstandardNth12alpha = INITVALUE; - CCTK_REAL PDstandardNth13alpha = INITVALUE; - CCTK_REAL PDstandardNth23alpha = INITVALUE; - CCTK_REAL PDstandardNth1beta1 = INITVALUE; - CCTK_REAL PDstandardNth2beta1 = INITVALUE; - CCTK_REAL PDstandardNth3beta1 = INITVALUE; - CCTK_REAL PDstandardNth1beta2 = INITVALUE; - CCTK_REAL PDstandardNth2beta2 = INITVALUE; - CCTK_REAL PDstandardNth3beta2 = INITVALUE; - CCTK_REAL PDstandardNth1beta3 = INITVALUE; - CCTK_REAL PDstandardNth2beta3 = INITVALUE; - CCTK_REAL PDstandardNth3beta3 = INITVALUE; - CCTK_REAL PDstandardNth1g11 = INITVALUE; - CCTK_REAL PDstandardNth2g11 = INITVALUE; - CCTK_REAL PDstandardNth3g11 = INITVALUE; - CCTK_REAL PDstandardNth22g11 = INITVALUE; - CCTK_REAL PDstandardNth33g11 = INITVALUE; - CCTK_REAL PDstandardNth23g11 = INITVALUE; - CCTK_REAL PDstandardNth1g12 = INITVALUE; - CCTK_REAL PDstandardNth2g12 = INITVALUE; - CCTK_REAL PDstandardNth3g12 = INITVALUE; - CCTK_REAL PDstandardNth33g12 = INITVALUE; - CCTK_REAL PDstandardNth12g12 = INITVALUE; - CCTK_REAL PDstandardNth13g12 = INITVALUE; - CCTK_REAL PDstandardNth23g12 = INITVALUE; - CCTK_REAL PDstandardNth1g13 = INITVALUE; - CCTK_REAL PDstandardNth2g13 = INITVALUE; - CCTK_REAL PDstandardNth3g13 = INITVALUE; - CCTK_REAL PDstandardNth22g13 = INITVALUE; - CCTK_REAL PDstandardNth12g13 = INITVALUE; - CCTK_REAL PDstandardNth13g13 = INITVALUE; - CCTK_REAL PDstandardNth23g13 = INITVALUE; - CCTK_REAL PDstandardNth1g22 = INITVALUE; - CCTK_REAL PDstandardNth2g22 = INITVALUE; - CCTK_REAL PDstandardNth3g22 = INITVALUE; - CCTK_REAL PDstandardNth11g22 = INITVALUE; - CCTK_REAL PDstandardNth33g22 = INITVALUE; - CCTK_REAL PDstandardNth13g22 = INITVALUE; - CCTK_REAL PDstandardNth1g23 = INITVALUE; - CCTK_REAL PDstandardNth2g23 = INITVALUE; - CCTK_REAL PDstandardNth3g23 = INITVALUE; - CCTK_REAL PDstandardNth11g23 = INITVALUE; - CCTK_REAL PDstandardNth12g23 = INITVALUE; - CCTK_REAL PDstandardNth13g23 = INITVALUE; - CCTK_REAL PDstandardNth23g23 = INITVALUE; - CCTK_REAL PDstandardNth1g33 = INITVALUE; - CCTK_REAL PDstandardNth2g33 = INITVALUE; - CCTK_REAL PDstandardNth3g33 = INITVALUE; - CCTK_REAL PDstandardNth11g33 = INITVALUE; - CCTK_REAL PDstandardNth22g33 = INITVALUE; - CCTK_REAL PDstandardNth12g33 = INITVALUE; - CCTK_REAL PDstandardNth1K11 = INITVALUE; - CCTK_REAL PDstandardNth2K11 = INITVALUE; - CCTK_REAL PDstandardNth3K11 = INITVALUE; - CCTK_REAL PDstandardNth1K12 = INITVALUE; - CCTK_REAL PDstandardNth2K12 = INITVALUE; - CCTK_REAL PDstandardNth3K12 = INITVALUE; - CCTK_REAL PDstandardNth1K13 = INITVALUE; - CCTK_REAL PDstandardNth2K13 = INITVALUE; - CCTK_REAL PDstandardNth3K13 = INITVALUE; - CCTK_REAL PDstandardNth1K22 = INITVALUE; - CCTK_REAL PDstandardNth2K22 = INITVALUE; - CCTK_REAL PDstandardNth3K22 = INITVALUE; - CCTK_REAL PDstandardNth1K23 = INITVALUE; - CCTK_REAL PDstandardNth2K23 = INITVALUE; - CCTK_REAL PDstandardNth3K23 = INITVALUE; - CCTK_REAL PDstandardNth1K33 = INITVALUE; - CCTK_REAL PDstandardNth2K33 = INITVALUE; - CCTK_REAL PDstandardNth3K33 = INITVALUE; + // CCTK_REAL PDstandardNth1alpha = INITVALUE; + // CCTK_REAL PDstandardNth2alpha = INITVALUE; + // CCTK_REAL PDstandardNth3alpha = INITVALUE; + // CCTK_REAL PDstandardNth11alpha = INITVALUE; + // CCTK_REAL PDstandardNth22alpha = INITVALUE; + // CCTK_REAL PDstandardNth33alpha = INITVALUE; + // CCTK_REAL PDstandardNth12alpha = INITVALUE; + // CCTK_REAL PDstandardNth13alpha = INITVALUE; + // CCTK_REAL PDstandardNth23alpha = INITVALUE; + // CCTK_REAL PDstandardNth1beta1 = INITVALUE; + // CCTK_REAL PDstandardNth2beta1 = INITVALUE; + // CCTK_REAL PDstandardNth3beta1 = INITVALUE; + // CCTK_REAL PDstandardNth1beta2 = INITVALUE; + // CCTK_REAL PDstandardNth2beta2 = INITVALUE; + // CCTK_REAL PDstandardNth3beta2 = INITVALUE; + // CCTK_REAL PDstandardNth1beta3 = INITVALUE; + // CCTK_REAL PDstandardNth2beta3 = INITVALUE; + // CCTK_REAL PDstandardNth3beta3 = INITVALUE; + // CCTK_REAL PDstandardNth1g11 = INITVALUE; + // CCTK_REAL PDstandardNth2g11 = INITVALUE; + // CCTK_REAL PDstandardNth3g11 = INITVALUE; + // CCTK_REAL PDstandardNth22g11 = INITVALUE; + // CCTK_REAL PDstandardNth33g11 = INITVALUE; + // CCTK_REAL PDstandardNth23g11 = INITVALUE; + // CCTK_REAL PDstandardNth1g12 = INITVALUE; + // CCTK_REAL PDstandardNth2g12 = INITVALUE; + // CCTK_REAL PDstandardNth3g12 = INITVALUE; + // CCTK_REAL PDstandardNth33g12 = INITVALUE; + // CCTK_REAL PDstandardNth12g12 = INITVALUE; + // CCTK_REAL PDstandardNth13g12 = INITVALUE; + // CCTK_REAL PDstandardNth23g12 = INITVALUE; + // CCTK_REAL PDstandardNth1g13 = INITVALUE; + // CCTK_REAL PDstandardNth2g13 = INITVALUE; + // CCTK_REAL PDstandardNth3g13 = INITVALUE; + // CCTK_REAL PDstandardNth22g13 = INITVALUE; + // CCTK_REAL PDstandardNth12g13 = INITVALUE; + // CCTK_REAL PDstandardNth13g13 = INITVALUE; + // CCTK_REAL PDstandardNth23g13 = INITVALUE; + // CCTK_REAL PDstandardNth1g22 = INITVALUE; + // CCTK_REAL PDstandardNth2g22 = INITVALUE; + // CCTK_REAL PDstandardNth3g22 = INITVALUE; + // CCTK_REAL PDstandardNth11g22 = INITVALUE; + // CCTK_REAL PDstandardNth33g22 = INITVALUE; + // CCTK_REAL PDstandardNth13g22 = INITVALUE; + // CCTK_REAL PDstandardNth1g23 = INITVALUE; + // CCTK_REAL PDstandardNth2g23 = INITVALUE; + // CCTK_REAL PDstandardNth3g23 = INITVALUE; + // CCTK_REAL PDstandardNth11g23 = INITVALUE; + // CCTK_REAL PDstandardNth12g23 = INITVALUE; + // CCTK_REAL PDstandardNth13g23 = INITVALUE; + // CCTK_REAL PDstandardNth23g23 = INITVALUE; + // CCTK_REAL PDstandardNth1g33 = INITVALUE; + // CCTK_REAL PDstandardNth2g33 = INITVALUE; + // CCTK_REAL PDstandardNth3g33 = INITVALUE; + // CCTK_REAL PDstandardNth11g33 = INITVALUE; + // CCTK_REAL PDstandardNth22g33 = INITVALUE; + // CCTK_REAL PDstandardNth12g33 = INITVALUE; + // CCTK_REAL PDstandardNth1K11 = INITVALUE; + // CCTK_REAL PDstandardNth2K11 = INITVALUE; + // CCTK_REAL PDstandardNth3K11 = INITVALUE; + // CCTK_REAL PDstandardNth1K12 = INITVALUE; + // CCTK_REAL PDstandardNth2K12 = INITVALUE; + // CCTK_REAL PDstandardNth3K12 = INITVALUE; + // CCTK_REAL PDstandardNth1K13 = INITVALUE; + // CCTK_REAL PDstandardNth2K13 = INITVALUE; + // CCTK_REAL PDstandardNth3K13 = INITVALUE; + // CCTK_REAL PDstandardNth1K22 = INITVALUE; + // CCTK_REAL PDstandardNth2K22 = INITVALUE; + // CCTK_REAL PDstandardNth3K22 = INITVALUE; + // CCTK_REAL PDstandardNth1K23 = INITVALUE; + // CCTK_REAL PDstandardNth2K23 = INITVALUE; + // CCTK_REAL PDstandardNth3K23 = INITVALUE; + // CCTK_REAL PDstandardNth1K33 = INITVALUE; + // CCTK_REAL PDstandardNth2K33 = INITVALUE; + // CCTK_REAL PDstandardNth3K33 = INITVALUE; /* Assign local copies of grid functions */ - alphaL = alpha[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - g11L = g11[index]; - g12L = g12[index]; - g13L = g13[index]; - g22L = g22[index]; - g23L = g23[index]; - g33L = g33[index]; - K11L = K11[index]; - K12L = K12[index]; - K13L = K13[index]; - K22L = K22[index]; - K23L = K23[index]; - K33L = K33[index]; + CCTK_REAL const alphaL = alpha[index]; + CCTK_REAL const beta1L = beta1[index]; + CCTK_REAL const beta2L = beta2[index]; + CCTK_REAL const beta3L = beta3[index]; + CCTK_REAL const g11L = g11[index]; + CCTK_REAL const g12L = g12[index]; + CCTK_REAL const g13L = g13[index]; + CCTK_REAL const g22L = g22[index]; + CCTK_REAL const g23L = g23[index]; + CCTK_REAL const g33L = g33[index]; + CCTK_REAL const K11L = K11[index]; + CCTK_REAL const K12L = K12[index]; + CCTK_REAL const K13L = K13[index]; + CCTK_REAL const K22L = K22[index]; + CCTK_REAL const K23L = K23[index]; + CCTK_REAL const K33L = K33[index]; /* Assign local copies of subblock grid functions */ /* Include user supplied include files */ /* Precompute derivatives (new style) */ - PDstandardNth1alpha = PDstandardNth1(alpha, i, j, k); - PDstandardNth2alpha = PDstandardNth2(alpha, i, j, k); - PDstandardNth3alpha = PDstandardNth3(alpha, i, j, k); - PDstandardNth11alpha = PDstandardNth11(alpha, i, j, k); - PDstandardNth22alpha = PDstandardNth22(alpha, i, j, k); - PDstandardNth33alpha = PDstandardNth33(alpha, i, j, k); - PDstandardNth12alpha = PDstandardNth12(alpha, i, j, k); - PDstandardNth13alpha = PDstandardNth13(alpha, i, j, k); - PDstandardNth23alpha = PDstandardNth23(alpha, i, j, k); - PDstandardNth1beta1 = PDstandardNth1(beta1, i, j, k); - PDstandardNth2beta1 = PDstandardNth2(beta1, i, j, k); - PDstandardNth3beta1 = PDstandardNth3(beta1, i, j, k); - PDstandardNth1beta2 = PDstandardNth1(beta2, i, j, k); - PDstandardNth2beta2 = PDstandardNth2(beta2, i, j, k); - PDstandardNth3beta2 = PDstandardNth3(beta2, i, j, k); - PDstandardNth1beta3 = PDstandardNth1(beta3, i, j, k); - PDstandardNth2beta3 = PDstandardNth2(beta3, i, j, k); - PDstandardNth3beta3 = PDstandardNth3(beta3, i, j, k); - PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); - PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); - PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); - PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); - PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); - PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); - PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); - PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); - PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); - PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); - PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); - PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); - PDstandardNth23g12 = PDstandardNth23(g12, i, j, k); - PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); - PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); - PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); - PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); - PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); - PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); - PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); - PDstandardNth1g22 = PDstandardNth1(g22, i, j, k); - PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); - PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); - PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); - PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); - PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); - PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); - PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); - PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); - PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); - PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); - PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); - PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); - PDstandardNth1g33 = PDstandardNth1(g33, i, j, k); - PDstandardNth2g33 = PDstandardNth2(g33, i, j, k); - PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); - PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); - PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); - PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); - PDstandardNth1K11 = PDstandardNth1(K11, i, j, k); - PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); - PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); - PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); - PDstandardNth2K12 = PDstandardNth2(K12, i, j, k); - PDstandardNth3K12 = PDstandardNth3(K12, i, j, k); - PDstandardNth1K13 = PDstandardNth1(K13, i, j, k); - PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); - PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); - PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); - PDstandardNth2K22 = PDstandardNth2(K22, i, j, k); - PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); - PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); - PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); - PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); - PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); - PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); - PDstandardNth3K33 = PDstandardNth3(K33, i, j, k); + CCTK_REAL const PDstandardNth1alpha = PDstandardNth1(alpha, i, j, k); + CCTK_REAL const PDstandardNth2alpha = PDstandardNth2(alpha, i, j, k); + CCTK_REAL const PDstandardNth3alpha = PDstandardNth3(alpha, i, j, k); + CCTK_REAL const PDstandardNth11alpha = PDstandardNth11(alpha, i, j, k); + CCTK_REAL const PDstandardNth22alpha = PDstandardNth22(alpha, i, j, k); + CCTK_REAL const PDstandardNth33alpha = PDstandardNth33(alpha, i, j, k); + CCTK_REAL const PDstandardNth12alpha = PDstandardNth12(alpha, i, j, k); + CCTK_REAL const PDstandardNth13alpha = PDstandardNth13(alpha, i, j, k); + CCTK_REAL const PDstandardNth23alpha = PDstandardNth23(alpha, i, j, k); + CCTK_REAL const PDstandardNth1beta1 = PDstandardNth1(beta1, i, j, k); + CCTK_REAL const PDstandardNth2beta1 = PDstandardNth2(beta1, i, j, k); + CCTK_REAL const PDstandardNth3beta1 = PDstandardNth3(beta1, i, j, k); + CCTK_REAL const PDstandardNth1beta2 = PDstandardNth1(beta2, i, j, k); + CCTK_REAL const PDstandardNth2beta2 = PDstandardNth2(beta2, i, j, k); + CCTK_REAL const PDstandardNth3beta2 = PDstandardNth3(beta2, i, j, k); + CCTK_REAL const PDstandardNth1beta3 = PDstandardNth1(beta3, i, j, k); + CCTK_REAL const PDstandardNth2beta3 = PDstandardNth2(beta3, i, j, k); + CCTK_REAL const PDstandardNth3beta3 = PDstandardNth3(beta3, i, j, k); + CCTK_REAL const PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); + CCTK_REAL const PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); + CCTK_REAL const PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); + CCTK_REAL const PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); + CCTK_REAL const PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); + CCTK_REAL const PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); + CCTK_REAL const PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); + CCTK_REAL const PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); + CCTK_REAL const PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); + CCTK_REAL const PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); + CCTK_REAL const PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); + CCTK_REAL const PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); + CCTK_REAL const PDstandardNth23g12 = PDstandardNth23(g12, i, j, k); + CCTK_REAL const PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); + CCTK_REAL const PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); + CCTK_REAL const PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); + CCTK_REAL const PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); + CCTK_REAL const PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); + CCTK_REAL const PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); + CCTK_REAL const PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); + CCTK_REAL const PDstandardNth1g22 = PDstandardNth1(g22, i, j, k); + CCTK_REAL const PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); + CCTK_REAL const PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); + CCTK_REAL const PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); + CCTK_REAL const PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); + CCTK_REAL const PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); + CCTK_REAL const PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); + CCTK_REAL const PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); + CCTK_REAL const PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); + CCTK_REAL const PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); + CCTK_REAL const PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); + CCTK_REAL const PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); + CCTK_REAL const PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); + CCTK_REAL const PDstandardNth1g33 = PDstandardNth1(g33, i, j, k); + CCTK_REAL const PDstandardNth2g33 = PDstandardNth2(g33, i, j, k); + CCTK_REAL const PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); + CCTK_REAL const PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); + CCTK_REAL const PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); + CCTK_REAL const PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); + CCTK_REAL const PDstandardNth1K11 = PDstandardNth1(K11, i, j, k); + CCTK_REAL const PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); + CCTK_REAL const PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); + CCTK_REAL const PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); + CCTK_REAL const PDstandardNth2K12 = PDstandardNth2(K12, i, j, k); + CCTK_REAL const PDstandardNth3K12 = PDstandardNth3(K12, i, j, k); + CCTK_REAL const PDstandardNth1K13 = PDstandardNth1(K13, i, j, k); + CCTK_REAL const PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); + CCTK_REAL const PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); + CCTK_REAL const PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); + CCTK_REAL const PDstandardNth2K22 = PDstandardNth2(K22, i, j, k); + CCTK_REAL const PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); + CCTK_REAL const PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); + CCTK_REAL const PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); + CCTK_REAL const PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); + CCTK_REAL const PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); + CCTK_REAL const PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); + CCTK_REAL const PDstandardNth3K33 = PDstandardNth3(K33, i, j, k); /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); + CCTK_REAL const detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); - gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); + CCTK_REAL const gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); - gu12 = (g13L*g23L - g12L*g33L)*INV(detg); + CCTK_REAL const gu12 = (g13L*g23L - g12L*g33L)*INV(detg); - gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + CCTK_REAL const gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); - gu21 = (g13L*g23L - g12L*g33L)*INV(detg); + CCTK_REAL const gu21 = (g13L*g23L - g12L*g33L)*INV(detg); - gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); + CCTK_REAL const gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); - gu23 = (g12L*g13L - g11L*g23L)*INV(detg); + CCTK_REAL const gu23 = (g12L*g13L - g11L*g23L)*INV(detg); - gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + CCTK_REAL const gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); - gu32 = (g12L*g13L - g11L*g23L)*INV(detg); + CCTK_REAL const gu32 = (g12L*g13L - g11L*g23L)*INV(detg); - gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); + CCTK_REAL const gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); - G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - gu12*PDstandardNth2g11 - + CCTK_REAL const G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11); - G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - gu22*PDstandardNth2g11 - + CCTK_REAL const G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11); - G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - + CCTK_REAL const G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11); - G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + + CCTK_REAL const G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + + CCTK_REAL const G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + + CCTK_REAL const G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + + CCTK_REAL const G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + + CCTK_REAL const G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + + CCTK_REAL const G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu12*PDstandardNth2g22 + + CCTK_REAL const G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + + CCTK_REAL const G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + + CCTK_REAL const G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G123 = khalf*(gu13*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G123 = khalf*(gu13*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu12*PDstandardNth3g22); - G223 = khalf*(gu23*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G223 = khalf*(gu23*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu22*PDstandardNth3g22); - G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu32*PDstandardNth3g22); - G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + + CCTK_REAL const G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + gu13*PDstandardNth3g33); - G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + + CCTK_REAL const G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + gu23*PDstandardNth3g33); - G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + + CCTK_REAL const G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33); - R11 = 2*(G112*G211 + G113*G311 + G213*G312) - G111*(G111 + G212 + G313) - G211*(G112 + G222 + G323) - + CCTK_REAL const R11 = 2*(G112*G211 + G113*G311 + G213*G312) - G111*(G111 + G212 + G313) - G211*(G112 + G222 + G323) - G311*(G113 + G223 + G333) + khalf*(-(gu22*(PDstandardNth11g22 - 2*PDstandardNth12g12 + PDstandardNth22g11)) + gu23*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + gu32*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) - gu33*(PDstandardNth11g33 - 2*PDstandardNth13g13 + PDstandardNth33g11)) + SQR(G111) + SQR(G212) + SQR(G313); - R12 = khalf*(2*(G122*G211 + G123*G311 + G213*G322 + G313*G323) - + CCTK_REAL const R12 = khalf*(2*(G122*G211 + G123*G311 + G213*G322 + G313*G323) - 2*(G112*G212 + G112*G313 + G212*G323 + G312*G333 + gu12*PDstandardNth12g12) - gu32*PDstandardNth12g23 - gu33*PDstandardNth12g33 + gu13*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + gu32*PDstandardNth13g22 + gu33*PDstandardNth13g23 + gu12*(PDstandardNth11g22 + PDstandardNth22g11) + gu32*PDstandardNth22g13 + gu13*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - gu33*PDstandardNth33g12); - R13 = khalf*(2*(G123*G211 + G212*G223 + G133*G311 + G233*G312) - + CCTK_REAL const R13 = khalf*(2*(G123*G211 + G212*G223 + G133*G311 + G233*G312) - 2*(G213*G222 + G223*G313 + G113*(G212 + G313) + gu13*PDstandardNth13g13) + gu12*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + gu13*(PDstandardNth11g33 + PDstandardNth33g11) + gu23*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12)); - R22 = -(G122*(G111 + G212 + G313)) + 2*(G122*G212 + G123*G312 + G223*G322) - G222*(G112 + G222 + G323) - + CCTK_REAL const R22 = -(G122*(G111 + G212 + G313)) + 2*(G122*G212 + G123*G312 + G223*G322) - G222*(G112 + G222 + G323) - G322*(G113 + G223 + G333) + khalf*(-(gu11*(PDstandardNth11g22 - 2*PDstandardNth12g12 + PDstandardNth22g11)) + gu13*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + gu31*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) - gu33*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) + SQR(G112) + SQR(G222) + SQR(G323); - R23 = khalf*(2*(G112*G113 + G122*G213 + G133*G312 + G233*G322) + + CCTK_REAL const R23 = khalf*(2*(G112*G113 + G122*G213 + G133*G312 + G233*G322) + gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + PDstandardNth22g13 - PDstandardNth23g12) - 2*(G111*G123 + G113*G323 + G223*(G112 + G323) + gu23*PDstandardNth23g23) + gu13*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + gu23*(PDstandardNth22g33 + PDstandardNth33g22)); - R33 = -(G133*(G111 + G212 + G313)) + 2*(G123*G213 + G133*G313) + 2*G233*G323 - G233*(G112 + G222 + G323) - + CCTK_REAL const R33 = -(G133*(G111 + G212 + G313)) + 2*(G123*G213 + G133*G313) + 2*G233*G323 - G233*(G112 + G222 + G323) - G333*(G113 + G223 + G333) + khalf*(-(gu11*(PDstandardNth11g33 - 2*PDstandardNth13g13 + PDstandardNth33g11)) + gu12*(-PDstandardNth12g33 + PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) + gu21*(-PDstandardNth12g33 + PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) - gu22*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) + SQR(G113) + SQR(G223) + SQR(G333); - Km11 = gu11*K11L + gu12*K12L + gu13*K13L; + CCTK_REAL const Km11 = gu11*K11L + gu12*K12L + gu13*K13L; - Km21 = gu21*K11L + gu22*K12L + gu23*K13L; + CCTK_REAL const Km21 = gu21*K11L + gu22*K12L + gu23*K13L; - Km31 = gu31*K11L + gu32*K12L + gu33*K13L; + CCTK_REAL const Km31 = gu31*K11L + gu32*K12L + gu33*K13L; - Km12 = gu11*K12L + gu12*K22L + gu13*K23L; + CCTK_REAL const Km12 = gu11*K12L + gu12*K22L + gu13*K23L; - Km22 = gu21*K12L + gu22*K22L + gu23*K23L; + CCTK_REAL const Km22 = gu21*K12L + gu22*K22L + gu23*K23L; - Km32 = gu31*K12L + gu32*K22L + gu33*K23L; + CCTK_REAL const Km32 = gu31*K12L + gu32*K22L + gu33*K23L; - Km13 = gu11*K13L + gu12*K23L + gu13*K33L; + CCTK_REAL const Km13 = gu11*K13L + gu12*K23L + gu13*K33L; - Km23 = gu21*K13L + gu22*K23L + gu23*K33L; + CCTK_REAL const Km23 = gu21*K13L + gu22*K23L + gu23*K33L; - Km33 = gu31*K13L + gu32*K23L + gu33*K33L; + CCTK_REAL const Km33 = gu31*K13L + gu32*K23L + gu33*K33L; - trK = Km11 + Km22 + Km33; + CCTK_REAL const trK = Km11 + Km22 + Km33; - g11rhsL = -2*alphaL*K11L + 2*(g11L*PDstandardNth1beta1 + g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + + CCTK_REAL const g11rhsL = -2*alphaL*K11L + 2*(g11L*PDstandardNth1beta1 + g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g11 + beta3L*PDstandardNth3g11; - g12rhsL = -2*alphaL*K12L + g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3 + beta1L*PDstandardNth1g12 + + CCTK_REAL const g12rhsL = -2*alphaL*K12L + g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3 + beta1L*PDstandardNth1g12 + g11L*PDstandardNth2beta1 + g12L*(PDstandardNth1beta1 + PDstandardNth2beta2) + g13L*PDstandardNth2beta3 + beta2L*PDstandardNth2g12 + beta3L*PDstandardNth3g12; - g13rhsL = -2*alphaL*K13L + g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3 + beta1L*PDstandardNth1g13 + + CCTK_REAL const g13rhsL = -2*alphaL*K13L + g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3 + beta1L*PDstandardNth1g13 + beta2L*PDstandardNth2g13 + g11L*PDstandardNth3beta1 + g12L*PDstandardNth3beta2 + g13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + beta3L*PDstandardNth3g13; - g22rhsL = -2*alphaL*K22L + beta1L*PDstandardNth1g22 + + CCTK_REAL const g22rhsL = -2*alphaL*K22L + beta1L*PDstandardNth1g22 + 2*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + beta3L*PDstandardNth3g22; - g23rhsL = -2*alphaL*K23L + beta1L*PDstandardNth1g23 + g13L*PDstandardNth2beta1 + g33L*PDstandardNth2beta3 + + CCTK_REAL const g23rhsL = -2*alphaL*K23L + beta1L*PDstandardNth1g23 + g13L*PDstandardNth2beta1 + g33L*PDstandardNth2beta3 + beta2L*PDstandardNth2g23 + g12L*PDstandardNth3beta1 + g22L*PDstandardNth3beta2 + g23L*(PDstandardNth2beta2 + PDstandardNth3beta3) + beta3L*PDstandardNth3g23; - g33rhsL = -2*alphaL*K33L + beta1L*PDstandardNth1g33 + beta2L*PDstandardNth2g33 + + CCTK_REAL const g33rhsL = -2*alphaL*K33L + beta1L*PDstandardNth1g33 + beta2L*PDstandardNth2g33 + 2*(g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + g33L*PDstandardNth3beta3) + beta3L*PDstandardNth3g33; - K11rhsL = -PDstandardNth11alpha + G111*PDstandardNth1alpha + + CCTK_REAL const K11rhsL = -PDstandardNth11alpha + G111*PDstandardNth1alpha + 2*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + G211*PDstandardNth2alpha + beta2L*PDstandardNth2K11 + G311*PDstandardNth3alpha + beta3L*PDstandardNth3K11 + alphaL*(-2*(K11L*Km11 + K12L*Km21 + K13L*Km31) + R11 + K11L*trK); - K12rhsL = -PDstandardNth12alpha + G112*PDstandardNth1alpha + K22L*PDstandardNth1beta2 + K23L*PDstandardNth1beta3 + + CCTK_REAL const K12rhsL = -PDstandardNth12alpha + G112*PDstandardNth1alpha + K22L*PDstandardNth1beta2 + K23L*PDstandardNth1beta3 + beta1L*PDstandardNth1K12 + G212*PDstandardNth2alpha + K11L*PDstandardNth2beta1 + K12L*(PDstandardNth1beta1 + PDstandardNth2beta2) + K13L*PDstandardNth2beta3 + beta2L*PDstandardNth2K12 + G312*PDstandardNth3alpha + beta3L*PDstandardNth3K12 + alphaL*(-2*(K11L*Km12 + K12L*Km22 + K13L*Km32) + R12 + K12L*trK); - K13rhsL = -PDstandardNth13alpha + G113*PDstandardNth1alpha + K23L*PDstandardNth1beta2 + K33L*PDstandardNth1beta3 + + CCTK_REAL const K13rhsL = -PDstandardNth13alpha + G113*PDstandardNth1alpha + K23L*PDstandardNth1beta2 + K33L*PDstandardNth1beta3 + beta1L*PDstandardNth1K13 + G213*PDstandardNth2alpha + beta2L*PDstandardNth2K13 + G313*PDstandardNth3alpha + K11L*PDstandardNth3beta1 + K12L*PDstandardNth3beta2 + K13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + beta3L*PDstandardNth3K13 + alphaL*(-2*(K11L*Km13 + K12L*Km23 + K13L*Km33) + R13 + K13L*trK); - K22rhsL = G122*PDstandardNth1alpha + beta1L*PDstandardNth1K22 - PDstandardNth22alpha + G222*PDstandardNth2alpha + + CCTK_REAL const K22rhsL = G122*PDstandardNth1alpha + beta1L*PDstandardNth1K22 - PDstandardNth22alpha + G222*PDstandardNth2alpha + 2*(K12L*PDstandardNth2beta1 + K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3) + beta2L*PDstandardNth2K22 + G322*PDstandardNth3alpha + beta3L*PDstandardNth3K22 + alphaL*(-2*(K12L*Km12 + K22L*Km22 + K23L*Km32) + R22 + K22L*trK); - K23rhsL = G123*PDstandardNth1alpha + beta1L*PDstandardNth1K23 - PDstandardNth23alpha + G223*PDstandardNth2alpha + + CCTK_REAL const K23rhsL = G123*PDstandardNth1alpha + beta1L*PDstandardNth1K23 - PDstandardNth23alpha + G223*PDstandardNth2alpha + K13L*PDstandardNth2beta1 + K33L*PDstandardNth2beta3 + beta2L*PDstandardNth2K23 + G323*PDstandardNth3alpha + K12L*PDstandardNth3beta1 + K22L*PDstandardNth3beta2 + K23L*(PDstandardNth2beta2 + PDstandardNth3beta3) + beta3L*PDstandardNth3K23 + alphaL*(-2*(K12L*Km13 + K22L*Km23 + K23L*Km33) + R23 + K23L*trK); - K33rhsL = G133*PDstandardNth1alpha + beta1L*PDstandardNth1K33 + G233*PDstandardNth2alpha + beta2L*PDstandardNth2K33 - + CCTK_REAL const K33rhsL = G133*PDstandardNth1alpha + beta1L*PDstandardNth1K33 + G233*PDstandardNth2alpha + beta2L*PDstandardNth2K33 - PDstandardNth33alpha + G333*PDstandardNth3alpha + 2*(K13L*PDstandardNth3beta1 + K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3) + beta3L*PDstandardNth3K33 + alphaL*(-2*(K13L*Km13 + K23L*Km23 + K33L*Km33) + R33 + K33L*trK); - alpharhsL = 0; + CCTK_REAL const alpharhsL = 0; - beta1rhsL = 0; + CCTK_REAL const beta1rhsL = 0; - beta2rhsL = 0; + CCTK_REAL const beta2rhsL = 0; - beta3rhsL = 0; + CCTK_REAL const beta3rhsL = 0; /* Copy local copies back to grid functions */ diff --git a/ML_ADM/src/ML_ADM_boundary.c b/ML_ADM/src/ML_ADM_boundary.c index 8270f22..a177e52 100644 --- a/ML_ADM/src/ML_ADM_boundary.c +++ b/ML_ADM/src/ML_ADM_boundary.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_boundary_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_boundary_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_boundary_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,18 +87,18 @@ void ML_ADM_boundary_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ /* Declare local copies of grid functions */ - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; - CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + // CCTK_REAL alphaL = INITVALUE; + // CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ @@ -116,37 +114,37 @@ void ML_ADM_boundary_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - g11L = 1; + CCTK_REAL const g11L = 1; - g12L = 0; + CCTK_REAL const g12L = 0; - g13L = 0; + CCTK_REAL const g13L = 0; - g22L = 1; + CCTK_REAL const g22L = 1; - g23L = 0; + CCTK_REAL const g23L = 0; - g33L = 1; + CCTK_REAL const g33L = 1; - K11L = 0; + CCTK_REAL const K11L = 0; - K12L = 0; + CCTK_REAL const K12L = 0; - K13L = 0; + CCTK_REAL const K13L = 0; - K22L = 0; + CCTK_REAL const K22L = 0; - K23L = 0; + CCTK_REAL const K23L = 0; - K33L = 0; + CCTK_REAL const K33L = 0; - alphaL = 1; + CCTK_REAL const alphaL = 1; - beta1L = 0; + CCTK_REAL const beta1L = 0; - beta2L = 0; + CCTK_REAL const beta2L = 0; - beta3L = 0; + CCTK_REAL const beta3L = 0; /* Copy local copies back to grid functions */ diff --git a/ML_ADM/src/ML_ADM_constraints.c b/ML_ADM/src/ML_ADM_constraints.c index f508db2..4b9338e 100644 --- a/ML_ADM/src/ML_ADM_constraints.c +++ b/ML_ADM/src/ML_ADM_constraints.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_constraints_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_ /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,303 +87,303 @@ void ML_ADM_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_ i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ - CCTK_REAL detg = INITVALUE; - CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; - CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; - CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; - CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; - CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; - CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; - CCTK_REAL trK = INITVALUE; - CCTK_REAL trR = INITVALUE; + // CCTK_REAL detg = INITVALUE; + // CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; + // CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; + // CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; + // CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; + // CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; + // CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; + // CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; + // CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; + // CCTK_REAL trK = INITVALUE; + // CCTK_REAL trR = INITVALUE; /* Declare local copies of grid functions */ - CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; - CCTK_REAL HL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; - CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + // CCTK_REAL HL = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + // CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ - CCTK_REAL PDstandardNth1g11 = INITVALUE; - CCTK_REAL PDstandardNth2g11 = INITVALUE; - CCTK_REAL PDstandardNth3g11 = INITVALUE; - CCTK_REAL PDstandardNth22g11 = INITVALUE; - CCTK_REAL PDstandardNth33g11 = INITVALUE; - CCTK_REAL PDstandardNth23g11 = INITVALUE; - CCTK_REAL PDstandardNth1g12 = INITVALUE; - CCTK_REAL PDstandardNth2g12 = INITVALUE; - CCTK_REAL PDstandardNth3g12 = INITVALUE; - CCTK_REAL PDstandardNth33g12 = INITVALUE; - CCTK_REAL PDstandardNth12g12 = INITVALUE; - CCTK_REAL PDstandardNth13g12 = INITVALUE; - CCTK_REAL PDstandardNth23g12 = INITVALUE; - CCTK_REAL PDstandardNth1g13 = INITVALUE; - CCTK_REAL PDstandardNth2g13 = INITVALUE; - CCTK_REAL PDstandardNth3g13 = INITVALUE; - CCTK_REAL PDstandardNth22g13 = INITVALUE; - CCTK_REAL PDstandardNth12g13 = INITVALUE; - CCTK_REAL PDstandardNth13g13 = INITVALUE; - CCTK_REAL PDstandardNth23g13 = INITVALUE; - CCTK_REAL PDstandardNth1g22 = INITVALUE; - CCTK_REAL PDstandardNth2g22 = INITVALUE; - CCTK_REAL PDstandardNth3g22 = INITVALUE; - CCTK_REAL PDstandardNth11g22 = INITVALUE; - CCTK_REAL PDstandardNth33g22 = INITVALUE; - CCTK_REAL PDstandardNth13g22 = INITVALUE; - CCTK_REAL PDstandardNth1g23 = INITVALUE; - CCTK_REAL PDstandardNth2g23 = INITVALUE; - CCTK_REAL PDstandardNth3g23 = INITVALUE; - CCTK_REAL PDstandardNth11g23 = INITVALUE; - CCTK_REAL PDstandardNth12g23 = INITVALUE; - CCTK_REAL PDstandardNth13g23 = INITVALUE; - CCTK_REAL PDstandardNth23g23 = INITVALUE; - CCTK_REAL PDstandardNth1g33 = INITVALUE; - CCTK_REAL PDstandardNth2g33 = INITVALUE; - CCTK_REAL PDstandardNth3g33 = INITVALUE; - CCTK_REAL PDstandardNth11g33 = INITVALUE; - CCTK_REAL PDstandardNth22g33 = INITVALUE; - CCTK_REAL PDstandardNth12g33 = INITVALUE; - CCTK_REAL PDstandardNth2K11 = INITVALUE; - CCTK_REAL PDstandardNth3K11 = INITVALUE; - CCTK_REAL PDstandardNth1K12 = INITVALUE; - CCTK_REAL PDstandardNth2K12 = INITVALUE; - CCTK_REAL PDstandardNth3K12 = INITVALUE; - CCTK_REAL PDstandardNth1K13 = INITVALUE; - CCTK_REAL PDstandardNth2K13 = INITVALUE; - CCTK_REAL PDstandardNth3K13 = INITVALUE; - CCTK_REAL PDstandardNth1K22 = INITVALUE; - CCTK_REAL PDstandardNth3K22 = INITVALUE; - CCTK_REAL PDstandardNth1K23 = INITVALUE; - CCTK_REAL PDstandardNth2K23 = INITVALUE; - CCTK_REAL PDstandardNth3K23 = INITVALUE; - CCTK_REAL PDstandardNth1K33 = INITVALUE; - CCTK_REAL PDstandardNth2K33 = INITVALUE; + // CCTK_REAL PDstandardNth1g11 = INITVALUE; + // CCTK_REAL PDstandardNth2g11 = INITVALUE; + // CCTK_REAL PDstandardNth3g11 = INITVALUE; + // CCTK_REAL PDstandardNth22g11 = INITVALUE; + // CCTK_REAL PDstandardNth33g11 = INITVALUE; + // CCTK_REAL PDstandardNth23g11 = INITVALUE; + // CCTK_REAL PDstandardNth1g12 = INITVALUE; + // CCTK_REAL PDstandardNth2g12 = INITVALUE; + // CCTK_REAL PDstandardNth3g12 = INITVALUE; + // CCTK_REAL PDstandardNth33g12 = INITVALUE; + // CCTK_REAL PDstandardNth12g12 = INITVALUE; + // CCTK_REAL PDstandardNth13g12 = INITVALUE; + // CCTK_REAL PDstandardNth23g12 = INITVALUE; + // CCTK_REAL PDstandardNth1g13 = INITVALUE; + // CCTK_REAL PDstandardNth2g13 = INITVALUE; + // CCTK_REAL PDstandardNth3g13 = INITVALUE; + // CCTK_REAL PDstandardNth22g13 = INITVALUE; + // CCTK_REAL PDstandardNth12g13 = INITVALUE; + // CCTK_REAL PDstandardNth13g13 = INITVALUE; + // CCTK_REAL PDstandardNth23g13 = INITVALUE; + // CCTK_REAL PDstandardNth1g22 = INITVALUE; + // CCTK_REAL PDstandardNth2g22 = INITVALUE; + // CCTK_REAL PDstandardNth3g22 = INITVALUE; + // CCTK_REAL PDstandardNth11g22 = INITVALUE; + // CCTK_REAL PDstandardNth33g22 = INITVALUE; + // CCTK_REAL PDstandardNth13g22 = INITVALUE; + // CCTK_REAL PDstandardNth1g23 = INITVALUE; + // CCTK_REAL PDstandardNth2g23 = INITVALUE; + // CCTK_REAL PDstandardNth3g23 = INITVALUE; + // CCTK_REAL PDstandardNth11g23 = INITVALUE; + // CCTK_REAL PDstandardNth12g23 = INITVALUE; + // CCTK_REAL PDstandardNth13g23 = INITVALUE; + // CCTK_REAL PDstandardNth23g23 = INITVALUE; + // CCTK_REAL PDstandardNth1g33 = INITVALUE; + // CCTK_REAL PDstandardNth2g33 = INITVALUE; + // CCTK_REAL PDstandardNth3g33 = INITVALUE; + // CCTK_REAL PDstandardNth11g33 = INITVALUE; + // CCTK_REAL PDstandardNth22g33 = INITVALUE; + // CCTK_REAL PDstandardNth12g33 = INITVALUE; + // CCTK_REAL PDstandardNth2K11 = INITVALUE; + // CCTK_REAL PDstandardNth3K11 = INITVALUE; + // CCTK_REAL PDstandardNth1K12 = INITVALUE; + // CCTK_REAL PDstandardNth2K12 = INITVALUE; + // CCTK_REAL PDstandardNth3K12 = INITVALUE; + // CCTK_REAL PDstandardNth1K13 = INITVALUE; + // CCTK_REAL PDstandardNth2K13 = INITVALUE; + // CCTK_REAL PDstandardNth3K13 = INITVALUE; + // CCTK_REAL PDstandardNth1K22 = INITVALUE; + // CCTK_REAL PDstandardNth3K22 = INITVALUE; + // CCTK_REAL PDstandardNth1K23 = INITVALUE; + // CCTK_REAL PDstandardNth2K23 = INITVALUE; + // CCTK_REAL PDstandardNth3K23 = INITVALUE; + // CCTK_REAL PDstandardNth1K33 = INITVALUE; + // CCTK_REAL PDstandardNth2K33 = INITVALUE; /* Assign local copies of grid functions */ - g11L = g11[index]; - g12L = g12[index]; - g13L = g13[index]; - g22L = g22[index]; - g23L = g23[index]; - g33L = g33[index]; - K11L = K11[index]; - K12L = K12[index]; - K13L = K13[index]; - K22L = K22[index]; - K23L = K23[index]; - K33L = K33[index]; + CCTK_REAL const g11L = g11[index]; + CCTK_REAL const g12L = g12[index]; + CCTK_REAL const g13L = g13[index]; + CCTK_REAL const g22L = g22[index]; + CCTK_REAL const g23L = g23[index]; + CCTK_REAL const g33L = g33[index]; + CCTK_REAL const K11L = K11[index]; + CCTK_REAL const K12L = K12[index]; + CCTK_REAL const K13L = K13[index]; + CCTK_REAL const K22L = K22[index]; + CCTK_REAL const K23L = K23[index]; + CCTK_REAL const K33L = K33[index]; /* Assign local copies of subblock grid functions */ /* Include user supplied include files */ /* Precompute derivatives (new style) */ - PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); - PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); - PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); - PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); - PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); - PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); - PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); - PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); - PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); - PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); - PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); - PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); - PDstandardNth23g12 = PDstandardNth23(g12, i, j, k); - PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); - PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); - PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); - PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); - PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); - PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); - PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); - PDstandardNth1g22 = PDstandardNth1(g22, i, j, k); - PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); - PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); - PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); - PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); - PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); - PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); - PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); - PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); - PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); - PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); - PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); - PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); - PDstandardNth1g33 = PDstandardNth1(g33, i, j, k); - PDstandardNth2g33 = PDstandardNth2(g33, i, j, k); - PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); - PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); - PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); - PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); - PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); - PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); - PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); - PDstandardNth2K12 = PDstandardNth2(K12, i, j, k); - PDstandardNth3K12 = PDstandardNth3(K12, i, j, k); - PDstandardNth1K13 = PDstandardNth1(K13, i, j, k); - PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); - PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); - PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); - PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); - PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); - PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); - PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); - PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); - PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); + CCTK_REAL const PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); + CCTK_REAL const PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); + CCTK_REAL const PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); + CCTK_REAL const PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); + CCTK_REAL const PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); + CCTK_REAL const PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); + CCTK_REAL const PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); + CCTK_REAL const PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); + CCTK_REAL const PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); + CCTK_REAL const PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); + CCTK_REAL const PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); + CCTK_REAL const PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); + CCTK_REAL const PDstandardNth23g12 = PDstandardNth23(g12, i, j, k); + CCTK_REAL const PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); + CCTK_REAL const PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); + CCTK_REAL const PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); + CCTK_REAL const PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); + CCTK_REAL const PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); + CCTK_REAL const PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); + CCTK_REAL const PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); + CCTK_REAL const PDstandardNth1g22 = PDstandardNth1(g22, i, j, k); + CCTK_REAL const PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); + CCTK_REAL const PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); + CCTK_REAL const PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); + CCTK_REAL const PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); + CCTK_REAL const PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); + CCTK_REAL const PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); + CCTK_REAL const PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); + CCTK_REAL const PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); + CCTK_REAL const PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); + CCTK_REAL const PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); + CCTK_REAL const PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); + CCTK_REAL const PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); + CCTK_REAL const PDstandardNth1g33 = PDstandardNth1(g33, i, j, k); + CCTK_REAL const PDstandardNth2g33 = PDstandardNth2(g33, i, j, k); + CCTK_REAL const PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); + CCTK_REAL const PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); + CCTK_REAL const PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); + CCTK_REAL const PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); + CCTK_REAL const PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); + CCTK_REAL const PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); + CCTK_REAL const PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); + CCTK_REAL const PDstandardNth2K12 = PDstandardNth2(K12, i, j, k); + CCTK_REAL const PDstandardNth3K12 = PDstandardNth3(K12, i, j, k); + CCTK_REAL const PDstandardNth1K13 = PDstandardNth1(K13, i, j, k); + CCTK_REAL const PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); + CCTK_REAL const PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); + CCTK_REAL const PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); + CCTK_REAL const PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); + CCTK_REAL const PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); + CCTK_REAL const PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); + CCTK_REAL const PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); + CCTK_REAL const PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); + CCTK_REAL const PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); + CCTK_REAL const detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); - gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); + CCTK_REAL const gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); - gu12 = (g13L*g23L - g12L*g33L)*INV(detg); + CCTK_REAL const gu12 = (g13L*g23L - g12L*g33L)*INV(detg); - gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + CCTK_REAL const gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); - gu21 = (g13L*g23L - g12L*g33L)*INV(detg); + CCTK_REAL const gu21 = (g13L*g23L - g12L*g33L)*INV(detg); - gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); + CCTK_REAL const gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); - gu23 = (g12L*g13L - g11L*g23L)*INV(detg); + CCTK_REAL const gu23 = (g12L*g13L - g11L*g23L)*INV(detg); - gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + CCTK_REAL const gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); - gu32 = (g12L*g13L - g11L*g23L)*INV(detg); + CCTK_REAL const gu32 = (g12L*g13L - g11L*g23L)*INV(detg); - gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); + CCTK_REAL const gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); - G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - gu12*PDstandardNth2g11 - + CCTK_REAL const G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11); - G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - gu22*PDstandardNth2g11 - + CCTK_REAL const G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11); - G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - + CCTK_REAL const G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11); - G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + + CCTK_REAL const G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + + CCTK_REAL const G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + + CCTK_REAL const G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + + CCTK_REAL const G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + + CCTK_REAL const G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + + CCTK_REAL const G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu12*PDstandardNth2g22 + + CCTK_REAL const G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + + CCTK_REAL const G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + + CCTK_REAL const G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G123 = khalf*(gu13*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G123 = khalf*(gu13*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu12*PDstandardNth3g22); - G223 = khalf*(gu23*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G223 = khalf*(gu23*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu22*PDstandardNth3g22); - G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + CCTK_REAL const G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu32*PDstandardNth3g22); - G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + + CCTK_REAL const G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + gu13*PDstandardNth3g33); - G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + + CCTK_REAL const G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + gu23*PDstandardNth3g33); - G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + + CCTK_REAL const G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33); - R11 = -(G111*(G111 + G122 + G133)) - G211*(G211 + G222 + G233) - G311*(G311 + G322 + G333) + + CCTK_REAL const R11 = -(G111*(G111 + G122 + G133)) - G211*(G211 + G222 + G233) - G311*(G311 + G322 + G333) + khalf*(-(gu22*(PDstandardNth11g22 - 2*PDstandardNth12g12 + PDstandardNth22g11)) + gu23*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + gu32*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) - gu33*(PDstandardNth11g33 - 2*PDstandardNth13g13 + PDstandardNth33g11)) + SQR(G111) + SQR(G112) + SQR(G113) + SQR(G211) + SQR(G212) + SQR(G213) + SQR(G311) + SQR(G312) + SQR(G313); - R12 = khalf*(2*(G113*G123 + G213*G223 + G313*G323) - 2*(G112*G133 + G212*G233 + G312*G333 + gu12*PDstandardNth12g12) - + CCTK_REAL const R12 = khalf*(2*(G113*G123 + G213*G223 + G313*G323) - 2*(G112*G133 + G212*G233 + G312*G333 + gu12*PDstandardNth12g12) - gu32*PDstandardNth12g23 - gu33*PDstandardNth12g33 + gu13*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + gu32*PDstandardNth13g22 + gu33*PDstandardNth13g23 + gu12*(PDstandardNth11g22 + PDstandardNth22g11) + gu32*PDstandardNth22g13 + gu13*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - gu33*PDstandardNth33g12); - R13 = khalf*(2*(G112*G123 + G212*G223 + G312*G323) - 2*(G113*G122 + G213*G222 + G313*G322 + gu13*PDstandardNth13g13) + + CCTK_REAL const R13 = khalf*(2*(G112*G123 + G212*G223 + G312*G323) - 2*(G113*G122 + G213*G222 + G313*G322 + gu13*PDstandardNth13g13) + gu12*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + gu13*(PDstandardNth11g33 + PDstandardNth33g11) + gu23*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12)); - R22 = -(G122*(G111 + G122 + G133)) - G222*(G211 + G222 + G233) - G322*(G311 + G322 + G333) + + CCTK_REAL const R22 = -(G122*(G111 + G122 + G133)) - G222*(G211 + G222 + G233) - G322*(G311 + G322 + G333) + khalf*(-(gu11*(PDstandardNth11g22 - 2*PDstandardNth12g12 + PDstandardNth22g11)) + gu13*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + gu31*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) - gu33*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) + SQR(G112) + SQR(G122) + SQR(G123) + SQR(G212) + SQR(G222) + SQR(G223) + SQR(G312) + SQR(G322) + SQR(G323); - R23 = khalf*(2*(G112*G113 + G212*G213 + G312*G313) + + CCTK_REAL const R23 = khalf*(2*(G112*G113 + G212*G213 + G312*G313) + gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + PDstandardNth22g13 - PDstandardNth23g12) - 2*(G111*G123 + G211*G223 + G311*G323 + gu23*PDstandardNth23g23) + gu13*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + gu23*(PDstandardNth22g33 + PDstandardNth33g22)); - R33 = -(G133*(G111 + G122 + G133)) - G233*(G211 + G222 + G233) - G333*(G311 + G322 + G333) + + CCTK_REAL const R33 = -(G133*(G111 + G122 + G133)) - G233*(G211 + G222 + G233) - G333*(G311 + G322 + G333) + khalf*(-(gu11*(PDstandardNth11g33 - 2*PDstandardNth13g13 + PDstandardNth33g11)) + gu12*(-PDstandardNth12g33 + PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) + gu21*(-PDstandardNth12g33 + PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) - gu22*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) + SQR(G113) + SQR(G123) + SQR(G133) + SQR(G213) + SQR(G223) + SQR(G233) + SQR(G313) + SQR(G323) + SQR(G333); - trR = gu11*R11 + (gu12 + gu21)*R12 + (gu13 + gu31)*R13 + gu22*R22 + (gu23 + gu32)*R23 + gu33*R33; + CCTK_REAL const trR = gu11*R11 + (gu12 + gu21)*R12 + (gu13 + gu31)*R13 + gu22*R22 + (gu23 + gu32)*R23 + gu33*R33; - Km11 = gu11*K11L + gu12*K12L + gu13*K13L; + CCTK_REAL const Km11 = gu11*K11L + gu12*K12L + gu13*K13L; - Km21 = gu21*K11L + gu22*K12L + gu23*K13L; + CCTK_REAL const Km21 = gu21*K11L + gu22*K12L + gu23*K13L; - Km31 = gu31*K11L + gu32*K12L + gu33*K13L; + CCTK_REAL const Km31 = gu31*K11L + gu32*K12L + gu33*K13L; - Km12 = gu11*K12L + gu12*K22L + gu13*K23L; + CCTK_REAL const Km12 = gu11*K12L + gu12*K22L + gu13*K23L; - Km22 = gu21*K12L + gu22*K22L + gu23*K23L; + CCTK_REAL const Km22 = gu21*K12L + gu22*K22L + gu23*K23L; - Km32 = gu31*K12L + gu32*K22L + gu33*K23L; + CCTK_REAL const Km32 = gu31*K12L + gu32*K22L + gu33*K23L; - Km13 = gu11*K13L + gu12*K23L + gu13*K33L; + CCTK_REAL const Km13 = gu11*K13L + gu12*K23L + gu13*K33L; - Km23 = gu21*K13L + gu22*K23L + gu23*K33L; + CCTK_REAL const Km23 = gu21*K13L + gu22*K23L + gu23*K33L; - Km33 = gu31*K13L + gu32*K23L + gu33*K33L; + CCTK_REAL const Km33 = gu31*K13L + gu32*K23L + gu33*K33L; - trK = Km11 + Km22 + Km33; + CCTK_REAL const trK = Km11 + Km22 + Km33; - HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK); + CCTK_REAL const HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK); - M1L = gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + G311*K23L - PDstandardNth1K12 + + CCTK_REAL const M1L = gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + G311*K23L - PDstandardNth1K12 + PDstandardNth2K11) + gu22*(-(G122*K11L) + G112*K12L - G222*K12L - G322*K13L + G212*K22L + G312*K23L - PDstandardNth1K22 + PDstandardNth2K12) + gu23* (-(G123*K11L) + G113*K12L - G223*K12L - G323*K13L + G213*K22L + G313*K23L - PDstandardNth1K23 + PDstandardNth2K13)\ @@ -394,7 +392,7 @@ void ML_ADM_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_ PDstandardNth1K23 + PDstandardNth3K12) + gu33* (-(G133*K11L) - G233*K12L + G113*K13L - G333*K13L + G213*K23L + G313*K33L - PDstandardNth1K33 + PDstandardNth3K13); - M2L = gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - G311*K23L + PDstandardNth1K12 - + CCTK_REAL const M2L = gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - G311*K23L + PDstandardNth1K12 - PDstandardNth2K11) + gu12*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - G312*K23L + PDstandardNth1K22 - PDstandardNth2K12) + gu13* (G123*K11L + (-G113 + G223)*K12L + G323*K13L - G213*K22L - G313*K23L + PDstandardNth1K23 - PDstandardNth2K13) + @@ -403,7 +401,7 @@ void ML_ADM_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_ PDstandardNth2K23 + PDstandardNth3K22) + gu33* (-(G133*K12L) + G123*K13L - G233*K22L + G223*K23L - G333*K23L + G323*K33L - PDstandardNth2K33 + PDstandardNth3K23); - M3L = gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - G311*K33L + PDstandardNth1K13 - + CCTK_REAL const M3L = gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - G311*K33L + PDstandardNth1K13 - PDstandardNth3K11) + gu12*(G123*K11L + G223*K12L + (-G112 + G323)*K13L - G212*K23L - G312*K33L + PDstandardNth1K23 - PDstandardNth3K12) + gu21* (G113*K12L - G112*K13L + G213*K22L - G212*K23L + G313*K23L - G312*K33L + PDstandardNth2K13 - PDstandardNth3K12) + diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.c b/ML_ADM/src/ML_ADM_constraints_boundary.c index 3a3afa0..fc7a5d4 100644 --- a/ML_ADM/src/ML_ADM_constraints_boundary.c +++ b/ML_ADM/src/ML_ADM_constraints_boundary.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_constraints_boundary_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_constraints_boundary_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_constraints_boundary_Body(cGH const * const cctkGH, CCTK_INT const d /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,16 +87,16 @@ void ML_ADM_constraints_boundary_Body(cGH const * const cctkGH, CCTK_INT const d i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ /* Declare local copies of grid functions */ - CCTK_REAL HL = INITVALUE; - CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; + // CCTK_REAL HL = INITVALUE; + // CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ @@ -114,13 +112,13 @@ void ML_ADM_constraints_boundary_Body(cGH const * const cctkGH, CCTK_INT const d /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - HL = 0; + CCTK_REAL const HL = 0; - M1L = 0; + CCTK_REAL const M1L = 0; - M2L = 0; + CCTK_REAL const M2L = 0; - M3L = 0; + CCTK_REAL const M3L = 0; /* Copy local copies back to grid functions */ diff --git a/ML_ADM/src/ML_ADM_convertFromADMBase.c b/ML_ADM/src/ML_ADM_convertFromADMBase.c index fd774c7..c2b7b5f 100644 --- a/ML_ADM/src/ML_ADM_convertFromADMBase.c +++ b/ML_ADM/src/ML_ADM_convertFromADMBase.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_convertFromADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_convertFromADMBase_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_convertFromADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,55 +87,55 @@ void ML_ADM_convertFromADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ /* Declare local copies of grid functions */ - CCTK_REAL alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; + // CCTK_REAL alpL = INITVALUE; + // CCTK_REAL alphaL = INITVALUE; + // CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; + // CCTK_REAL betaxL = INITVALUE; + // CCTK_REAL betayL = INITVALUE; + // CCTK_REAL betazL = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + // CCTK_REAL gxxL = INITVALUE; + // CCTK_REAL gxyL = INITVALUE; + // CCTK_REAL gxzL = INITVALUE; + // CCTK_REAL gyyL = INITVALUE; + // CCTK_REAL gyzL = INITVALUE; + // CCTK_REAL gzzL = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + // CCTK_REAL kxxL = INITVALUE; + // CCTK_REAL kxyL = INITVALUE; + // CCTK_REAL kxzL = INITVALUE; + // CCTK_REAL kyyL = INITVALUE; + // CCTK_REAL kyzL = INITVALUE; + // CCTK_REAL kzzL = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ /* Assign local copies of grid functions */ - alpL = alp[index]; - betaxL = betax[index]; - betayL = betay[index]; - betazL = betaz[index]; - gxxL = gxx[index]; - gxyL = gxy[index]; - gxzL = gxz[index]; - gyyL = gyy[index]; - gyzL = gyz[index]; - gzzL = gzz[index]; - kxxL = kxx[index]; - kxyL = kxy[index]; - kxzL = kxz[index]; - kyyL = kyy[index]; - kyzL = kyz[index]; - kzzL = kzz[index]; + CCTK_REAL const alpL = alp[index]; + CCTK_REAL const betaxL = betax[index]; + CCTK_REAL const betayL = betay[index]; + CCTK_REAL const betazL = betaz[index]; + CCTK_REAL const gxxL = gxx[index]; + CCTK_REAL const gxyL = gxy[index]; + CCTK_REAL const gxzL = gxz[index]; + CCTK_REAL const gyyL = gyy[index]; + CCTK_REAL const gyzL = gyz[index]; + CCTK_REAL const gzzL = gzz[index]; + CCTK_REAL const kxxL = kxx[index]; + CCTK_REAL const kxyL = kxy[index]; + CCTK_REAL const kxzL = kxz[index]; + CCTK_REAL const kyyL = kyy[index]; + CCTK_REAL const kyzL = kyz[index]; + CCTK_REAL const kzzL = kzz[index]; /* Assign local copies of subblock grid functions */ @@ -148,37 +146,37 @@ void ML_ADM_convertFromADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - g11L = gxxL; + CCTK_REAL const g11L = gxxL; - g12L = gxyL; + CCTK_REAL const g12L = gxyL; - g13L = gxzL; + CCTK_REAL const g13L = gxzL; - g22L = gyyL; + CCTK_REAL const g22L = gyyL; - g23L = gyzL; + CCTK_REAL const g23L = gyzL; - g33L = gzzL; + CCTK_REAL const g33L = gzzL; - K11L = kxxL; + CCTK_REAL const K11L = kxxL; - K12L = kxyL; + CCTK_REAL const K12L = kxyL; - K13L = kxzL; + CCTK_REAL const K13L = kxzL; - K22L = kyyL; + CCTK_REAL const K22L = kyyL; - K23L = kyzL; + CCTK_REAL const K23L = kyzL; - K33L = kzzL; + CCTK_REAL const K33L = kzzL; - alphaL = alpL; + CCTK_REAL const alphaL = alpL; - beta1L = betaxL; + CCTK_REAL const beta1L = betaxL; - beta2L = betayL; + CCTK_REAL const beta2L = betayL; - beta3L = betazL; + CCTK_REAL const beta3L = betazL; /* Copy local copies back to grid functions */ diff --git a/ML_ADM/src/ML_ADM_convertToADMBase.c b/ML_ADM/src/ML_ADM_convertToADMBase.c index 7aadbff..8c9ea1c 100644 --- a/ML_ADM/src/ML_ADM_convertToADMBase.c +++ b/ML_ADM/src/ML_ADM_convertToADMBase.c @@ -20,29 +20,24 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_ADM_convertToADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], CCTK_INT const min[3], CCTK_INT const max[3], CCTK_INT const n_subblock_gfs, CCTK_REAL * const subblock_gfs[]) +void ML_ADM_convertToADMBase_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; + // CCTK_REAL p1o12dx = INITVALUE; + // CCTK_REAL p1o12dy = INITVALUE; + // CCTK_REAL p1o12dz = INITVALUE; + // CCTK_REAL p1o144dxdy = INITVALUE; + // CCTK_REAL p1o144dxdz = INITVALUE; + // CCTK_REAL p1o144dydz = INITVALUE; + // CCTK_REAL pm1o12dx2 = INITVALUE; + // CCTK_REAL pm1o12dy2 = INITVALUE; + // CCTK_REAL pm1o12dz2 = INITVALUE; if (verbose > 1) { @@ -57,31 +52,34 @@ void ML_ADM_convertToADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir, /* Include user-supplied include files */ /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; + CCTK_REAL const dx = CCTK_DELTA_SPACE(0); + CCTK_REAL const dy = CCTK_DELTA_SPACE(1); + CCTK_REAL const dz = CCTK_DELTA_SPACE(2); + int const di = 1; + int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + CCTK_REAL const dxi = 1.0 / dx; + CCTK_REAL const dyi = 1.0 / dy; + CCTK_REAL const dzi = 1.0 / dz; + CCTK_REAL const khalf = 0.5; + CCTK_REAL const kthird = 1/3.0; + CCTK_REAL const ktwothird = 2.0/3.0; + CCTK_REAL const kfourthird = 4.0/3.0; + CCTK_REAL const keightthird = 8.0/3.0; + CCTK_REAL const hdxi = 0.5 * dxi; + CCTK_REAL const hdyi = 0.5 * dyi; + CCTK_REAL const hdzi = 0.5 * dzi; /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; + CCTK_REAL const p1o12dx = INV(dx)/12.; + CCTK_REAL const p1o12dy = INV(dy)/12.; + CCTK_REAL const p1o12dz = INV(dz)/12.; + CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; + CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; + CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; + CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; + CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; + CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ #pragma omp parallel @@ -89,59 +87,59 @@ void ML_ADM_convertToADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + // int index = INITVALUE; + // int subblock_index = INITVALUE; + int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); + int const subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ /* Declare local copies of grid functions */ - CCTK_REAL alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL dtalpL = INITVALUE; - CCTK_REAL dtbetaxL = INITVALUE; - CCTK_REAL dtbetayL = INITVALUE; - CCTK_REAL dtbetazL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; + // CCTK_REAL alpL = INITVALUE; + // CCTK_REAL alphaL = INITVALUE; + // CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; + // CCTK_REAL betaxL = INITVALUE; + // CCTK_REAL betayL = INITVALUE; + // CCTK_REAL betazL = INITVALUE; + // CCTK_REAL dtalpL = INITVALUE; + // CCTK_REAL dtbetaxL = INITVALUE; + // CCTK_REAL dtbetayL = INITVALUE; + // CCTK_REAL dtbetazL = INITVALUE; + // CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + // CCTK_REAL gxxL = INITVALUE; + // CCTK_REAL gxyL = INITVALUE; + // CCTK_REAL gxzL = INITVALUE; + // CCTK_REAL gyyL = INITVALUE; + // CCTK_REAL gyzL = INITVALUE; + // CCTK_REAL gzzL = INITVALUE; + // CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + // CCTK_REAL kxxL = INITVALUE; + // CCTK_REAL kxyL = INITVALUE; + // CCTK_REAL kxzL = INITVALUE; + // CCTK_REAL kyyL = INITVALUE; + // CCTK_REAL kyzL = INITVALUE; + // CCTK_REAL kzzL = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ /* Assign local copies of grid functions */ - alphaL = alpha[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - g11L = g11[index]; - g12L = g12[index]; - g13L = g13[index]; - g22L = g22[index]; - g23L = g23[index]; - g33L = g33[index]; - K11L = K11[index]; - K12L = K12[index]; - K13L = K13[index]; - K22L = K22[index]; - K23L = K23[index]; - K33L = K33[index]; + CCTK_REAL const alphaL = alpha[index]; + CCTK_REAL const beta1L = beta1[index]; + CCTK_REAL const beta2L = beta2[index]; + CCTK_REAL const beta3L = beta3[index]; + CCTK_REAL const g11L = g11[index]; + CCTK_REAL const g12L = g12[index]; + CCTK_REAL const g13L = g13[index]; + CCTK_REAL const g22L = g22[index]; + CCTK_REAL const g23L = g23[index]; + CCTK_REAL const g33L = g33[index]; + CCTK_REAL const K11L = K11[index]; + CCTK_REAL const K12L = K12[index]; + CCTK_REAL const K13L = K13[index]; + CCTK_REAL const K22L = K22[index]; + CCTK_REAL const K23L = K23[index]; + CCTK_REAL const K33L = K33[index]; /* Assign local copies of subblock grid functions */ @@ -152,45 +150,45 @@ void ML_ADM_convertToADMBase_Body(cGH const * const cctkGH, CCTK_INT const dir, /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - gxxL = g11L; + CCTK_REAL const gxxL = g11L; - gxyL = g12L; + CCTK_REAL const gxyL = g12L; - gxzL = g13L; + CCTK_REAL const gxzL = g13L; - gyyL = g22L; + CCTK_REAL const gyyL = g22L; - gyzL = g23L; + CCTK_REAL const gyzL = g23L; - gzzL = g33L; + CCTK_REAL const gzzL = g33L; - kxxL = K11L; + CCTK_REAL const kxxL = K11L; - kxyL = K12L; + CCTK_REAL const kxyL = K12L; - kxzL = K13L; + CCTK_REAL const kxzL = K13L; - kyyL = K22L; + CCTK_REAL const kyyL = K22L; - kyzL = K23L; + CCTK_REAL const kyzL = K23L; - kzzL = K33L; + CCTK_REAL const kzzL = K33L; - alpL = alphaL; + CCTK_REAL const alpL = alphaL; - dtalpL = 0; + CCTK_REAL const dtalpL = 0; - betaxL = beta1L; + CCTK_REAL const betaxL = beta1L; - betayL = beta2L; + CCTK_REAL const betayL = beta2L; - betazL = beta3L; + CCTK_REAL const betazL = beta3L; - dtbetaxL = 0; + CCTK_REAL const dtbetaxL = 0; - dtbetayL = 0; + CCTK_REAL const dtbetayL = 0; - dtbetazL = 0; + CCTK_REAL const dtbetazL = 0; /* Copy local copies back to grid functions */ |