From 4c70ec64ffd2a662b2d6bf0b017f41711971a960 Mon Sep 17 00:00:00 2001 From: Maxwell Rizzo Date: Wed, 11 Mar 2026 16:04:22 -0500 Subject: [PATCH 1/3] VolumeIntegrals_vacuumX: fixing ADM_Momentum and ADM_Angular integrands/ls --- .../src/evaluate_integrands_local.cxx | 30 +- VolumeIntegrals_vacuumX/src/integrands.hxx | 337 ++++++++---------- .../src/number_of_reductions.cxx | 4 +- 3 files changed, 163 insertions(+), 208 deletions(-) diff --git a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx index 337e2e4..0774f8b 100644 --- a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx +++ b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx @@ -265,20 +265,13 @@ extern "C" void VI_vacuumX_ComputeIntegrand(CCTK_ARGUMENTS) { const CCTK_REAL idz = 1.0 / CCTK_DELTA_SPACE(2); grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 1, + grid.nghostzones, 2, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - VI_vacuumX_ADM_Momentum_integrand_eval_derivs( - VolIntegrand2, VolIntegrand3, VolIntegrand4, p, idx, idy, idz, + VI_vacuumX_ADM_Momentum_integrand( + VolIntegrand1, VolIntegrand2, VolIntegrand3, p, idx, idy, idz, alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, kzz); - }); - - grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 2, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - VI_vacuumX_ADM_Momentum_integrand(VolIntegrand1, p, idx, idy, idz, - VolIntegrand2, VolIntegrand3, - VolIntegrand4); + VolIntegrand4(p.I) = 0.0; }); } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], @@ -288,21 +281,14 @@ extern "C" void VI_vacuumX_ComputeIntegrand(CCTK_ARGUMENTS) { const CCTK_REAL idy = 1.0 / CCTK_DELTA_SPACE(1); const CCTK_REAL idz = 1.0 / CCTK_DELTA_SPACE(2); - grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 1, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - VI_vacuumX_ADM_Angular_Momentum_integrand_eval_derivs( - VolIntegrand2, VolIntegrand3, VolIntegrand4, p, idx, idy, idz, - alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, - kzz); - }); - grid.loop_allmn_device<1, 1, 1>( grid.nghostzones, 2, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { VI_vacuumX_ADM_Angular_Momentum_integrand( - VolIntegrand1, p, idx, idy, idz, VolIntegrand2, VolIntegrand3, - VolIntegrand4); + VolIntegrand1, VolIntegrand2, VolIntegrand3, p, idx, idy, idz, + alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, + kzz); + VolIntegrand4(p.I) = 0.0; }); } else { diff --git a/VolumeIntegrals_vacuumX/src/integrands.hxx b/VolumeIntegrals_vacuumX/src/integrands.hxx index 22f5c5b..775506d 100644 --- a/VolumeIntegrals_vacuumX/src/integrands.hxx +++ b/VolumeIntegrals_vacuumX/src/integrands.hxx @@ -4,6 +4,7 @@ #include #include +#include #include using namespace Loop; @@ -237,23 +238,17 @@ VI_vacuumX_ADM_Mass_integrand( } /* ADM Momentum */ -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -VI_vacuumX_ADM_Momentum_integrand_eval_derivs( - const GF3D2 ADM_Py_integrand_x, - const GF3D2 ADM_Py_integrand_y, - const GF3D2 ADM_Py_integrand_z, const PointDesc &p, - const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +VI_vacuumX_ADM_surface_data_at( + const PointDesc &p, const vect &index, const GF3D2 alp, const GF3D2 gxx, const GF3D2 gxy, const GF3D2 gxz, const GF3D2 gyy, const GF3D2 gyz, const GF3D2 gzz, const GF3D2 kxx, const GF3D2 kxy, const GF3D2 kxz, const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz) { - - const auto index = p.I; - - // Read in gamma_{i j} (vertex-centered -> cell-centered average) + const GF3D2 kzz, CCTK_REAL &alp_cc, CCTK_REAL &detgL, + CCTK_REAL ginv[3][3], CCTK_REAL Kup[3][3], CCTK_REAL &K) { const CCTK_REAL g11L = VI_vacuumX_avg_v2c_at(gxx, p, index); const CCTK_REAL g12L = VI_vacuumX_avg_v2c_at(gxy, p, index); const CCTK_REAL g13L = VI_vacuumX_avg_v2c_at(gxz, p, index); @@ -261,20 +256,12 @@ VI_vacuumX_ADM_Momentum_integrand_eval_derivs( const CCTK_REAL g23L = VI_vacuumX_avg_v2c_at(gyz, p, index); const CCTK_REAL g33L = VI_vacuumX_avg_v2c_at(gzz, p, index); - // Metric determinant - const CCTK_REAL detgL = - -g13L * g13L * g22L + 2 * g12L * g13L * g23L - g11L * g23L * g23L - - g12L * g12L * g33L + g11L * g22L * g33L; + detgL = -g13L * g13L * g22L + 2 * g12L * g13L * g23L - + g11L * g23L * g23L - g12L * g12L * g33L + g11L * g22L * g33L; if (!std::isfinite(detgL) || detgL <= 0.0) { - ADM_Py_integrand_x(index) = 0.0; - ADM_Py_integrand_y(index) = 0.0; - ADM_Py_integrand_z(index) = 0.0; - return; + return false; } - CCTK_REAL ginv[3][3]; - - // Calculate inverse metric gamma^{i j} ginv[0][0] = (g22L * g33L - g23L * g23L) / detgL; ginv[0][1] = (g13L * g23L - g12L * g33L) / detgL; ginv[0][2] = (g12L * g23L - g13L * g22L) / detgL; @@ -287,8 +274,6 @@ VI_vacuumX_ADM_Momentum_integrand_eval_derivs( ginv[2][1] = ginv[1][2]; CCTK_REAL Kdown[3][3]; - - // Read in covariant extrinsic curvature K_{i j} Kdown[0][0] = VI_vacuumX_avg_v2c_at(kxx, p, index); Kdown[0][1] = VI_vacuumX_avg_v2c_at(kxy, p, index); Kdown[0][2] = VI_vacuumX_avg_v2c_at(kxz, p, index); @@ -300,73 +285,56 @@ VI_vacuumX_ADM_Momentum_integrand_eval_derivs( Kdown[2][0] = Kdown[0][2]; Kdown[2][1] = Kdown[1][2]; - CCTK_REAL Kup[3][3]; - - // Calculate contravariant extrinsic curvature K^{i j} - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { - Kup[i2][j2] = 0; - - for (int k2 = 0; k2 < 3; k2++) { - for (int l2 = 0; l2 < 3; l2++) { + for (int i2 = 0; i2 < 3; ++i2) { + for (int j2 = 0; j2 < 3; ++j2) { + Kup[i2][j2] = 0.0; + for (int k2 = 0; k2 < 3; ++k2) { + for (int l2 = 0; l2 < 3; ++l2) { Kup[i2][j2] += ginv[i2][k2] * ginv[j2][l2] * Kdown[k2][l2]; } } } } - CCTK_REAL K = 0; - - // Calculate mean curvature K = g^{i j} K_{i j} - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { + K = 0.0; + for (int i2 = 0; i2 < 3; ++i2) { + for (int j2 = 0; j2 < 3; ++j2) { K += ginv[i2][j2] * Kdown[i2][j2]; } } - CCTK_REAL surface_integrand[3][3]; + alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); + return true; +} - const CCTK_REAL alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { - surface_integrand[i2][j2] = - alp_cc * sqrt(detgL) * (Kup[i2][j2] - ginv[i2][j2] * K); - } +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::array +VI_vacuumX_ADM_Momentum_surface_integrand( + const PointDesc &p, const vect &index, + const GF3D2 alp, const GF3D2 gxx, + const GF3D2 gxy, const GF3D2 gxz, + const GF3D2 gyy, const GF3D2 gyz, + const GF3D2 gzz, const GF3D2 kxx, + const GF3D2 kxy, const GF3D2 kxz, + const GF3D2 kyy, const GF3D2 kyz, + const GF3D2 kzz, const int component) { + CCTK_REAL alp_cc, detgL, ginv[3][3], Kup[3][3], K; + if (!VI_vacuumX_ADM_surface_data_at(p, index, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, + alp_cc, detgL, ginv, Kup, K)) { + return {0.0, 0.0, 0.0}; } - ADM_Py_integrand_x(index) = surface_integrand[1][0]; - ADM_Py_integrand_y(index) = surface_integrand[1][1]; - ADM_Py_integrand_z(index) = surface_integrand[1][2]; + const CCTK_REAL prefactor = alp_cc * sqrt(detgL); + return {prefactor * (Kup[component][0] - ginv[component][0] * K), + prefactor * (Kup[component][1] - ginv[component][1] * K), + prefactor * (Kup[component][2] - ginv[component][2] * K)}; } CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void VI_vacuumX_ADM_Momentum_integrand( - const GF3D2 ADM_Py_integrand, const PointDesc &p, - const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, - const GF3D2 ADM_Py_integrand_x, - const GF3D2 ADM_Py_integrand_y, - const GF3D2 ADM_Py_integrand_z) { - const CCTK_REAL cm1 = -0.5; - - ADM_Py_integrand(p.I) = - 0.125 / M_PI * - ((cm1 * (ADM_Py_integrand_x(p.I - p.DI[0]) - - ADM_Py_integrand_x(p.I + p.DI[0]))) * - idx + - (cm1 * (ADM_Py_integrand_y(p.I - p.DI[1]) - - ADM_Py_integrand_y(p.I + p.DI[1]))) * - idy + - (cm1 * (ADM_Py_integrand_z(p.I - p.DI[2]) - - ADM_Py_integrand_z(p.I + p.DI[2]))) * - idz); -} - -/* ADM Angular Momentum */ -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -VI_vacuumX_ADM_Angular_Momentum_integrand_eval_derivs( - const GF3D2 ADM_Jz_integrand_x, - const GF3D2 ADM_Jz_integrand_y, - const GF3D2 ADM_Jz_integrand_z, const PointDesc &p, + const GF3D2 ADM_Px_integrand, + const GF3D2 ADM_Py_integrand, + const GF3D2 ADM_Pz_integrand, const PointDesc &p, const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, const GF3D2 alp, const GF3D2 gxx, const GF3D2 gxy, const GF3D2 gxz, @@ -375,132 +343,133 @@ VI_vacuumX_ADM_Angular_Momentum_integrand_eval_derivs( const GF3D2 kxy, const GF3D2 kxz, const GF3D2 kyy, const GF3D2 kyz, const GF3D2 kzz) { + const CCTK_REAL cm1 = -0.5; - const auto index = p.I; - - // Read in gamma_{i j} (vertex-centered -> cell-centered average) - const CCTK_REAL g11L = VI_vacuumX_avg_v2c_at(gxx, p, index); - const CCTK_REAL g12L = VI_vacuumX_avg_v2c_at(gxy, p, index); - const CCTK_REAL g13L = VI_vacuumX_avg_v2c_at(gxz, p, index); - const CCTK_REAL g22L = VI_vacuumX_avg_v2c_at(gyy, p, index); - const CCTK_REAL g23L = VI_vacuumX_avg_v2c_at(gyz, p, index); - const CCTK_REAL g33L = VI_vacuumX_avg_v2c_at(gzz, p, index); - - // Metric determinant - const CCTK_REAL detgL = - -g13L * g13L * g22L + 2 * g12L * g13L * g23L - g11L * g23L * g23L - - g12L * g12L * g33L + g11L * g22L * g33L; - if (!std::isfinite(detgL) || detgL <= 0.0) { - ADM_Jz_integrand_x(index) = 0.0; - ADM_Jz_integrand_y(index) = 0.0; - ADM_Jz_integrand_z(index) = 0.0; - return; - } - - CCTK_REAL ginv[3][3]; - - // Calculate inverse metric gamma^{i j} - ginv[0][0] = (g22L * g33L - g23L * g23L) / detgL; - ginv[0][1] = (g13L * g23L - g12L * g33L) / detgL; - ginv[0][2] = (g12L * g23L - g13L * g22L) / detgL; - ginv[1][1] = (g11L * g33L - g13L * g13L) / detgL; - ginv[1][2] = (g12L * g13L - g11L * g23L) / detgL; - ginv[2][2] = (g11L * g22L - g12L * g12L) / detgL; - - ginv[1][0] = ginv[0][1]; - ginv[2][0] = ginv[0][2]; - ginv[2][1] = ginv[1][2]; - - CCTK_REAL Kdown[3][3]; - - // Read in covariant extrinsic curvature K_{i j} - Kdown[0][0] = VI_vacuumX_avg_v2c_at(kxx, p, index); - Kdown[0][1] = VI_vacuumX_avg_v2c_at(kxy, p, index); - Kdown[0][2] = VI_vacuumX_avg_v2c_at(kxz, p, index); - Kdown[1][1] = VI_vacuumX_avg_v2c_at(kyy, p, index); - Kdown[1][2] = VI_vacuumX_avg_v2c_at(kyz, p, index); - Kdown[2][2] = VI_vacuumX_avg_v2c_at(kzz, p, index); - - Kdown[1][0] = Kdown[0][1]; - Kdown[2][0] = Kdown[0][2]; - Kdown[2][1] = Kdown[1][2]; - - CCTK_REAL Kup[3][3]; - - // Calculate contravariant extrinsic curvature K^{i j} - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { - Kup[i2][j2] = 0; - - for (int k2 = 0; k2 < 3; k2++) { - for (int l2 = 0; l2 < 3; l2++) { - Kup[i2][j2] += ginv[i2][k2] * ginv[j2][l2] * Kdown[k2][l2]; - } - } + for (int component = 0; component < 3; ++component) { + const auto sxm = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I - p.DI[0], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + const auto sxp = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I + p.DI[0], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + const auto sym = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I - p.DI[1], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + const auto syp = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I + p.DI[1], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + const auto szm = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I - p.DI[2], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + const auto szp = VI_vacuumX_ADM_Momentum_surface_integrand( + p, p.I + p.DI[2], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, + kyy, kyz, kzz, component); + + const CCTK_REAL value = + 0.125 / M_PI * + ((cm1 * (sxm[0] - sxp[0])) * idx + (cm1 * (sym[1] - syp[1])) * idy + + (cm1 * (szm[2] - szp[2])) * idz); + + if (component == 0) { + ADM_Px_integrand(p.I) = value; + } else if (component == 1) { + ADM_Py_integrand(p.I) = value; + } else { + ADM_Pz_integrand(p.I) = value; } } +} - CCTK_REAL K = 0; - - // Calculate mean curvature K = g^{i j} K_{i j} - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { - K += ginv[i2][j2] * Kdown[i2][j2]; - } +/* ADM Angular Momentum */ +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::array +VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + const PointDesc &p, const vect &index, const CCTK_REAL xcoord, + const CCTK_REAL ycoord, const CCTK_REAL zcoord, + const GF3D2 alp, const GF3D2 gxx, + const GF3D2 gxy, const GF3D2 gxz, + const GF3D2 gyy, const GF3D2 gyz, + const GF3D2 gzz, const GF3D2 kxx, + const GF3D2 kxy, const GF3D2 kxz, + const GF3D2 kyy, const GF3D2 kyz, + const GF3D2 kzz, const int component) { + CCTK_REAL alp_cc, detgL, ginv[3][3], Kup[3][3], K; + if (!VI_vacuumX_ADM_surface_data_at(p, index, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, + alp_cc, detgL, ginv, Kup, K)) { + return {0.0, 0.0, 0.0}; } - // Levi-Civita tensor - CCTK_REAL LCT[3][3][3] = { - {{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, - {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, - {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; - - // Coordinate vector - CCTK_REAL xx[3] = {p.x, p.y, p.z}; - - CCTK_REAL surface_integrand[3][3]; - - const CCTK_REAL alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); - for (int i2 = 0; i2 < 3; i2++) { - for (int j2 = 0; j2 < 3; j2++) { - surface_integrand[i2][j2] = 0; - - for (int k2 = 0; k2 < 3; k2++) { - for (int l2 = 0; l2 < 3; l2++) { - surface_integrand[i2][j2] += - LCT[i2][k2][l2] * xx[k2] * (Kup[l2][j2] - ginv[l2][j2] * K); - } + const CCTK_REAL LCT[3][3][3] = { + {{0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, -1.0, 0.0}}, + {{0.0, 0.0, -1.0}, {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}}, + {{0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}}; + const CCTK_REAL xx[3] = {xcoord, ycoord, zcoord}; + + std::array result{0.0, 0.0, 0.0}; + const CCTK_REAL prefactor = alp_cc * sqrt(detgL); + for (int j2 = 0; j2 < 3; ++j2) { + for (int k2 = 0; k2 < 3; ++k2) { + for (int l2 = 0; l2 < 3; ++l2) { + result[j2] += + LCT[component][k2][l2] * xx[k2] * (Kup[l2][j2] - ginv[l2][j2] * K); } - - surface_integrand[i2][j2] *= alp_cc * sqrt(detgL); } + result[j2] *= prefactor; } - - ADM_Jz_integrand_x(index) = surface_integrand[2][0]; - ADM_Jz_integrand_y(index) = surface_integrand[2][1]; - ADM_Jz_integrand_z(index) = surface_integrand[2][2]; + return result; } CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void VI_vacuumX_ADM_Angular_Momentum_integrand( + const GF3D2 ADM_Jx_integrand, + const GF3D2 ADM_Jy_integrand, const GF3D2 ADM_Jz_integrand, const PointDesc &p, const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, - const GF3D2 ADM_Jz_integrand_x, - const GF3D2 ADM_Jz_integrand_y, - const GF3D2 ADM_Jz_integrand_z) { + const GF3D2 alp, const GF3D2 gxx, + const GF3D2 gxy, const GF3D2 gxz, + const GF3D2 gyy, const GF3D2 gyz, + const GF3D2 gzz, const GF3D2 kxx, + const GF3D2 kxy, const GF3D2 kxz, + const GF3D2 kyy, const GF3D2 kyz, + const GF3D2 kzz) { const CCTK_REAL cm1 = -0.5; - - ADM_Jz_integrand(p.I) = - 0.125 / M_PI * - ((cm1 * (ADM_Jz_integrand_x(p.I - p.DI[0]) - - ADM_Jz_integrand_x(p.I + p.DI[0]))) * - idx + - (cm1 * (ADM_Jz_integrand_y(p.I - p.DI[1]) - - ADM_Jz_integrand_y(p.I + p.DI[1]))) * - idy + - (cm1 * (ADM_Jz_integrand_z(p.I - p.DI[2]) - - ADM_Jz_integrand_z(p.I + p.DI[2]))) * - idz); + const CCTK_REAL dx = 1.0 / idx; + const CCTK_REAL dy = 1.0 / idy; + const CCTK_REAL dz = 1.0 / idz; + + for (int component = 0; component < 3; ++component) { + const auto sxm = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I - p.DI[0], p.x - dx, p.y, p.z, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + const auto sxp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I + p.DI[0], p.x + dx, p.y, p.z, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + const auto sym = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I - p.DI[1], p.x, p.y - dy, p.z, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + const auto syp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I + p.DI[1], p.x, p.y + dy, p.z, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + const auto szm = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I - p.DI[2], p.x, p.y, p.z - dz, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + const auto szp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( + p, p.I + p.DI[2], p.x, p.y, p.z + dz, alp, gxx, gxy, gxz, gyy, gyz, + gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); + + const CCTK_REAL value = + 0.125 / M_PI * + ((cm1 * (sxm[0] - sxp[0])) * idx + (cm1 * (sym[1] - syp[1])) * idy + + (cm1 * (szm[2] - szp[2])) * idz); + + if (component == 0) { + ADM_Jx_integrand(p.I) = value; + } else if (component == 1) { + ADM_Jy_integrand(p.I) = value; + } else { + ADM_Jz_integrand(p.I) = value; + } + } } #endif // INTEGRANDS__VOLUMEINTEGRALS_VACUUMX__ diff --git a/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx b/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx index 931e890..c10206d 100644 --- a/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx +++ b/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx @@ -32,10 +32,10 @@ VI_vacuumX_number_of_reductions(int which_integral) { return 1; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "ADM_Momentum")) - return 1; + return 3; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "ADM_Angular_Momentum")) - return 1; + return 3; CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, "VolumeIntegrals_vacuumX ERROR: You forgot to specify the number of reductions for Integration_quantity_keyword=%s!\n", From 56e519039668f2dfb3422cdc44c7557c599fe2ec Mon Sep 17 00:00:00 2001 From: Maxwell Rizzo Date: Thu, 9 Apr 2026 10:59:13 -0500 Subject: [PATCH 2/3] VolumeIntegrals_vacuumX: Schedule change and M_ADM integrand modification. Main schedule group is now enforced to run after z4c analysis group, to ensure that the Hamiltonian constraint is updated. Additionally, the M_ADM integrand now manuallay excludes the boundary/ghost zone contribution using 'grid.box_int'. --- VolumeIntegrals_vacuumX/schedule.ccl | 2 +- .../src/evaluate_integrands_local.cxx | 25 ++++++++++++++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/VolumeIntegrals_vacuumX/schedule.ccl b/VolumeIntegrals_vacuumX/schedule.ccl index ec0c820..3aceb54 100644 --- a/VolumeIntegrals_vacuumX/schedule.ccl +++ b/VolumeIntegrals_vacuumX/schedule.ccl @@ -60,7 +60,7 @@ SCHEDULE VI_vacuumX_InitializeIntegralCounter before VI_vacuumX_VolumeIntegralGr WRITES: IntegralCounter(everywhere) } "Initialize IntegralCounter variable" ################## -SCHEDULE GROUP VI_vacuumX_VolumeIntegralGroup AT CCTK_ANALYSIS WHILE VolumeIntegrals_vacuumX::IntegralCounter +SCHEDULE GROUP VI_vacuumX_VolumeIntegralGroup AT CCTK_ANALYSIS WHILE VolumeIntegrals_vacuumX::IntegralCounter AFTER Z4c_AnalysisGroup { } "Evaluate all volume integrals" diff --git a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx index 0774f8b..f435280 100644 --- a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx +++ b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx @@ -240,18 +240,35 @@ extern "C" void VI_vacuumX_ComputeIntegrand(CCTK_ARGUMENTS) { const CCTK_REAL idx = 1.0 / CCTK_DELTA_SPACE(0); const CCTK_REAL idy = 1.0 / CCTK_DELTA_SPACE(1); const CCTK_REAL idz = 1.0 / CCTK_DELTA_SPACE(2); + vect imin, imax; + grid.box_int<1, 1, 1>(grid.nghostzones, imin, imax); - grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 1, + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + VolIntegrand1(p.I) = 0.0; + VolIntegrand2(p.I) = 0.0; + VolIntegrand3(p.I) = 0.0; + VolIntegrand4(p.I) = 0.0; + }); + + grid.loop_int_device<1, 1, 1>( + grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + for (int d = 0; d < dim; ++d) + if (p.I[d] < imin[d] + 1 || p.I[d] >= imax[d] - 1) + return; VI_vacuumX_ADM_Mass_integrand_eval_derivs( VolIntegrand2, VolIntegrand3, VolIntegrand4, p, idx, idy, idz, alp, gxx, gxy, gxz, gyy, gyz, gzz); }); - grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 2, + grid.loop_int_device<1, 1, 1>( + grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + for (int d = 0; d < dim; ++d) + if (p.I[d] < imin[d] + 2 || p.I[d] >= imax[d] - 2) + return; VI_vacuumX_ADM_Mass_integrand(VolIntegrand1, p, idx, idy, idz, VolIntegrand2, VolIntegrand3, VolIntegrand4); From 3e6f9cc5b9b52c1ed9bd0f8ce06f4696d682c95d Mon Sep 17 00:00:00 2001 From: Maxwell Rizzo Date: Thu, 9 Apr 2026 11:36:07 -0500 Subject: [PATCH 3/3] VolumeIntegrals_vacuumX: Reverting previous ADM Linear/Angular momentum vector fix --- .../src/evaluate_integrands_local.cxx | 30 +- VolumeIntegrals_vacuumX/src/integrands.hxx | 337 ++++++++++-------- .../src/number_of_reductions.cxx | 4 +- 3 files changed, 208 insertions(+), 163 deletions(-) diff --git a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx index f435280..7780df7 100644 --- a/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx +++ b/VolumeIntegrals_vacuumX/src/evaluate_integrands_local.cxx @@ -282,13 +282,20 @@ extern "C" void VI_vacuumX_ComputeIntegrand(CCTK_ARGUMENTS) { const CCTK_REAL idz = 1.0 / CCTK_DELTA_SPACE(2); grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 2, + grid.nghostzones, 1, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - VI_vacuumX_ADM_Momentum_integrand( - VolIntegrand1, VolIntegrand2, VolIntegrand3, p, idx, idy, idz, + VI_vacuumX_ADM_Momentum_integrand_eval_derivs( + VolIntegrand2, VolIntegrand3, VolIntegrand4, p, idx, idy, idz, alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, kzz); - VolIntegrand4(p.I) = 0.0; + }); + + grid.loop_allmn_device<1, 1, 1>( + grid.nghostzones, 2, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + VI_vacuumX_ADM_Momentum_integrand(VolIntegrand1, p, idx, idy, idz, + VolIntegrand2, VolIntegrand3, + VolIntegrand4); }); } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], @@ -299,13 +306,20 @@ extern "C" void VI_vacuumX_ComputeIntegrand(CCTK_ARGUMENTS) { const CCTK_REAL idz = 1.0 / CCTK_DELTA_SPACE(2); grid.loop_allmn_device<1, 1, 1>( - grid.nghostzones, 2, + grid.nghostzones, 1, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - VI_vacuumX_ADM_Angular_Momentum_integrand( - VolIntegrand1, VolIntegrand2, VolIntegrand3, p, idx, idy, idz, + VI_vacuumX_ADM_Angular_Momentum_integrand_eval_derivs( + VolIntegrand2, VolIntegrand3, VolIntegrand4, p, idx, idy, idz, alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, kzz); - VolIntegrand4(p.I) = 0.0; + }); + + grid.loop_allmn_device<1, 1, 1>( + grid.nghostzones, 2, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + VI_vacuumX_ADM_Angular_Momentum_integrand( + VolIntegrand1, p, idx, idy, idz, VolIntegrand2, VolIntegrand3, + VolIntegrand4); }); } else { diff --git a/VolumeIntegrals_vacuumX/src/integrands.hxx b/VolumeIntegrals_vacuumX/src/integrands.hxx index 775506d..22f5c5b 100644 --- a/VolumeIntegrals_vacuumX/src/integrands.hxx +++ b/VolumeIntegrals_vacuumX/src/integrands.hxx @@ -4,7 +4,6 @@ #include #include -#include #include using namespace Loop; @@ -238,17 +237,23 @@ VI_vacuumX_ADM_Mass_integrand( } /* ADM Momentum */ -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool -VI_vacuumX_ADM_surface_data_at( - const PointDesc &p, const vect &index, +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +VI_vacuumX_ADM_Momentum_integrand_eval_derivs( + const GF3D2 ADM_Py_integrand_x, + const GF3D2 ADM_Py_integrand_y, + const GF3D2 ADM_Py_integrand_z, const PointDesc &p, + const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, const GF3D2 alp, const GF3D2 gxx, const GF3D2 gxy, const GF3D2 gxz, const GF3D2 gyy, const GF3D2 gyz, const GF3D2 gzz, const GF3D2 kxx, const GF3D2 kxy, const GF3D2 kxz, const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz, CCTK_REAL &alp_cc, CCTK_REAL &detgL, - CCTK_REAL ginv[3][3], CCTK_REAL Kup[3][3], CCTK_REAL &K) { + const GF3D2 kzz) { + + const auto index = p.I; + + // Read in gamma_{i j} (vertex-centered -> cell-centered average) const CCTK_REAL g11L = VI_vacuumX_avg_v2c_at(gxx, p, index); const CCTK_REAL g12L = VI_vacuumX_avg_v2c_at(gxy, p, index); const CCTK_REAL g13L = VI_vacuumX_avg_v2c_at(gxz, p, index); @@ -256,12 +261,20 @@ VI_vacuumX_ADM_surface_data_at( const CCTK_REAL g23L = VI_vacuumX_avg_v2c_at(gyz, p, index); const CCTK_REAL g33L = VI_vacuumX_avg_v2c_at(gzz, p, index); - detgL = -g13L * g13L * g22L + 2 * g12L * g13L * g23L - - g11L * g23L * g23L - g12L * g12L * g33L + g11L * g22L * g33L; + // Metric determinant + const CCTK_REAL detgL = + -g13L * g13L * g22L + 2 * g12L * g13L * g23L - g11L * g23L * g23L - + g12L * g12L * g33L + g11L * g22L * g33L; if (!std::isfinite(detgL) || detgL <= 0.0) { - return false; + ADM_Py_integrand_x(index) = 0.0; + ADM_Py_integrand_y(index) = 0.0; + ADM_Py_integrand_z(index) = 0.0; + return; } + CCTK_REAL ginv[3][3]; + + // Calculate inverse metric gamma^{i j} ginv[0][0] = (g22L * g33L - g23L * g23L) / detgL; ginv[0][1] = (g13L * g23L - g12L * g33L) / detgL; ginv[0][2] = (g12L * g23L - g13L * g22L) / detgL; @@ -274,6 +287,8 @@ VI_vacuumX_ADM_surface_data_at( ginv[2][1] = ginv[1][2]; CCTK_REAL Kdown[3][3]; + + // Read in covariant extrinsic curvature K_{i j} Kdown[0][0] = VI_vacuumX_avg_v2c_at(kxx, p, index); Kdown[0][1] = VI_vacuumX_avg_v2c_at(kxy, p, index); Kdown[0][2] = VI_vacuumX_avg_v2c_at(kxz, p, index); @@ -285,191 +300,207 @@ VI_vacuumX_ADM_surface_data_at( Kdown[2][0] = Kdown[0][2]; Kdown[2][1] = Kdown[1][2]; - for (int i2 = 0; i2 < 3; ++i2) { - for (int j2 = 0; j2 < 3; ++j2) { - Kup[i2][j2] = 0.0; - for (int k2 = 0; k2 < 3; ++k2) { - for (int l2 = 0; l2 < 3; ++l2) { + CCTK_REAL Kup[3][3]; + + // Calculate contravariant extrinsic curvature K^{i j} + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { + Kup[i2][j2] = 0; + + for (int k2 = 0; k2 < 3; k2++) { + for (int l2 = 0; l2 < 3; l2++) { Kup[i2][j2] += ginv[i2][k2] * ginv[j2][l2] * Kdown[k2][l2]; } } } } - K = 0.0; - for (int i2 = 0; i2 < 3; ++i2) { - for (int j2 = 0; j2 < 3; ++j2) { + CCTK_REAL K = 0; + + // Calculate mean curvature K = g^{i j} K_{i j} + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { K += ginv[i2][j2] * Kdown[i2][j2]; } } - alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); - return true; -} + CCTK_REAL surface_integrand[3][3]; -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::array -VI_vacuumX_ADM_Momentum_surface_integrand( - const PointDesc &p, const vect &index, - const GF3D2 alp, const GF3D2 gxx, - const GF3D2 gxy, const GF3D2 gxz, - const GF3D2 gyy, const GF3D2 gyz, - const GF3D2 gzz, const GF3D2 kxx, - const GF3D2 kxy, const GF3D2 kxz, - const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz, const int component) { - CCTK_REAL alp_cc, detgL, ginv[3][3], Kup[3][3], K; - if (!VI_vacuumX_ADM_surface_data_at(p, index, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, - alp_cc, detgL, ginv, Kup, K)) { - return {0.0, 0.0, 0.0}; + const CCTK_REAL alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { + surface_integrand[i2][j2] = + alp_cc * sqrt(detgL) * (Kup[i2][j2] - ginv[i2][j2] * K); + } } - const CCTK_REAL prefactor = alp_cc * sqrt(detgL); - return {prefactor * (Kup[component][0] - ginv[component][0] * K), - prefactor * (Kup[component][1] - ginv[component][1] * K), - prefactor * (Kup[component][2] - ginv[component][2] * K)}; + ADM_Py_integrand_x(index) = surface_integrand[1][0]; + ADM_Py_integrand_y(index) = surface_integrand[1][1]; + ADM_Py_integrand_z(index) = surface_integrand[1][2]; } CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void VI_vacuumX_ADM_Momentum_integrand( - const GF3D2 ADM_Px_integrand, - const GF3D2 ADM_Py_integrand, - const GF3D2 ADM_Pz_integrand, const PointDesc &p, + const GF3D2 ADM_Py_integrand, const PointDesc &p, const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, - const GF3D2 alp, const GF3D2 gxx, - const GF3D2 gxy, const GF3D2 gxz, - const GF3D2 gyy, const GF3D2 gyz, - const GF3D2 gzz, const GF3D2 kxx, - const GF3D2 kxy, const GF3D2 kxz, - const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz) { + const GF3D2 ADM_Py_integrand_x, + const GF3D2 ADM_Py_integrand_y, + const GF3D2 ADM_Py_integrand_z) { const CCTK_REAL cm1 = -0.5; - for (int component = 0; component < 3; ++component) { - const auto sxm = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I - p.DI[0], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - const auto sxp = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I + p.DI[0], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - const auto sym = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I - p.DI[1], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - const auto syp = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I + p.DI[1], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - const auto szm = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I - p.DI[2], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - const auto szp = VI_vacuumX_ADM_Momentum_surface_integrand( - p, p.I + p.DI[2], alp, gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, - kyy, kyz, kzz, component); - - const CCTK_REAL value = - 0.125 / M_PI * - ((cm1 * (sxm[0] - sxp[0])) * idx + (cm1 * (sym[1] - syp[1])) * idy + - (cm1 * (szm[2] - szp[2])) * idz); - - if (component == 0) { - ADM_Px_integrand(p.I) = value; - } else if (component == 1) { - ADM_Py_integrand(p.I) = value; - } else { - ADM_Pz_integrand(p.I) = value; - } - } + ADM_Py_integrand(p.I) = + 0.125 / M_PI * + ((cm1 * (ADM_Py_integrand_x(p.I - p.DI[0]) - + ADM_Py_integrand_x(p.I + p.DI[0]))) * + idx + + (cm1 * (ADM_Py_integrand_y(p.I - p.DI[1]) - + ADM_Py_integrand_y(p.I + p.DI[1]))) * + idy + + (cm1 * (ADM_Py_integrand_z(p.I - p.DI[2]) - + ADM_Py_integrand_z(p.I + p.DI[2]))) * + idz); } /* ADM Angular Momentum */ -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::array -VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - const PointDesc &p, const vect &index, const CCTK_REAL xcoord, - const CCTK_REAL ycoord, const CCTK_REAL zcoord, +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +VI_vacuumX_ADM_Angular_Momentum_integrand_eval_derivs( + const GF3D2 ADM_Jz_integrand_x, + const GF3D2 ADM_Jz_integrand_y, + const GF3D2 ADM_Jz_integrand_z, const PointDesc &p, + const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, const GF3D2 alp, const GF3D2 gxx, const GF3D2 gxy, const GF3D2 gxz, const GF3D2 gyy, const GF3D2 gyz, const GF3D2 gzz, const GF3D2 kxx, const GF3D2 kxy, const GF3D2 kxz, const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz, const int component) { - CCTK_REAL alp_cc, detgL, ginv[3][3], Kup[3][3], K; - if (!VI_vacuumX_ADM_surface_data_at(p, index, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, - alp_cc, detgL, ginv, Kup, K)) { - return {0.0, 0.0, 0.0}; + const GF3D2 kzz) { + + const auto index = p.I; + + // Read in gamma_{i j} (vertex-centered -> cell-centered average) + const CCTK_REAL g11L = VI_vacuumX_avg_v2c_at(gxx, p, index); + const CCTK_REAL g12L = VI_vacuumX_avg_v2c_at(gxy, p, index); + const CCTK_REAL g13L = VI_vacuumX_avg_v2c_at(gxz, p, index); + const CCTK_REAL g22L = VI_vacuumX_avg_v2c_at(gyy, p, index); + const CCTK_REAL g23L = VI_vacuumX_avg_v2c_at(gyz, p, index); + const CCTK_REAL g33L = VI_vacuumX_avg_v2c_at(gzz, p, index); + + // Metric determinant + const CCTK_REAL detgL = + -g13L * g13L * g22L + 2 * g12L * g13L * g23L - g11L * g23L * g23L - + g12L * g12L * g33L + g11L * g22L * g33L; + if (!std::isfinite(detgL) || detgL <= 0.0) { + ADM_Jz_integrand_x(index) = 0.0; + ADM_Jz_integrand_y(index) = 0.0; + ADM_Jz_integrand_z(index) = 0.0; + return; + } + + CCTK_REAL ginv[3][3]; + + // Calculate inverse metric gamma^{i j} + ginv[0][0] = (g22L * g33L - g23L * g23L) / detgL; + ginv[0][1] = (g13L * g23L - g12L * g33L) / detgL; + ginv[0][2] = (g12L * g23L - g13L * g22L) / detgL; + ginv[1][1] = (g11L * g33L - g13L * g13L) / detgL; + ginv[1][2] = (g12L * g13L - g11L * g23L) / detgL; + ginv[2][2] = (g11L * g22L - g12L * g12L) / detgL; + + ginv[1][0] = ginv[0][1]; + ginv[2][0] = ginv[0][2]; + ginv[2][1] = ginv[1][2]; + + CCTK_REAL Kdown[3][3]; + + // Read in covariant extrinsic curvature K_{i j} + Kdown[0][0] = VI_vacuumX_avg_v2c_at(kxx, p, index); + Kdown[0][1] = VI_vacuumX_avg_v2c_at(kxy, p, index); + Kdown[0][2] = VI_vacuumX_avg_v2c_at(kxz, p, index); + Kdown[1][1] = VI_vacuumX_avg_v2c_at(kyy, p, index); + Kdown[1][2] = VI_vacuumX_avg_v2c_at(kyz, p, index); + Kdown[2][2] = VI_vacuumX_avg_v2c_at(kzz, p, index); + + Kdown[1][0] = Kdown[0][1]; + Kdown[2][0] = Kdown[0][2]; + Kdown[2][1] = Kdown[1][2]; + + CCTK_REAL Kup[3][3]; + + // Calculate contravariant extrinsic curvature K^{i j} + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { + Kup[i2][j2] = 0; + + for (int k2 = 0; k2 < 3; k2++) { + for (int l2 = 0; l2 < 3; l2++) { + Kup[i2][j2] += ginv[i2][k2] * ginv[j2][l2] * Kdown[k2][l2]; + } + } + } } - const CCTK_REAL LCT[3][3][3] = { - {{0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, -1.0, 0.0}}, - {{0.0, 0.0, -1.0}, {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}}, - {{0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}}; - const CCTK_REAL xx[3] = {xcoord, ycoord, zcoord}; - - std::array result{0.0, 0.0, 0.0}; - const CCTK_REAL prefactor = alp_cc * sqrt(detgL); - for (int j2 = 0; j2 < 3; ++j2) { - for (int k2 = 0; k2 < 3; ++k2) { - for (int l2 = 0; l2 < 3; ++l2) { - result[j2] += - LCT[component][k2][l2] * xx[k2] * (Kup[l2][j2] - ginv[l2][j2] * K); + CCTK_REAL K = 0; + + // Calculate mean curvature K = g^{i j} K_{i j} + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { + K += ginv[i2][j2] * Kdown[i2][j2]; + } + } + + // Levi-Civita tensor + CCTK_REAL LCT[3][3][3] = { + {{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, + {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, + {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; + + // Coordinate vector + CCTK_REAL xx[3] = {p.x, p.y, p.z}; + + CCTK_REAL surface_integrand[3][3]; + + const CCTK_REAL alp_cc = VI_vacuumX_avg_v2c_at(alp, p, index); + for (int i2 = 0; i2 < 3; i2++) { + for (int j2 = 0; j2 < 3; j2++) { + surface_integrand[i2][j2] = 0; + + for (int k2 = 0; k2 < 3; k2++) { + for (int l2 = 0; l2 < 3; l2++) { + surface_integrand[i2][j2] += + LCT[i2][k2][l2] * xx[k2] * (Kup[l2][j2] - ginv[l2][j2] * K); + } } + + surface_integrand[i2][j2] *= alp_cc * sqrt(detgL); } - result[j2] *= prefactor; } - return result; + + ADM_Jz_integrand_x(index) = surface_integrand[2][0]; + ADM_Jz_integrand_y(index) = surface_integrand[2][1]; + ADM_Jz_integrand_z(index) = surface_integrand[2][2]; } CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void VI_vacuumX_ADM_Angular_Momentum_integrand( - const GF3D2 ADM_Jx_integrand, - const GF3D2 ADM_Jy_integrand, const GF3D2 ADM_Jz_integrand, const PointDesc &p, const CCTK_REAL idx, const CCTK_REAL idy, const CCTK_REAL idz, - const GF3D2 alp, const GF3D2 gxx, - const GF3D2 gxy, const GF3D2 gxz, - const GF3D2 gyy, const GF3D2 gyz, - const GF3D2 gzz, const GF3D2 kxx, - const GF3D2 kxy, const GF3D2 kxz, - const GF3D2 kyy, const GF3D2 kyz, - const GF3D2 kzz) { + const GF3D2 ADM_Jz_integrand_x, + const GF3D2 ADM_Jz_integrand_y, + const GF3D2 ADM_Jz_integrand_z) { const CCTK_REAL cm1 = -0.5; - const CCTK_REAL dx = 1.0 / idx; - const CCTK_REAL dy = 1.0 / idy; - const CCTK_REAL dz = 1.0 / idz; - - for (int component = 0; component < 3; ++component) { - const auto sxm = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I - p.DI[0], p.x - dx, p.y, p.z, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - const auto sxp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I + p.DI[0], p.x + dx, p.y, p.z, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - const auto sym = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I - p.DI[1], p.x, p.y - dy, p.z, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - const auto syp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I + p.DI[1], p.x, p.y + dy, p.z, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - const auto szm = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I - p.DI[2], p.x, p.y, p.z - dz, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - const auto szp = VI_vacuumX_ADM_Angular_Momentum_surface_integrand( - p, p.I + p.DI[2], p.x, p.y, p.z + dz, alp, gxx, gxy, gxz, gyy, gyz, - gzz, kxx, kxy, kxz, kyy, kyz, kzz, component); - - const CCTK_REAL value = - 0.125 / M_PI * - ((cm1 * (sxm[0] - sxp[0])) * idx + (cm1 * (sym[1] - syp[1])) * idy + - (cm1 * (szm[2] - szp[2])) * idz); - - if (component == 0) { - ADM_Jx_integrand(p.I) = value; - } else if (component == 1) { - ADM_Jy_integrand(p.I) = value; - } else { - ADM_Jz_integrand(p.I) = value; - } - } + + ADM_Jz_integrand(p.I) = + 0.125 / M_PI * + ((cm1 * (ADM_Jz_integrand_x(p.I - p.DI[0]) - + ADM_Jz_integrand_x(p.I + p.DI[0]))) * + idx + + (cm1 * (ADM_Jz_integrand_y(p.I - p.DI[1]) - + ADM_Jz_integrand_y(p.I + p.DI[1]))) * + idy + + (cm1 * (ADM_Jz_integrand_z(p.I - p.DI[2]) - + ADM_Jz_integrand_z(p.I + p.DI[2]))) * + idz); } #endif // INTEGRANDS__VOLUMEINTEGRALS_VACUUMX__ diff --git a/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx b/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx index c10206d..931e890 100644 --- a/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx +++ b/VolumeIntegrals_vacuumX/src/number_of_reductions.cxx @@ -32,10 +32,10 @@ VI_vacuumX_number_of_reductions(int which_integral) { return 1; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "ADM_Momentum")) - return 3; + return 1; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "ADM_Angular_Momentum")) - return 3; + return 1; CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, "VolumeIntegrals_vacuumX ERROR: You forgot to specify the number of reductions for Integration_quantity_keyword=%s!\n",