aboutsummaryrefslogtreecommitdiff
path: root/src/volume_integral.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/volume_integral.c')
-rw-r--r--src/volume_integral.c211
1 files changed, 140 insertions, 71 deletions
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 03:41:22 +0000