aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl2
-rw-r--r--schedule.ccl10
-rw-r--r--src/surface_integral.c72
-rw-r--r--src/volume_integral.c211
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]);
}
-}
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 00:11:03 +0000