Do not use ADMMAcros (C version), since they are bugged (for Carpet). Finish changes that ensure compatibility with Carpet and multiprocessor runs. Not tested extensively.

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

3 files changed, 205 insertions, 101 deletions

diff --git a/src/loop.c b/src/loop.c
index 7cc36bb..34122a0 100644
--- a/src/loop.c
+++ b/src/loop.c
@@ -2,13 +2,23 @@
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
-/* Initialises the loop counter to the last radius/distance set */
+
+/* Initialise the loop counter */
 void ADMMass_InitLoopCounter(CCTK_ARGUMENTS)
 {
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    *ADMMass_LoopCounter = ADMMass_number;
+          *ADMMass_LoopCounter = 0;
+}
+
+/* Set the loop counter to the value of the parameter ADMMass:ADMMass_number */
+void ADMMass_SetLoopCounter(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+      *ADMMass_LoopCounter = ADMMass_number;
 }
 
 /* Decrements the counter to loop over all radii/distances set */
@@ -18,4 +28,5 @@ void ADMMass_Loop(CCTK_ARGUMENTS)
     DECLARE_CCTK_PARAMETERS
 
     (*ADMMass_LoopCounter)--;
+
 }
diff --git a/src/surface_integral.c b/src/surface_integral.c
index 7927a06..fad504b 100644
--- a/src/surface_integral.c
+++ b/src/surface_integral.c
@@ -1,7 +1,17 @@
+ /*@@
+   @file      surface_integral.c
+   @date      Wed Oct 5 10:02:12 2005
+   @author    Frank Loeffler and Luca Baiotti
+   @desc 
+          Computes the ADM mass as a surface integral.
+   @enddesc 
+ @@*/
+
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <assert.h>
+#include <limits.h>
 
 #include <cctk.h>
 #include <cctk_Arguments.h>
@@ -35,7 +45,7 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
 
     CCTK_INT i,j,k, ijk, ierr, ghost;
     CCTK_INT  x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k;
-    CCTK_REAL detg, idetg, ds[3];
+    CCTK_REAL detg, idetg, oneDX, oneDY, oneDZ, twoDX, twoDY, twoDZ, ds[3];
     CCTK_REAL u[3][3], dg[3][3][3];
 
     CCTK_REAL physical_min[3];
@@ -52,13 +62,21 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     const CCTK_INT di = 1;
     const CCTK_INT dj = cctk_lsh[0];
     const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];
-    
-    x_min_i = -10000000;
-    x_max_i =  10000000; 
-    y_min_j = -10000000; 
-    y_max_j =  10000000;
-    z_min_k = -10000000;
-    z_max_k =  10000000;
+
+    /* 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);
+
+    x_min_i = -INT_MAX;
+    x_max_i =  INT_MAX; 
+    y_min_j = -INT_MAX; 
+    y_max_j =  INT_MAX;
+    z_min_k = -INT_MAX;
+    z_max_k =  INT_MAX;
 
     if (ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter] > 0.0)
       { 
@@ -139,8 +157,6 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
       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"
-#include "CactusEinstein/ADMMacros/src/macro/DG_declare.h"
-#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h"
 
     for(i=ghost; i<cctk_lsh[0]-ghost; i++)
      for(j=ghost; j<cctk_lsh[1]-ghost; j++)
@@ -148,6 +164,8 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
           ADMMass_SurfaceMass_GF[ijk] = 0.0;
+	  
+	  /* delimit the cube on whose surface we want to integrate */
 	  if ((i >= x_min_i) &&
               (i <= x_max_i) &&
               (j >= y_min_j) &&
@@ -158,23 +176,36 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
               ds[0] = 0.0;
               ds[1] = 0.0;
               ds[2] = 0.0;
-#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing.h"
-              if ((i == x_min_i) && (j != y_max_j) && (k != z_max_k))
-                  ds[0] = -dy*dz;
-              if ((i == x_max_i) && (j != y_max_j) && (k != z_max_k))
-                  ds[0] = dy*dz;
-              if ((j == y_min_j) && (i != x_max_i) && (k != z_max_k))
-                  ds[1] = -dx*dz;
-              if ((j == y_max_j) && (i != x_max_i) && (k != z_max_k))
-                  ds[1] = dx*dz;
-              if ((k == z_min_k) && (j != y_max_j) && (i != x_max_i))
-                  ds[2] = -dx*dy;
-              if ((k == z_max_k) && (j != y_max_j) && (i != x_max_i))
-                  ds[2] = dx*dy;
+
+	      /* Select the points on the surfaces of the requested cube, 
+		 paying attention not to compute the corners (and the verteces) more than once. */
+	     
+	      /* the whole face of the cube */
+              if (i == x_min_i)  ds[0] = -oneDY*oneDZ;
+
+	      /* the whole face of the cube */
+	      if (i == x_max_i)  ds[0] = oneDY*oneDZ;
+
+	      /* face, excluding the two corners (and the four verteces) that were already computed above*/
+              if ((j == y_min_j) && (i != x_min_i) && (i != x_max_i))
+                  ds[1] = -oneDX*oneDZ;
+
+	      /* face, excluding the other two corners (and the four verteces) 
+		 that were already computed above*/
+              if ((j == y_max_j) && (i != x_min_i) && (i != x_max_i))
+                  ds[1] = oneDX*oneDZ;
+
+	      /* face, excluding all corners (and verteces), since they were already computed above */
+              if ((k == z_min_k) && (j != y_min_j) && (i != x_min_i) && (j != y_max_j) && (i != x_max_i))
+                  ds[2] = -oneDX*oneDY;
+
+	      /* face, excluding all corners (and verteces), since they were already computed above */
+              if ((k == z_max_k) && (j != y_min_j) && (i != x_min_i) && (j != y_max_j) && (i != x_max_i))
+                  ds[2] = oneDX*oneDY;
+
               if ((ds[0] != 0.0) || (ds[1] != 0.0) || (ds[2] != 0.0))
               {
 #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_guts.h"
-#include "CactusEinstein/ADMMacros/src/macro/DG_guts.h"
                   u[0][0] = UPPERMET_UXX;
                   u[0][1] = UPPERMET_UXY;
                   u[0][2] = UPPERMET_UXZ;
@@ -184,52 +215,57 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
                   u[2][0] = UPPERMET_UXZ;
                   u[2][1] = UPPERMET_UYZ;
                   u[2][2] = UPPERMET_UZZ;
-                  dg[0][0][0] = DXDG_DXDGXX;
-                  dg[0][0][1] = DYDG_DYDGXX;
-                  dg[0][0][2] = DZDG_DZDGXX;
-                  dg[0][1][0] = DXDG_DXDGXY;
-                  dg[0][1][1] = DYDG_DYDGXY;
-                  dg[0][1][2] = DZDG_DZDGXY;
-                  dg[0][2][0] = DXDG_DXDGXZ;
-                  dg[0][2][1] = DYDG_DYDGXZ;
-                  dg[0][2][2] = DZDG_DZDGXZ;
-                  dg[1][0][0] = DXDG_DXDGXY;
-                  dg[1][0][1] = DYDG_DYDGXY;
-                  dg[1][0][2] = DZDG_DZDGXY;
-                  dg[1][1][0] = DXDG_DXDGYY;
-                  dg[1][1][1] = DYDG_DYDGYY;
-                  dg[1][1][2] = DZDG_DZDGYY;
-                  dg[1][2][0] = DXDG_DXDGYZ;
-                  dg[1][2][1] = DYDG_DYDGYZ;
-                  dg[1][2][2] = DZDG_DZDGYZ;
-                  dg[2][0][0] = DXDG_DXDGXZ;
-                  dg[2][0][1] = DYDG_DYDGXZ;
-                  dg[2][0][2] = DZDG_DZDGXZ;
-                  dg[2][1][0] = DXDG_DXDGYZ;
-                  dg[2][1][1] = DYDG_DYDGYZ;
-                  dg[2][1][2] = DZDG_DZDGYZ;
-                  dg[2][2][0] = DXDG_DXDGZZ;
-                  dg[2][2][1] = DYDG_DYDGZZ;
-                  dg[2][2][2] = DZDG_DZDGZZ;
+
+		  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][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][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[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][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[2][0][0] = dg[0][2][0];
+		  dg[2][0][1] = dg[0][2][1];
+		  dg[2][0][2] = dg[0][2][2];
+		  
+		  dg[2][1][0] = dg[1][2][0];
+		  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++)
+		       {
                          ADMMass_SurfaceMass_GF[ijk] +=
-                             u[ti][tj] * u[tk][tl] *
-                             ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
+			   u[ti][tj] * u[tk][tl] * ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
+		       }
                   ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG);
-                  /* local_sum += ADMMass_SurfaceMass_GF[ijk];*/
               }
           }
       }
 
-    #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
-    #include "CactusEinstein/ADMMacros/src/macro/DG_undefine.h"
-    #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_undefine.h"
-    /*CCTK_VInfo(CCTK_THORNSTRING,
-     *           "ADM mass local to one processor (surface): %g\n", local_sum);
-     */
+#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
 }
 
 void ADMMass_Surface_Global(CCTK_ARGUMENTS)
@@ -265,11 +301,24 @@ void ADMMass_Surface_Lapse(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk;
-    
-    for(i=3; i<cctk_lsh[0]-3; i++)
-     for(j=3; j<cctk_lsh[1]-3; j++)
-      for(k=3; k<cctk_lsh[2]-3; k++)
+    CCTK_INT i,j,k, ijk, ghost;
+
+    if (CCTK_IsThornActive("PUGH"))
+      {
+	const void *ghost_ptr = CCTK_ParameterGet("ghost_size","PUGH",NULL);
+	assert( ghost_ptr != NULL );
+	ghost = *(const int *)ghost_ptr;
+      }
+    else /* carpet */ 
+      {
+	const void *ghost_ptr = CCTK_ParameterGet("ghost_size","Carpet",NULL);
+	assert( ghost_ptr != NULL );
+	ghost = *(const int *)ghost_ptr;
+      }
+
+    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++)
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
           ADMMass_SurfaceMass_GF[ijk] *= alp[ijk];
diff --git a/src/volume_integral.c b/src/volume_integral.c
index b80b215..ee4009d 100644
--- a/src/volume_integral.c
+++ b/src/volume_integral.c
@@ -1,3 +1,12 @@
+ /*@@
+   @file      volume_integral.c
+   @date      Wed Oct 5 10:18:34 2005
+   @author    Frank Loeffler and Luca Baiotti
+   @desc 
+          Computes the ADM mass as a volume integral.
+   @enddesc 
+ @@*/
+
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -19,19 +28,23 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
     CCTK_INT i,j,k, ijk, ghost;
     CCTK_REAL detg, idetg;
     CCTK_REAL u[3][3], dg[3][3][3];
-    /* CCTK_REAL local_sum = 0.0;*/
-    CCTK_REAL radius;
+    CCTK_REAL radius, twoDX, twoDY, twoDZ;
 
     CCTK_INT mask_descriptor, state_descriptor_outside;
 
     mask_descriptor = SpaceMask_GetTypeBits("OutsideMask");
     state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside");
 
-    /* grid-function strides for ADMMacros */
+    /* grid-function strides */
     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);
+
     if (ADMMass_use_surface_distance_as_volume_radius &&
         (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0))
         radius = ADMMass_surface_distance[*ADMMass_LoopCounter];
@@ -59,22 +72,18 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
       }
 
 #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
-#include "CactusEinstein/ADMMacros/src/macro/DG_declare.h"
-#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h"
 
-    /* Compute everywhere. These ADMMacros in C only have second order differencing. */
-    for(i=2; i<cctk_lsh[0]-2; i++)
-     for(j=2; j<cctk_lsh[1]-2; j++)
-      for(k=2; k<cctk_lsh[2]-2; k++)   
+    for(i=1; i<cctk_lsh[0]-1; i++)
+     for(j=1; j<cctk_lsh[1]-1; j++)
+      for(k=1; k<cctk_lsh[2]-1; k++)   
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
           ADMMass_VolumeMass_GF[ijk] = 0.0;
           ADMMass_VolumeMass_pot_x[ijk] = 0.0;
           ADMMass_VolumeMass_pot_y[ijk] = 0.0;
           ADMMass_VolumeMass_pot_z[ijk] = 0.0;
-#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing.h"
+
 #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_guts.h"
-#include "CactusEinstein/ADMMacros/src/macro/DG_guts.h"
           u[0][0] = UPPERMET_UXX;
           u[0][1] = UPPERMET_UXY;
           u[0][2] = UPPERMET_UXZ;
@@ -84,7 +93,44 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
           u[2][0] = UPPERMET_UXZ;
           u[2][1] = UPPERMET_UYZ;
           u[2][2] = UPPERMET_UZZ;
-          dg[0][0][0] = DXDG_DXDGXX;
+
+	  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][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][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[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][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[2][0][0] = dg[0][2][0];
+          dg[2][0][1] = dg[0][2][1];
+          dg[2][0][2] = dg[0][2][2];
+
+          dg[2][1][0] = dg[1][2][0];
+          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;
+
+	  /*	  dg[0][0][0] = DXDG_DXDGXX;
           dg[0][0][1] = DYDG_DYDGXX;
           dg[0][0][2] = DZDG_DZDGXX;
           dg[0][1][0] = DXDG_DXDGXY;
@@ -110,20 +156,21 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
           dg[2][1][2] = DZDG_DZDGYZ;
           dg[2][2][0] = DXDG_DXDGZZ;
           dg[2][2][1] = DYDG_DYDGZZ;
-          dg[2][2][2] = DZDG_DZDGZZ;
+          dg[2][2][2] = DZDG_DZDGZZ;*/
+
           for (int ti = 0; ti < 3; ti++)
            for (int tj = 0; tj < 3; tj++)
             for (int tk = 0; tk < 3; tk++)
             {
                 ADMMass_VolumeMass_pot_x[ijk] +=
                     u[ti][tj] * u[tk][0] *
-                    ( dg[ti][tk][tj] - dg[ti][tj][tk] );
+		  ( dg[ti][tk][tj] - dg[ti][tj][tk] );
                 ADMMass_VolumeMass_pot_y[ijk] +=
                     u[ti][tj] * u[tk][1] *
-                    ( dg[ti][tk][tj] - dg[ti][tj][tk] );
+		  ( dg[ti][tk][tj] - dg[ti][tj][tk] );
                 ADMMass_VolumeMass_pot_z[ijk] +=
                     u[ti][tj] * u[tk][2] *
-                    ( dg[ti][tk][tj] - dg[ti][tj][tk] );
+		  ( dg[ti][tk][tj] - dg[ti][tj][tk] );
             }
           ADMMass_VolumeMass_pot_x[ijk] *= alp[ijk] * sqrt(DETG_DETG);
           ADMMass_VolumeMass_pot_y[ijk] *= alp[ijk] * sqrt(DETG_DETG);
@@ -136,8 +183,8 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
       for(k=ghost; k<cctk_lsh[2]-ghost; k++)
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
-          ADMMass_VolumeMass_GF[ijk] = 0.0;
-          if ((ADMMass_use_all_volume_as_volume_radius)||((!ADMMass_Excise_Horizons ||
+
+	  if ((ADMMass_use_all_volume_as_volume_radius)||((!ADMMass_Excise_Horizons ||
                SpaceMask_CheckStateBits(space_mask, ijk,
                                         mask_descriptor,
                                         state_descriptor_outside)) &&
@@ -149,21 +196,17 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
                (z[ijk]-ADMMass_z_pos[*ADMMass_LoopCounter]) <=
                radius * radius)))
           {
-              ADMMass_VolumeMass_GF[ijk] +=
-               (i2dx*(ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]-
-                      ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i-1,j,k)])+
-                i2dy*(ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]-
-                      ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j-1,k)])+
-                i2dz*(ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]-
-                      ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]))
-               *dx*dy*dz;
-              /* local_sum += ADMMass_VolumeMass_GF[ijk]; */
-          }
+              ADMMass_VolumeMass_GF[ijk] =
+               ((ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]-
+                      ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i-1,j,k)])/twoDX+
+                (ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]-
+                      ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j-1,k)])/twoDY+
+                (ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]-
+                      ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)])/twoDZ);
+
+	      }
       }
-    #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
-    #include "CactusEinstein/ADMMacros/src/macro/DG_undefine.h"
-    #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_undefine.h"
-    /*printf("ADM mass local to one processor (volume): %g\n", local_sum);*/
+#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
 }
 
 void ADMMass_Volume_Global(CCTK_ARGUMENTS)
@@ -176,7 +219,7 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS)
 
     reduction_handle = CCTK_ReductionHandle("sum");
     if (reduction_handle < 0)
-        CCTK_WARN(0, "Unable to ge reduction handle.");
+        CCTK_WARN(0, "Unable to get reduction handle.");
 
     if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1,
                     CCTK_VARIABLE_REAL,
@@ -184,7 +227,8 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS)
                     CCTK_VarIndex("ADMMass::ADMMass_VolumeMass_GF")))
         CCTK_WARN(0, "Error while reducing ADMMass_VolumeMass_GF");
     
-    ADMMass_VolumeMass[*ADMMass_LoopCounter] /=  16.0*PI;
+    ADMMass_VolumeMass[*ADMMass_LoopCounter] *=
+      cctk_delta_space[0]*cctk_delta_space[1]*cctk_delta_space[2] / (16.0*PI);
 
     if (ADMMass_use_surface_distance_as_volume_radius)
       {