diff options
author | baiotti <baiotti@54511f98-0e4f-0410-826e-eb8b393f5a1e> | 2005-10-27 15:34:30 +0000 |
---|---|---|
committer | baiotti <baiotti@54511f98-0e4f-0410-826e-eb8b393f5a1e> | 2005-10-27 15:34:30 +0000 |
commit | 049ce8a3bf87bbebd787dc781b76e4320e5c7231 (patch) | |
tree | 4b3513efe67bf0c93f00aaa8f24cecce5a6c823c | |
parent | a770662a2f1becd61ab39302aa5636b7a9771427 (diff) |
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
-rw-r--r-- | interface.ccl | 2 | ||||
-rw-r--r-- | schedule.ccl | 44 | ||||
-rw-r--r-- | src/loop.c | 15 | ||||
-rw-r--r-- | src/surface_integral.c | 179 | ||||
-rw-r--r-- | src/volume_integral.c | 112 |
5 files changed, 236 insertions, 116 deletions
diff --git a/interface.ccl b/interface.ccl index 05bed81..c71e1f8 100644 --- a/interface.ccl +++ b/interface.ccl @@ -29,7 +29,7 @@ CCTK_REAL ADMMass_Masses[ADMMass_number] type = SCALAR tags='checkpoint="no"' ADMMass_VolumeMass } "ADMMass Scalars" -CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar" checkpoint="no"' +CCTK_REAL ADMMass_GFs type = GF Timelevels = 3 tags='Prolongation="none" tensortypealias="Scalar" checkpoint="no"' { ADMMass_SurfaceMass_GF ADMMass_VolumeMass_pot_x diff --git a/schedule.ccl b/schedule.ccl index 437ebeb..42ebdd0 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -3,50 +3,66 @@ STORAGE:ADMMass_LoopCounterG[1] STORAGE:ADMMass_Masses[1] -STORAGE:ADMMass_GFs[1] -schedule ADMMass_InitLoopCounter AT POSTSTEP AFTER OutsideMask_UpdateMask +schedule ADMMass_InitLoopCounter AT INITIAL { LANG: C OPTIONS: global -} "Initialise loop counter" +} "Initialise the loop counter for ADMMass" -schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter +schedule ADMMass_SetLoopCounter AT POSTSTEP AFTER OutsideMask_UpdateMask { - STORAGE:ADMMass_box -} "ADMMass loop" + LANG: C + OPTIONS: global +} "Set the loop counter to the value of the parameter ADMMass:ADMMass_number" -schedule ADMMass_Loop IN ADMMass +schedule ADMMass_Loop IN ADMMass BEFORE ADMMass_Local { LANG: C OPTIONS: global } "Decrement loop counter" + +schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter +{ + STORAGE:ADMMass_GFs[3] + STORAGE:ADMMass_box +} "ADMMass loop" + schedule ADMMass_Surface IN ADMMass { LANG: C -} "Calculate the ADMmass using a surface integral" + OPTIONS: global loop-local +} "Calculate the ADMmass using a surface integral: local routine" + schedule ADMMass_Surface_Global IN ADMMass AFTER ADMMass_Surface { LANG: C OPTIONS: global -} "Calculate the ADMmass using a surface integral" +} "Calculate the ADMmass using a surface integral: global routine" + schedule ADMMass_Surface_Lapse IN ADMMass AFTER ADMMass_Surface_Global { LANG: C -} "Calculate the ADMmass using a surface integral" + OPTIONS: global loop-local +} "Calculate the ADMmass*lapse using a surface integral: local routine" + schedule ADMMass_Surface_Lapse_Global IN ADMMass AFTER ADMMass_Surface_Lapse { LANG: C OPTIONS: global -} "Calculate the ADMmass using a surface integral" +} "Calculate the ADMmass*lapse using a surface integral: global routine" + schedule ADMMass_Volume IN ADMMass { LANG: C -} "Calculate the ADMmass using a volume integral" -schedule ADMMass_Volume_Global IN ADMMass AFTER ADMMass_Volume + OPTIONS: global loop-local + SYNC: ADMMass_GFs +} "Calculate the ADMmass using a volume integral: local routine" + +schedule ADMMass_Volume_Global IN ADMMass AFTER ADMMass_Volume { LANG: C OPTIONS: global -} "Calculate the ADMmass using a volume integral" +} "Calculate the ADMmass using a volume integral: global routine" @@ -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) { |