- 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

1 files changed, 56 insertions, 26 deletions

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]);
+      }
 }