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

1 files changed, 178 insertions, 0 deletions

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