diff options
-rw-r--r-- | interface.ccl | 2 | ||||
-rw-r--r-- | schedule.ccl | 10 | ||||
-rw-r--r-- | src/surface_integral.c | 72 | ||||
-rw-r--r-- | src/volume_integral.c | 211 |
4 files changed, 188 insertions, 107 deletions
diff --git a/interface.ccl b/interface.ccl index c71e1f8..384c54d 100644 --- a/interface.ccl +++ b/interface.ccl @@ -47,3 +47,5 @@ CCTK_REAL ADMMass_box type = scalar tags='checkpoint="no"' ADMMass_box_z_min ADMMass_box_z_max } "Physical coordinates of the surface on which the integral is computed" + +real grid_spacing_product type=SCALAR tags='checkpoint="no"' "product of cctk_delta_space, to be computed in local mode and later used in global mode (carpet problems)" diff --git a/schedule.ccl b/schedule.ccl index 42ebdd0..fc54cb1 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -22,11 +22,19 @@ schedule ADMMass_Loop IN ADMMass BEFORE ADMMass_Local OPTIONS: global } "Decrement loop counter" - +######################################################################################### +# We must schedule the local routines to compute the integrals in global-loop-local +# (as opposed to local) mode, in order to make sure that the scheduling condition +# "AFTER", which describes the ADMMass group, is respected. This group may depend +# on local routines (excision, emask) and must be run after all local routines. +# If we had no such a dependence, we could have scheduled the integral computations +# simply in local mode. +########################################################################################## schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter { STORAGE:ADMMass_GFs[3] STORAGE:ADMMass_box + STORAGE:grid_spacing_product } "ADMMass loop" schedule ADMMass_Surface IN ADMMass diff --git a/src/surface_integral.c b/src/surface_integral.c index fad504b..b8a3998 100644 --- a/src/surface_integral.c +++ b/src/surface_integral.c @@ -43,9 +43,9 @@ void ADMMass_Surface(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT i,j,k, ijk, ierr, ghost; + CCTK_INT i,j,k, ijk, ierr, ghost, ti, tj, tk, tl; CCTK_INT x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k; - CCTK_REAL detg, idetg, oneDX, oneDY, oneDZ, twoDX, twoDY, twoDZ, ds[3]; + CCTK_REAL detg, idetg, ds[3]; CCTK_REAL u[3][3], dg[3][3][3]; CCTK_REAL physical_min[3]; @@ -58,18 +58,22 @@ void ADMMass_Surface(CCTK_ARGUMENTS) char coordinate_parameter_type[9]; +#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" + /* 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]; /* 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); + const CCTK_REAL OneOverTwoDX = 1.0 / (2.0 * CCTK_DELTA_SPACE(0)); + const CCTK_REAL OneOverTwoDY = 1.0 / (2.0 * CCTK_DELTA_SPACE(1)); + const CCTK_REAL OneOverTwoDZ = 1.0 / (2.0 * CCTK_DELTA_SPACE(2)); + + const CCTK_REAL oneDX = CCTK_DELTA_SPACE(0); + const CCTK_REAL oneDY = CCTK_DELTA_SPACE(1); + const CCTK_REAL oneDZ = CCTK_DELTA_SPACE(2); + x_min_i = -INT_MAX; x_max_i = INT_MAX; @@ -156,8 +160,6 @@ void ADMMass_Surface(CCTK_ARGUMENTS) 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" - 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++) @@ -216,29 +218,29 @@ void ADMMass_Surface(CCTK_ARGUMENTS) u[2][1] = UPPERMET_UYZ; u[2][2] = UPPERMET_UZZ; - 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][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) * OneOverTwoDX; + dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) * OneOverTwoDY; + dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) * OneOverTwoDZ; - 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][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) * OneOverTwoDX; + dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) * OneOverTwoDY; + dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) * OneOverTwoDZ; - 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[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) * OneOverTwoDX; + dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) * OneOverTwoDY; + dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) * OneOverTwoDZ; 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][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) * OneOverTwoDX; + dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) * OneOverTwoDY; + dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) * OneOverTwoDZ; - 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[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) * OneOverTwoDX; + dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) * OneOverTwoDY; + dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) * OneOverTwoDZ; dg[2][0][0] = dg[0][2][0]; dg[2][0][1] = dg[0][2][1]; @@ -248,14 +250,14 @@ void ADMMass_Surface(CCTK_ARGUMENTS) 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++) + dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) * OneOverTwoDX; + dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) * OneOverTwoDY; + dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) * OneOverTwoDZ; + + for (ti = 0; ti < 3; ti++) + for (tj = 0; tj < 3; tj++) + for (tk = 0; tk < 3; tk++) + for (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]; @@ -288,7 +290,7 @@ void ADMMass_Surface_Global(CCTK_ARGUMENTS) ADMMass_SurfaceMass[*ADMMass_LoopCounter] /= 16.0*PI; CCTK_VInfo(CCTK_THORNSTRING, - "detector %d: ADM mass %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n", + "detector %d: ADM mass %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g", *ADMMass_LoopCounter, ADMMass_SurfaceMass[*ADMMass_LoopCounter], *ADMMass_box_x_min,*ADMMass_box_x_max, @@ -345,7 +347,7 @@ void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS) ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter] /= 16.0*PI; CCTK_VInfo(CCTK_THORNSTRING, - "detector %d: ADM mass with lapse %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n", + "detector %d: ADM mass with lapse %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g", *ADMMass_LoopCounter, ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter], *ADMMass_box_x_min,*ADMMass_box_x_max, diff --git a/src/volume_integral.c b/src/volume_integral.c index ee4009d..b64feaa 100644 --- a/src/volume_integral.c +++ b/src/volume_integral.c @@ -25,15 +25,16 @@ void ADMMass_Volume(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT i,j,k, ijk, ghost; + CCTK_INT i,j,k, ijk, ghost, ti, tj, tk, tl; CCTK_REAL detg, idetg; CCTK_REAL u[3][3], dg[3][3][3]; - CCTK_REAL radius, twoDX, twoDY, twoDZ; + CCTK_REAL radius; CCTK_INT mask_descriptor, state_descriptor_outside; - mask_descriptor = SpaceMask_GetTypeBits("OutsideMask"); - state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside"); + int avoid_origin_parameter; + +#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" /* grid-function strides */ const CCTK_INT di = 1; @@ -41,9 +42,12 @@ void ADMMass_Volume(CCTK_ARGUMENTS) 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); + const CCTK_REAL OneOverTwoDX = 1.0 / (2.0 * CCTK_DELTA_SPACE(0)); + const CCTK_REAL OneOverTwoDY = 1.0 / (2.0 * CCTK_DELTA_SPACE(1)); + const CCTK_REAL OneOverTwoDZ = 1.0 / (2.0 * CCTK_DELTA_SPACE(2)); + + mask_descriptor = SpaceMask_GetTypeBits("OutsideMask"); + state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside"); if (ADMMass_use_surface_distance_as_volume_radius && (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0)) @@ -71,8 +75,6 @@ void ADMMass_Volume(CCTK_ARGUMENTS) ghost = *(const int *)ghost_ptr; } -#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" - 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++) @@ -94,29 +96,29 @@ void ADMMass_Volume(CCTK_ARGUMENTS) u[2][1] = UPPERMET_UYZ; u[2][2] = UPPERMET_UZZ; - 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][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) * OneOverTwoDX; + dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) * OneOverTwoDY; + dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) * OneOverTwoDZ; - 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][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) * OneOverTwoDX; + dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) * OneOverTwoDY; + dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) * OneOverTwoDZ; - 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[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) * OneOverTwoDX; + dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) * OneOverTwoDY; + dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) * OneOverTwoDZ; 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][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) * OneOverTwoDX; + dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) * OneOverTwoDY; + dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) * OneOverTwoDZ; - 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[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) * OneOverTwoDX; + dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) * OneOverTwoDY; + dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) * OneOverTwoDZ; dg[2][0][0] = dg[0][2][0]; dg[2][0][1] = dg[0][2][1]; @@ -125,42 +127,14 @@ void ADMMass_Volume(CCTK_ARGUMENTS) 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; - 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++) + + dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) * OneOverTwoDX; + dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) * OneOverTwoDY; + dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) * OneOverTwoDZ; + + for (ti = 0; ti < 3; ti++) + for (tj = 0; tj < 3; tj++) + for (tk = 0; tk < 3; tk++) { ADMMass_VolumeMass_pot_x[ijk] += u[ti][tj] * u[tk][0] * @@ -198,16 +172,111 @@ void ADMMass_Volume(CCTK_ARGUMENTS) { 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_x[CCTK_GFINDEX3D(cctkGH,i-1,j,k)])*OneOverTwoDX+ (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_y[CCTK_GFINDEX3D(cctkGH,i,j-1,k)])*OneOverTwoDY+ (ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]- - ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)])/twoDZ); + ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)])*OneOverTwoDZ); + } + } + + /* Carpet does the following itself, but pugh does not. */ + + if (CCTK_IsThornActive("PUGH")) + { + /*Recompute the integrand on the symmetry and physical boundaries: devide by 2 on the boundary faces + cctk_lbnd values start from zero, so must add one to get the total number of points. + If cctk_ubnd[2]=n, it means that it is the (n+1)-st point (C notation) */ + + /* for short hand, right and left modified stencil values */ + const int lst = ghost; + const int rst1 = cctk_lsh[1] - ghost; + const int rst2 = cctk_lsh[2] - ghost; + const int rst3 = cctk_lsh[3] - ghost; + + /*find the value of the "avoid_origin" parameter */ + const void *avoid_origin_parameter_ptr = CCTK_ParameterValString("avoid_origin","CartGrid3D"); + assert( avoid_origin_parameter_ptr != NULL ); + avoid_origin_parameter = *(const CCTK_INT *)avoid_origin_parameter_ptr; + + printf("avoid origin %d\n",avoid_origin_parameter); + + /* if avoid_origin = yes (default), then do not devide by two. + This gives the correct result for the trapezoidal integration rule on the symmetry boundaries, + while (even if it is wrong not to devide by two the points on the faces that are + physical outer boundaries) it should not spoil the result. */ + + if (! avoid_origin_parameter) + { + if (cctk_lbnd[2] == 0) + { + for(j = lst; j <= rst2; j++) + for(i = lst; i <= rst1; i++) + { + ijk = CCTK_GFINDEX3D(cctkGH, i, j, lst); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } } - } + + if (cctk_ubnd[2]+1 == cctk_gsh[2]) + { + for(j = lst; j <= rst2; j++) + for(i = lst; i <= rst1; i++) + { + ijk = CCTK_GFINDEX3D(cctkGH, i, j, rst3); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } + } + + if (cctk_lbnd[1] == 0) + { + for(k = lst; k<= rst3; k++) + for(i = lst; i<= rst1; i++) + { + ijk = CCTK_GFINDEX3D(cctkGH, i, lst, k); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } + } + + if (cctk_ubnd[1]+1 == cctk_gsh[1]) + { + for(k = lst; k<= rst3; k++) + for(i = lst; i<= rst1; i++) + { + ijk = CCTK_GFINDEX3D(cctkGH, i, rst2, k); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } + } + + if (cctk_lbnd[0] == 0) + { + for(k = lst; k <= rst3; k++) + for(j = lst; j <= rst2; j++) + { + ijk = CCTK_GFINDEX3D(cctkGH, lst, j, k); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } + } + + if (cctk_ubnd[0]+1 == cctk_gsh[0]) + { + for(k = lst; k <= rst3; k++) + for(j = lst; j <= rst2; j++) + { + ijk = CCTK_GFINDEX3D(cctkGH, rst1, j, k); + ADMMass_VolumeMass_GF[ijk] = ADMMass_VolumeMass_GF[ijk] * 0.5; + } + } + + } /* end if avoid_origin */ + + } /* end if PUGH */ + + *grid_spacing_product = cctk_delta_space[0]*cctk_delta_space[1]*cctk_delta_space[2]; + #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h" -} +} void ADMMass_Volume_Global(CCTK_ARGUMENTS) { @@ -226,14 +295,14 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS) &ADMMass_VolumeMass[*ADMMass_LoopCounter], 1, CCTK_VarIndex("ADMMass::ADMMass_VolumeMass_GF"))) CCTK_WARN(0, "Error while reducing ADMMass_VolumeMass_GF"); - + ADMMass_VolumeMass[*ADMMass_LoopCounter] *= - cctk_delta_space[0]*cctk_delta_space[1]*cctk_delta_space[2] / (16.0*PI); + *grid_spacing_product / (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", + "detector %d: ADM mass %g, volume (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g", *ADMMass_LoopCounter, ADMMass_SurfaceMass[*ADMMass_LoopCounter], *ADMMass_box_x_min,*ADMMass_box_x_max, @@ -242,14 +311,14 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS) } else if (ADMMass_use_all_volume_as_volume_radius) { - CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume: the whole grid\n", + CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume: the whole grid", *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter]); } else { - CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume (sphere of radius) %g\n", + CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume (sphere of radius) %g", *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter], ADMMass_volume_radius[*ADMMass_LoopCounter]); } -} +} |