The carpet patch timed Sat Mar 4 18:38:10 CET 2006, not allowing any longer the use of CCTK_DELTA_SPACE in global mode, had broken our way of computing volume integrals. Here is the upgrade.

Further changes are the transformation of divisions into moltiplications and the boundary corrections in the computation of the integral when PUGH is used. This last part is not tested (but will be), anyway who uses PUGH nowadays?


git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/ADMMass/trunk@9 54511f98-0e4f-0410-826e-eb8b393f5a1e

1 files changed, 37 insertions, 35 deletions

diff --git a/src/surface_integral.c b/src/surface_integral.c
index fad504b..b8a3998 100644
--- a/src/surface_integral.c
+++ b/src/surface_integral.c
@@ -43,9 +43,9 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk, ierr, ghost;
+    CCTK_INT i,j,k, ijk, ierr, ghost, ti, tj, tk, tl;
     CCTK_INT  x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k;
-    CCTK_REAL detg, idetg, oneDX, oneDY, oneDZ, twoDX, twoDY, twoDZ, ds[3];
+    CCTK_REAL detg, idetg, ds[3];
     CCTK_REAL u[3][3], dg[3][3][3];
 
     CCTK_REAL physical_min[3];
@@ -58,18 +58,22 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
 
     char coordinate_parameter_type[9];
 
+#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
+
     /* grid-function strides for ADMMacros */
     const CCTK_INT di = 1;
     const CCTK_INT dj = cctk_lsh[0];
     const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];
 
     /* deonminators for derivatives */
-    twoDX = 2.0 * CCTK_DELTA_SPACE(0);
-    twoDY = 2.0 * CCTK_DELTA_SPACE(1);
-    twoDZ = 2.0 * CCTK_DELTA_SPACE(2);
-    oneDX =  CCTK_DELTA_SPACE(0);
-    oneDY =  CCTK_DELTA_SPACE(1);
-    oneDZ =  CCTK_DELTA_SPACE(2);
+    const CCTK_REAL OneOverTwoDX = 1.0 / (2.0 * CCTK_DELTA_SPACE(0));
+    const CCTK_REAL OneOverTwoDY = 1.0 / (2.0 * CCTK_DELTA_SPACE(1));
+    const CCTK_REAL OneOverTwoDZ = 1.0 / (2.0 * CCTK_DELTA_SPACE(2));
+
+    const CCTK_REAL oneDX =  CCTK_DELTA_SPACE(0);
+    const CCTK_REAL oneDY =  CCTK_DELTA_SPACE(1);
+    const CCTK_REAL oneDZ =  CCTK_DELTA_SPACE(2);
+
 
     x_min_i = -INT_MAX;
     x_max_i =  INT_MAX; 
@@ -156,8 +160,6 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     if ( (z[0] < *ADMMass_box_z_max) && (*ADMMass_box_z_max < z[ijk]) )
       z_max_k = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, z, *ADMMass_box_z_max, 2);
 
-#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
-
     for(i=ghost; i<cctk_lsh[0]-ghost; i++)
      for(j=ghost; j<cctk_lsh[1]-ghost; j++)
       for(k=ghost; k<cctk_lsh[2]-ghost; k++)
@@ -216,29 +218,29 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
                   u[2][1] = UPPERMET_UYZ;
                   u[2][2] = UPPERMET_UZZ;
 
-		  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) / twoDX;
-		  dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) / twoDY;
-		  dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) / twoDZ;
+		  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) / twoDX;
-		  dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) / twoDY;
-		  dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) / twoDZ;
+		  dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) / twoDX;
-		  dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) / twoDY;
-		  dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) / twoDZ;
+		  dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) * OneOverTwoDZ;
 		  
 		  dg[1][0][0] = dg[0][1][0];
 		  dg[1][0][1] = dg[0][1][1];
 		  dg[1][0][2] = dg[0][1][2];
 		  
-		  dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) / twoDX;
-		  dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) / twoDY;
-		  dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) / twoDZ;
+		  dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) * OneOverTwoDX;
+		  dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) * OneOverTwoDY;
+		  dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) / twoDX;
-		  dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) / twoDY;
-		  dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) / twoDZ;
+		  dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) * OneOverTwoDX;
+		  dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) * OneOverTwoDZ;
 		  
 		  dg[2][0][0] = dg[0][2][0];
 		  dg[2][0][1] = dg[0][2][1];
@@ -248,14 +250,14 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
 		  dg[2][1][1] = dg[1][2][1];
 		  dg[2][1][2] = dg[1][2][2];
 		  
-		  dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) / twoDX;
-		  dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) / twoDY;
-		  dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) / twoDZ;
-
-                  for (int ti = 0; ti < 3; ti++)
-                   for (int tj = 0; tj < 3; tj++)
-                    for (int tk = 0; tk < 3; tk++)
-                     for (int tl = 0; tl < 3; tl++)
+		  dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) * OneOverTwoDX;
+		  dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) * OneOverTwoDZ;
+
+                  for (ti = 0; ti < 3; ti++)
+                   for (tj = 0; tj < 3; tj++)
+                    for (tk = 0; tk < 3; tk++)
+                     for (tl = 0; tl < 3; tl++)
 		       {
                          ADMMass_SurfaceMass_GF[ijk] +=
 			   u[ti][tj] * u[tk][tl] * ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
@@ -288,7 +290,7 @@ void ADMMass_Surface_Global(CCTK_ARGUMENTS)
     ADMMass_SurfaceMass[*ADMMass_LoopCounter] /=  16.0*PI;
 
     CCTK_VInfo(CCTK_THORNSTRING,
-	       "detector %d: ADM mass  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       "detector %d: ADM mass  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g",
 	       *ADMMass_LoopCounter,
 	       ADMMass_SurfaceMass[*ADMMass_LoopCounter],
 	       *ADMMass_box_x_min,*ADMMass_box_x_max,
@@ -345,7 +347,7 @@ void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS)
     ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter] /=  16.0*PI;    
 
     CCTK_VInfo(CCTK_THORNSTRING,
-	       "detector %d: ADM mass with lapse  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       "detector %d: ADM mass with lapse  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g",
 	       *ADMMass_LoopCounter,
 	       ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter],
 	       *ADMMass_box_x_min,*ADMMass_box_x_max,