- make results independent from the number of processors (not yet finished)

- change some parameter names, so that, to me, they are more clear (debatable)

- add another way of setting the integration surface position

- move constant divisions out of the loops

- modify info messages


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

6 files changed, 238 insertions, 99 deletions

diff --git a/README b/README
index 9c4e62b..883b9c8 100644
--- a/README
+++ b/README
@@ -4,7 +4,7 @@ This thorn can be used to calculate different approximations of the ADM mass
 on a finite grid.
 
 Warning: This thorn is still _very_ experimental. Do not expect it to work,
-         especially not the way you might think it is.
+         especially not the way you might think it should.
 
 Copyright
 ---------
diff --git a/interface.ccl b/interface.ccl
index bb2660d..05bed81 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -6,19 +6,30 @@ inherits: ADMBase ADMMacros StaticConformal SpaceMask
 
 USES INCLUDE: SpaceMask.h
 
+CCTK_INT FUNCTION GetDomainSpecification \
+  (CCTK_INT IN size, \
+   CCTK_REAL OUT ARRAY physical_min, \
+   CCTK_REAL OUT ARRAY physical_max, \
+   CCTK_REAL OUT ARRAY interior_min, \
+   CCTK_REAL OUT ARRAY interior_max, \
+   CCTK_REAL OUT ARRAY exterior_min, \
+   CCTK_REAL OUT ARRAY exterior_max, \
+   CCTK_REAL OUT ARRAY spacing)
+REQUIRES FUNCTION GetDomainSpecification
+
 CCTK_INT ADMMass_LoopCounterG type = SCALAR
 {
   ADMMass_LoopCounter
 } "ADMMass LoopCounter"
 
-CCTK_REAL ADMMass_Masses[ADMMass_number] type = SCALAR
+CCTK_REAL ADMMass_Masses[ADMMass_number] type = SCALAR tags='checkpoint="no"'
 {
   ADMMass_SurfaceMass
   ADMMass_SurfaceMass_Lapse
   ADMMass_VolumeMass
 } "ADMMass Scalars"
 
-CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"'
+CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar" checkpoint="no"'
 {
   ADMMass_SurfaceMass_GF
   ADMMass_VolumeMass_pot_x
@@ -26,3 +37,13 @@ CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensort
   ADMMass_VolumeMass_pot_z
   ADMMass_VolumeMass_GF
 } "ADMMass gridfunctions"
+
+CCTK_REAL  ADMMass_box type = scalar tags='checkpoint="no"'
+{
+ ADMMass_box_x_min
+ ADMMass_box_x_max
+ ADMMass_box_y_min
+ ADMMass_box_y_max
+ ADMMass_box_z_min
+ ADMMass_box_z_max
+} "Physical coordinates of the surface on which the integral is computed"
diff --git a/param.ccl b/param.ccl
index 51a01e2..f2337a7 100644
--- a/param.ccl
+++ b/param.ccl
@@ -6,59 +6,69 @@ CCTK_INT ADMMass_number "number of measurements" STEERABLE=ALWAYS
   0: :: "0 or positive"
 } 1
 
-CCTK_REAL ADMMass_x_min[100] "x_min" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_x_min[100] "x position of the leftmost yz-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } -100.0
-CCTK_REAL ADMMass_x_max[100] "x_max" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_x_max[100] "x position of the righttmost yz-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 100.0
 
-CCTK_REAL ADMMass_y_min[100] "y_min" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_y_min[100] "y position of the leftmost xz-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } -100.0
-CCTK_REAL ADMMass_y_max[100] "y_max" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_y_max[100] "y position of the rightmost xz-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 100.0
 
-CCTK_REAL ADMMass_z_min[100] "z_min" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_z_min[100] "z position of the leftmost xy-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } -100.0
-CCTK_REAL ADMMass_z_max[100] "z_max" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_z_max[100] "z position of the rightmost xy-plane for the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 100.0
 
-CCTK_REAL ADMMass_x_pos[100] "x_pos" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_x_pos[100] "x position of the center of the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 0.0
-CCTK_REAL ADMMass_y_pos[100] "y_pos" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_y_pos[100] "y position of the center of the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 0.0
-CCTK_REAL ADMMass_z_pos[100] "z_pos" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_z_pos[100] "z position of the center of the integration box" STEERABLE=ALWAYS
 {
   :	:: "anything"
 } 0.0
 
-CCTK_REAL ADMMass_boundary_distance[100] "distance to boundary" STEERABLE=ALWAYS
+CCTK_REAL ADMMass_distance_from_grid_boundary[100] "distance between the grid boundaries and the surface of the integration box" STEERABLE=ALWAYS
 {
   :	:: "<=0 for disable, positive otherwise"
 } -1.0
-CCTK_REAL ADMMass_volume_radius[100] "radius of volume extraction" STEERABLE=ALWAYS
+
+CCTK_REAL ADMMass_surface_distance[100] "distance between the above-defined center of the integration (cubic) box and its surface" STEERABLE=ALWAYS
+{
+  :	:: "<=0 for disable, positive otherwise"
+} -1.0
+
+CCTK_REAL ADMMass_volume_radius[100] "radius of the sphere inside which the volume integral is computed" STEERABLE=ALWAYS
 {
   :	:: "<=0 for disable, positive otherwise"
 } -1.0
 
-BOOLEAN ADMMass_use_boundary_distance_as_volume_radius "Use ADMMass_boundary_distance insteat of ADMMass_volume_radius" STEERABLE=ALWAYS
+BOOLEAN ADMMass_use_surface_distance_as_volume_radius "Use ADMMass_surface_distance instead of ADMMass_volume_radius" STEERABLE=ALWAYS
 {
 } "yes"
 
-BOOLEAN ADMMass_Excise_Horizons "Do not calculate the volume integral inside AHs" STEERABLE=ALWAYS
+BOOLEAN ADMMass_use_all_volume_as_volume_radius "Use ADMMass_distance_from_grid_boundary instead of ADMMass_volume_radius" STEERABLE=ALWAYS
+{
+} "no"
+
+BOOLEAN ADMMass_Excise_Horizons "Should we calculate the volume integral inside AHs?" STEERABLE=ALWAYS
 {
 } "no"
diff --git a/schedule.ccl b/schedule.ccl
index 9196694..437ebeb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -13,6 +13,7 @@ schedule ADMMass_InitLoopCounter AT POSTSTEP AFTER OutsideMask_UpdateMask
 
 schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter
 {
+  STORAGE:ADMMass_box
 } "ADMMass loop"
 
 schedule ADMMass_Loop IN ADMMass
diff --git a/src/surface_integral.c b/src/surface_integral.c
index d84347d..7927a06 100644
--- a/src/surface_integral.c
+++ b/src/surface_integral.c
@@ -1,20 +1,27 @@
 #include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
 
 #include <cctk.h>
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
-CCTK_INT find_closest(const cGH *cctkGH, CCTK_INT *cctk_lsh,
-                      CCTK_REAL *x, CCTK_REAL x_min, CCTK_INT dir)
+#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
+
+CCTK_INT find_closest(const cGH *cctkGH, CCTK_INT *cctk_lsh, CCTK_REAL *cctk_delta_space,
+                      CCTK_INT ghost, CCTK_REAL *coord, CCTK_REAL coord_min, CCTK_INT dir)
 {
     CCTK_INT i, ijk, min_i = -1;
     CCTK_REAL min = 1.e100;
-    for(i=3; i<cctk_lsh[dir]-3; i++)
+    
+    for(i=ghost; i<cctk_lsh[dir]-ghost; i++)
     {
         ijk = CCTK_GFINDEX3D(cctkGH, (dir==0)?i:0, (dir==1)?i:0, (dir==2)?i:0);
-        if (fabs(x[ijk] - x_min) < min)
+
+        if (fabs(coord[ijk] - coord_min) < min)
         {
-            min = fabs(x[ijk] - x_min);
+            min = fabs(coord[ijk] - coord_min);
             min_i = i;
         }
     }
@@ -26,60 +33,122 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk;
-    CCTK_REAL x_min, x_max, y_min, y_max, z_min, z_max;
-    CCTK_INT  x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k, min;
+    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 u[3][3], dg[3][3][3];
-    CCTK_REAL local_sum = 0.0;
+
+    CCTK_REAL physical_min[3];
+    CCTK_REAL physical_max[3];
+    CCTK_REAL interior_min[3];
+    CCTK_REAL interior_max[3];
+    CCTK_REAL exterior_min[3];
+    CCTK_REAL exterior_max[3];
+    CCTK_REAL spacing[3];
+
+    char coordinate_parameter_type[9];
 
     /* 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];
+    
+    x_min_i = -10000000;
+    x_max_i =  10000000; 
+    y_min_j = -10000000; 
+    y_max_j =  10000000;
+    z_min_k = -10000000;
+    z_max_k =  10000000;
 
-    if (ADMMass_boundary_distance[*ADMMass_LoopCounter] > 0.0)
-    {
-        x_min = ADMMass_x_pos[*ADMMass_LoopCounter] -
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-        x_max = ADMMass_x_pos[*ADMMass_LoopCounter] +
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-        y_min = ADMMass_y_pos[*ADMMass_LoopCounter] -
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-        y_max = ADMMass_y_pos[*ADMMass_LoopCounter] +
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-        z_min = ADMMass_z_pos[*ADMMass_LoopCounter] -
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-        z_max = ADMMass_z_pos[*ADMMass_LoopCounter] +
-                ADMMass_boundary_distance[*ADMMass_LoopCounter];
-    }
+    if (ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter] > 0.0)
+      { 
+	if (strcmp(CCTK_ParameterValString("type", "cartgrid3d"), "coordbase"))
+	  {	    
+	    CCTK_WARN(0,"this thorns used with the ADMMass_distance_from_grid_boundar parameter requires to set coordinates through coordbase.");
+	    }
+
+	/* Find the physical coordinates of the boundaries */
+	ierr = GetDomainSpecification
+	  (3, physical_min, physical_max, interior_min, interior_max,
+	   exterior_min, exterior_max, spacing);
+	if (ierr)
+	  CCTK_VWarn ( 0, __LINE__, __FILE__, "CartGrid3D", "error returned from function GetDomainSpecification");
+	*ADMMass_box_x_min = physical_min[0] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter];
+  	*ADMMass_box_x_max = physical_max[0] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter];
+	*ADMMass_box_y_min = physical_min[1] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; 
+	*ADMMass_box_y_max = physical_max[1] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter];
+	*ADMMass_box_z_min = physical_min[2] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter];
+	*ADMMass_box_z_max = physical_max[2] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter];
+      }
+    else if (ADMMass_surface_distance[*ADMMass_LoopCounter] > 0.0)
+      {
+        *ADMMass_box_x_min = ADMMass_x_pos[*ADMMass_LoopCounter] -
+	                    ADMMass_surface_distance[*ADMMass_LoopCounter];
+        *ADMMass_box_x_max = ADMMass_x_pos[*ADMMass_LoopCounter] +
+                            ADMMass_surface_distance[*ADMMass_LoopCounter];
+        *ADMMass_box_y_min = ADMMass_y_pos[*ADMMass_LoopCounter] -
+                            ADMMass_surface_distance[*ADMMass_LoopCounter];
+        *ADMMass_box_y_max = ADMMass_y_pos[*ADMMass_LoopCounter] +
+                            ADMMass_surface_distance[*ADMMass_LoopCounter];
+        *ADMMass_box_z_min = ADMMass_z_pos[*ADMMass_LoopCounter] -
+                            ADMMass_surface_distance[*ADMMass_LoopCounter];
+        *ADMMass_box_z_max = ADMMass_z_pos[*ADMMass_LoopCounter] +
+                            ADMMass_surface_distance[*ADMMass_LoopCounter];
+      }
     else
-    {
-        x_min = ADMMass_x_min[*ADMMass_LoopCounter];
-        x_max = ADMMass_x_max[*ADMMass_LoopCounter];
-        y_min = ADMMass_y_min[*ADMMass_LoopCounter];
-        y_max = ADMMass_y_max[*ADMMass_LoopCounter];
-        z_min = ADMMass_z_min[*ADMMass_LoopCounter];
-        z_max = ADMMass_z_max[*ADMMass_LoopCounter];
-    }
-    x_min_i = find_closest(cctkGH, cctk_lsh, x, x_min, 0);
-    x_max_i = find_closest(cctkGH, cctk_lsh, x, x_max, 0);
-    y_min_j = find_closest(cctkGH, cctk_lsh, y, y_min, 1);
-    y_max_j = find_closest(cctkGH, cctk_lsh, y, y_max, 1);
-    z_min_k = find_closest(cctkGH, cctk_lsh, z, z_min, 2);
-    z_max_k = find_closest(cctkGH, cctk_lsh, z, z_max, 2);
+      {
+        *ADMMass_box_x_min = ADMMass_x_min[*ADMMass_LoopCounter];
+        *ADMMass_box_x_max = ADMMass_x_max[*ADMMass_LoopCounter];
+        *ADMMass_box_y_min = ADMMass_y_min[*ADMMass_LoopCounter];
+        *ADMMass_box_y_max = ADMMass_y_max[*ADMMass_LoopCounter];
+        *ADMMass_box_z_min = ADMMass_z_min[*ADMMass_LoopCounter];
+        *ADMMass_box_z_max = ADMMass_z_max[*ADMMass_LoopCounter];
+      }
+
+    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;
+      }
+
+    ijk = CCTK_GFINDEX3D(cctkGH, cctk_ubnd[0]-cctk_lbnd[0], cctk_ubnd[1]-cctk_lbnd[1],cctk_ubnd[2]-cctk_lbnd[2] );
+
+    if ( (x[0] < *ADMMass_box_x_min) && (*ADMMass_box_x_min < x[ijk]) )
+      x_min_i = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, x, *ADMMass_box_x_min, 0);
+
+    if ( (x[0] < *ADMMass_box_x_max) && (*ADMMass_box_x_max < x[ijk]) )
+      x_max_i = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, x, *ADMMass_box_x_max, 0);
+
+    if ( (y[0] < *ADMMass_box_y_min) && (*ADMMass_box_y_min < y[ijk]) )
+      y_min_j = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, y, *ADMMass_box_y_min, 1);
+
+    if ( (y[0] < *ADMMass_box_y_max) && (*ADMMass_box_y_max < y[ijk]) )
+      y_max_j = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, y, *ADMMass_box_y_max, 1);
+
+    if ( (z[0] < *ADMMass_box_z_min) && (*ADMMass_box_z_min < z[ijk]) )
+      z_min_k = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, z, *ADMMass_box_z_min, 2);
+
+    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"
 #include "CactusEinstein/ADMMacros/src/macro/DG_declare.h"
 #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h"
 
-    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++)
+    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] = 0.0;
-          if ((i >= x_min_i) &&
+	  if ((i >= x_min_i) &&
               (i <= x_max_i) &&
               (j >= y_min_j) &&
               (j <= y_max_j) &&
@@ -149,11 +218,12 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
                          ADMMass_SurfaceMass_GF[ijk] +=
                              u[ti][tj] * u[tk][tl] *
                              ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
-                  ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG) / 16 / 3.14159265;
-                  local_sum += ADMMass_SurfaceMass_GF[ijk];
+                  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"
@@ -171,18 +241,23 @@ void ADMMass_Surface_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,
                     &ADMMass_SurfaceMass[*ADMMass_LoopCounter], 1,
                     CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF")))
         CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF");
-    
+
+    ADMMass_SurfaceMass[*ADMMass_LoopCounter] /=  16.0*PI;
+
     CCTK_VInfo(CCTK_THORNSTRING,
-               "ADM mass, surface, detector %d: %g %g\n", *ADMMass_LoopCounter,
-               ADMMass_boundary_distance[*ADMMass_LoopCounter],
-               ADMMass_SurfaceMass[*ADMMass_LoopCounter]);
+	       "detector %d: ADM mass  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       *ADMMass_LoopCounter,
+	       ADMMass_SurfaceMass[*ADMMass_LoopCounter],
+	       *ADMMass_box_x_min,*ADMMass_box_x_max,
+	       *ADMMass_box_y_min,*ADMMass_box_y_max,
+	       *ADMMass_box_z_min,*ADMMass_box_z_max);
 }
 
 void ADMMass_Surface_Lapse(CCTK_ARGUMENTS)
@@ -191,15 +266,13 @@ void ADMMass_Surface_Lapse(CCTK_ARGUMENTS)
     DECLARE_CCTK_PARAMETERS
 
     CCTK_INT i,j,k, ijk;
-    CCTK_REAL local_sum = 0.0;
-
+    
     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++)
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
           ADMMass_SurfaceMass_GF[ijk] *= alp[ijk];
-          local_sum += ADMMass_SurfaceMass_GF[ijk];
       }
 }
 
@@ -219,11 +292,15 @@ void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS)
                     &ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter], 1,
                     CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF")))
         CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF");
-    
+
+    ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter] /=  16.0*PI;    
+
     CCTK_VInfo(CCTK_THORNSTRING,
-           "ADM mass, surface lapse, detector %d: %g %g\n",
-           *ADMMass_LoopCounter,
-           ADMMass_boundary_distance[*ADMMass_LoopCounter],
-           ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter]);
+	       "detector %d: ADM mass with lapse  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       *ADMMass_LoopCounter,
+	       ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter],
+	       *ADMMass_box_x_min,*ADMMass_box_x_max,
+	       *ADMMass_box_y_min,*ADMMass_box_y_max,
+	       *ADMMass_box_z_min,*ADMMass_box_z_max);
 }
 
diff --git a/src/volume_integral.c b/src/volume_integral.c
index 7db9d6d..b80b215 100644
--- a/src/volume_integral.c
+++ b/src/volume_integral.c
@@ -1,7 +1,7 @@
-
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <assert.h>
 
 #include <cctk.h>
 #include <cctk_Arguments.h>
@@ -16,10 +16,10 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk;
+    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 local_sum = 0.0;*/
     CCTK_REAL radius;
 
     CCTK_INT mask_descriptor, state_descriptor_outside;
@@ -32,26 +32,40 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
     const CCTK_INT dj = cctk_lsh[0];
     const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];
 
-    if (ADMMass_use_boundary_distance_as_volume_radius &&
+    if (ADMMass_use_surface_distance_as_volume_radius &&
         (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0))
-        radius = ADMMass_boundary_distance[*ADMMass_LoopCounter];
+        radius = ADMMass_surface_distance[*ADMMass_LoopCounter];
     else
         radius = ADMMass_volume_radius[*ADMMass_LoopCounter];
 
-    if (radius <= 0.0)
+    if ((radius <= 0.0)&&(!ADMMass_use_all_volume_as_volume_radius))
     {
         CCTK_WARN(1, "radius < 0 / not set, not calculating "
                      "the volume integral to get the ADM mass.");
         return;
     }
 
+    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;
+      }
+
 #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=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++)
+    /* 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++)   
       {
           ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
           ADMMass_VolumeMass_GF[ijk] = 0.0;
@@ -111,20 +125,19 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
                     u[ti][tj] * u[tk][2] *
                     ( dg[ti][tk][tj] - dg[ti][tj][tk] );
             }
-          ADMMass_VolumeMass_pot_x[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
-                                           16 / PI;
-          ADMMass_VolumeMass_pot_y[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
-                                           16 / PI;
-          ADMMass_VolumeMass_pot_z[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
-                                           16 / PI;
+          ADMMass_VolumeMass_pot_x[ijk] *= alp[ijk] * sqrt(DETG_DETG);
+          ADMMass_VolumeMass_pot_y[ijk] *= alp[ijk] * sqrt(DETG_DETG);
+          ADMMass_VolumeMass_pot_z[ijk] *= alp[ijk] * sqrt(DETG_DETG);
       }
-    for(i=4; i<cctk_lsh[0]-4; i++)
-     for(j=4; j<cctk_lsh[1]-4; j++)
-      for(k=4; k<cctk_lsh[2]-4; k++)
+
+    /* Do not compute in ghost zones */
+    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_VolumeMass_GF[ijk] = 0.0;
-          if ((!ADMMass_Excise_Horizons ||
+          if ((ADMMass_use_all_volume_as_volume_radius)||((!ADMMass_Excise_Horizons ||
                SpaceMask_CheckStateBits(space_mask, ijk,
                                         mask_descriptor,
                                         state_descriptor_outside)) &&
@@ -134,7 +147,7 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
                (y[ijk]-ADMMass_y_pos[*ADMMass_LoopCounter]) +
                (z[ijk]-ADMMass_z_pos[*ADMMass_LoopCounter])*
                (z[ijk]-ADMMass_z_pos[*ADMMass_LoopCounter]) <=
-               radius * radius))
+               radius * radius)))
           {
               ADMMass_VolumeMass_GF[ijk] +=
                (i2dx*(ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]-
@@ -144,7 +157,7 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
                 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];
+              /* local_sum += ADMMass_VolumeMass_GF[ijk]; */
           }
       }
     #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
@@ -171,11 +184,28 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS)
                     CCTK_VarIndex("ADMMass::ADMMass_VolumeMass_GF")))
         CCTK_WARN(0, "Error while reducing ADMMass_VolumeMass_GF");
     
-    if (ADMMass_use_boundary_distance_as_volume_radius)
-        radius = ADMMass_boundary_distance[*ADMMass_LoopCounter];
+    ADMMass_VolumeMass[*ADMMass_LoopCounter] /=  16.0*PI;
+
+    if (ADMMass_use_surface_distance_as_volume_radius)
+      {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "detector %d: ADM mass  %g, volume (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+		   *ADMMass_LoopCounter,
+		   ADMMass_SurfaceMass[*ADMMass_LoopCounter],
+		   *ADMMass_box_x_min,*ADMMass_box_x_max,
+		   *ADMMass_box_y_min,*ADMMass_box_y_max,
+		   *ADMMass_box_z_min,*ADMMass_box_z_max);
+      }
+    else if (ADMMass_use_all_volume_as_volume_radius)
+      {
+	    CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume: the whole grid\n", 
+	       *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter]);
+      }
     else
-        radius = ADMMass_volume_radius[*ADMMass_LoopCounter];
-    printf("ADM mass, volume, detector %d: %g %g\n", *ADMMass_LoopCounter,
-           radius, ADMMass_VolumeMass[*ADMMass_LoopCounter]);
+      {
+	CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume (sphere of radius) %g\n", 
+	       *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter],
+	       ADMMass_volume_radius[*ADMMass_LoopCounter]);
+      }
 }