thorn to comute the ADM mass

_very_ experimental at the moment


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

4 files changed, 436 insertions, 0 deletions

diff --git a/src/loop.c b/src/loop.c
new file mode 100644
index 0000000..7cc36bb
--- /dev/null
+++ b/src/loop.c
@@ -0,0 +1,21 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+/* Initialises the loop counter to the last radius/distance set */
+void ADMMass_InitLoopCounter(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    *ADMMass_LoopCounter = ADMMass_number;
+}
+
+/* Decrements the counter to loop over all radii/distances set */
+void ADMMass_Loop(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    (*ADMMass_LoopCounter)--;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..400543f
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn ADMMass
+# $Header$
+
+# Source files in this directory
+SRCS = surface_integral.c volume_integral.c loop.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/surface_integral.c b/src/surface_integral.c
new file mode 100644
index 0000000..d84347d
--- /dev/null
+++ b/src/surface_integral.c
@@ -0,0 +1,229 @@
+#include <math.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)
+{
+    CCTK_INT i, ijk, min_i = -1;
+    CCTK_REAL min = 1.e100;
+    for(i=3; i<cctk_lsh[dir]-3; i++)
+    {
+        ijk = CCTK_GFINDEX3D(cctkGH, (dir==0)?i:0, (dir==1)?i:0, (dir==2)?i:0);
+        if (fabs(x[ijk] - x_min) < min)
+        {
+            min = fabs(x[ijk] - x_min);
+            min_i = i;
+        }
+    }
+    return min_i;
+}
+
+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_REAL detg, idetg, ds[3];
+    CCTK_REAL u[3][3], dg[3][3][3];
+    CCTK_REAL local_sum = 0.0;
+
+    /* 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];
+
+    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];
+    }
+    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);
+
+#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++)
+      {
+          ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
+          ADMMass_SurfaceMass_GF[ijk] = 0.0;
+          if ((i >= x_min_i) &&
+              (i <= x_max_i) &&
+              (j >= y_min_j) &&
+              (j <= y_max_j) &&
+              (k >= z_min_k) &&
+              (k <= z_max_k))
+          {
+              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;
+              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;
+                  u[1][0] = UPPERMET_UXY;
+                  u[1][1] = UPPERMET_UYY;
+                  u[1][2] = UPPERMET_UYZ;
+                  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;
+                  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];
+                  ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG) / 16 / 3.14159265;
+                  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);
+     */
+}
+
+void ADMMass_Surface_Global(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    CCTK_INT reduction_handle;
+
+    reduction_handle = CCTK_ReductionHandle("sum");
+    if (reduction_handle < 0)
+        CCTK_WARN(0, "Unable to ge 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");
+    
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "ADM mass, surface, detector %d: %g %g\n", *ADMMass_LoopCounter,
+               ADMMass_boundary_distance[*ADMMass_LoopCounter],
+               ADMMass_SurfaceMass[*ADMMass_LoopCounter]);
+}
+
+void ADMMass_Surface_Lapse(CCTK_ARGUMENTS)
+{
+    DECLARE_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];
+      }
+}
+
+void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    CCTK_INT reduction_handle;
+
+    reduction_handle = CCTK_ReductionHandle("sum");
+    if (reduction_handle < 0)
+        CCTK_WARN(0, "Unable to ge reduction handle.");
+
+    if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1,
+                    CCTK_VARIABLE_REAL,
+                    &ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter], 1,
+                    CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF")))
+        CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF");
+    
+    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]);
+}
+
diff --git a/src/volume_integral.c b/src/volume_integral.c
new file mode 100644
index 0000000..0a9a1dd
--- /dev/null
+++ b/src/volume_integral.c
@@ -0,0 +1,178 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "SpaceMask.h"
+
+void ADMMass_Volume(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    CCTK_INT i,j,k, ijk;
+    CCTK_REAL detg, idetg;
+    CCTK_REAL u[3][3], dg[3][3][3];
+    CCTK_REAL local_sum = 0.0;
+    CCTK_REAL radius;
+
+    CCTK_INT mask_descriptor, state_descriptor_outside;
+
+    mask_descriptor = SpaceMask_GetTypeBits("OutsideMask");
+    state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside");
+
+    /* 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];
+
+    if (ADMMass_use_boundary_distance_as_volume_radius &&
+        (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0))
+        radius = ADMMass_boundary_distance[*ADMMass_LoopCounter];
+    else
+        radius = ADMMass_volume_radius[*ADMMass_LoopCounter];
+
+    if (radius <= 0.0)
+    {
+        CCTK_WARN(1, "radius < 0 / not set, not calculating "
+                     "the volume integral to get the ADM mass.");
+        return;
+    }
+
+#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++)
+      {
+          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;
+          u[1][0] = UPPERMET_UXY;
+          u[1][1] = UPPERMET_UYY;
+          u[1][2] = UPPERMET_UYZ;
+          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;
+          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] );
+                ADMMass_VolumeMass_pot_y[ijk] +=
+                    u[ti][tj] * u[tk][1] *
+                    ( 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] );
+            }
+          ADMMass_VolumeMass_pot_x[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
+                                           16 / 3.14159265;
+          ADMMass_VolumeMass_pot_y[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
+                                           16 / 3.14159265;
+          ADMMass_VolumeMass_pot_z[ijk] *= alp[ijk] * sqrt(DETG_DETG) /
+                                           16 / 3.14159265;
+      }
+    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++)
+      {
+          ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
+          ADMMass_VolumeMass_GF[ijk] = 0.0;
+          if ((!ADMMass_Excise_Horizons ||
+               SpaceMask_CheckStateBits(space_mask, ijk,
+                                        mask_descriptor,
+                                        state_descriptor_outside)) &&
+              ((x[ijk]-ADMMass_x_pos[*ADMMass_LoopCounter])*
+               (x[ijk]-ADMMass_x_pos[*ADMMass_LoopCounter]) +
+               (y[ijk]-ADMMass_y_pos[*ADMMass_LoopCounter])*
+               (y[ijk]-ADMMass_y_pos[*ADMMass_LoopCounter]) +
+               (z[ijk]-ADMMass_z_pos[*ADMMass_LoopCounter])*
+               (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];
+          }
+      }
+    #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);*/
+}
+
+void ADMMass_Volume_Global(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    CCTK_INT reduction_handle;
+    CCTK_REAL radius;
+
+    reduction_handle = CCTK_ReductionHandle("sum");
+    if (reduction_handle < 0)
+        CCTK_WARN(0, "Unable to ge reduction handle.");
+
+    if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1,
+                    CCTK_VARIABLE_REAL,
+                    &ADMMass_VolumeMass[*ADMMass_LoopCounter], 1,
+                    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];
+    else
+        radius = ADMMass_volume_radius[*ADMMass_LoopCounter];
+    printf("ADM mass, volume, detector %d: %g %g\n", *ADMMass_LoopCounter,
+           radius, ADMMass_VolumeMass[*ADMMass_LoopCounter]);
+}
+