diff --git a/MHD_Diagnostics/param.ccl b/MHD_Diagnostics/param.ccl index 0ba6d3d..f9ccd47 100644 --- a/MHD_Diagnostics/param.ccl +++ b/MHD_Diagnostics/param.ccl @@ -13,6 +13,9 @@ # USES CCTK_INT metric_timelevels ############################################################################# +shares: AsterX +USES CCTK_INT mag_correction_order + private: CCTK_INT diagnostics_compute_every "How often to compute diagnostics?" STEERABLE=ALWAYS diff --git a/MHD_Diagnostics/src/master.cxx b/MHD_Diagnostics/src/master.cxx index 823ce88..befda0a 100644 --- a/MHD_Diagnostics/src/master.cxx +++ b/MHD_Diagnostics/src/master.cxx @@ -387,29 +387,11 @@ extern "C" void compute_mhd_derivs(CCTK_ARGUMENTS) { double invdet = 1.0 / det; // Calculate also divB from AsterX staggered field - // Use fourth-order ECHO expressions - // TODO: check if face index is correct, how to handle densitized B? - - double Bx_stagg = 13. * dBx_stag(p.I + p.DI[0]) / 12. - - (dBx_stag(p.I) + dBx_stag(p.I + 2 * p.DI[0])) / 24.; - double By_stagg = 13. * dBy_stag(p.I + p.DI[1]) / 12. - - (dBy_stag(p.I) + dBy_stag(p.I + 2 * p.DI[1])) / 24.; - double Bz_stagg = 13. * dBz_stag(p.I + p.DI[2]) / 12. - - (dBz_stag(p.I) + dBz_stag(p.I + 2 * p.DI[2])) / 24.; - - double Bx_stagg_m1 = - 13. * dBx_stag(p.I) / 12. - - (dBx_stag(p.I + p.DI[0]) + dBx_stag(p.I - p.DI[0])) / 24.; - double By_stagg_m1 = - 13. * dBy_stag(p.I) / 12. - - (dBy_stag(p.I + p.DI[1]) + dBy_stag(p.I - p.DI[1])) / 24.; - double Bz_stagg_m1 = - 13. * dBz_stag(p.I) / 12. - - (dBz_stag(p.I + p.DI[2]) + dBz_stag(p.I - p.DI[2])) / 24.; - - divB(p.I) = dxi[1] * (Bx_stagg - Bx_stagg_m1) + - dxi[2] * (By_stagg - By_stagg_m1) + - dxi[3] * (Bz_stagg - Bz_stagg_m1); + // using fourth-order ECHO expressions + + divB(p.I) = calc_fd_forward_midpoint<0>(dBx_stag, p, mag_correction_order) + + calc_fd_forward_midpoint<1>(dBy_stag, p, mag_correction_order) + + calc_fd_forward_midpoint<2>(dBz_stag, p, mag_correction_order); // =========================================== diff --git a/MeudonBNSX/src/Bin_NS.cxx b/MeudonBNSX/src/Bin_NS.cxx index c0bc9c1..9901aef 100644 --- a/MeudonBNSX/src/Bin_NS.cxx +++ b/MeudonBNSX/src/Bin_NS.cxx @@ -28,7 +28,7 @@ using namespace Arith; using namespace Loop; // define namespace here for old versions of Lorene that don't do so -namespace Lorene {} +// namespace Lorene {} using namespace Lorene; namespace { @@ -75,9 +75,6 @@ extern "C" void MeudonBNSX_initialise(CCTK_ARGUMENTS) { int const npoints = nx * ny * nz; vector xx(npoints), yy(npoints), zz(npoints); - // TODO: currently works only with polytropic EOS - CCTK_INFO("MeudonBNSX will use the polytropic equation of state."); - grid.loop_all<1, 1, 1>(grid.nghostzones, [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { CCTK_INT idx = CCC_layout.linear(p.I); @@ -210,7 +207,9 @@ extern "C" void MeudonBNSX_initialise(CCTK_ARGUMENTS) { rho(p.I) = bin_ns.nbar[idx] / rho_unit; if (!recalculate_eps) eps(p.I) = bin_ns.ener_spec[idx]; - // TODO: currently works only with polytropic EOS + // TODO: this only works for simple polytropic, or tabulated EOS. + // The expressions are valid for polytropes. Eps and press can and must be + // recalculated in a separate thorn for tabulated EOS. else eps(p.I) = K * pow(rho(p.I), bin_ns.gamma_poly1 - 1.) / (bin_ns.gamma_poly1 - 1.); diff --git a/VolumeIntegrals_GRMHDX/interface.ccl b/VolumeIntegrals_GRMHDX/interface.ccl index c54058d..683729b 100644 --- a/VolumeIntegrals_GRMHDX/interface.ccl +++ b/VolumeIntegrals_GRMHDX/interface.ccl @@ -6,6 +6,22 @@ inherits: ADMBaseX HydroBaseX AsterX BoxInBox USES INCLUDE HEADER: loop_device.hxx +public: +CCTK_REAL GlobalCoM type=SCALAR tags='checkpoint="no"' +{ + comx, comy, comz +} "Record CoM to set origin for toroidal/poloidal decompositions" + +CCTK_REAL CoM1 type=SCALAR tags='checkpoint="no"' +{ + comx1, comy1, comz1 +} "Record CoM of NS 1" + +CCTK_REAL CoM2 type=SCALAR tags='checkpoint="no"' +{ + comx2, comy2, comz2 +} "Record CoM of NS 2" + private: CCTK_REAL VolIntegrands TYPE=GF CENTERING={ccc} tags='checkpoint="no"' { @@ -33,7 +49,3 @@ CCTK_REAL VolIntegrals_vacuum_time type=SCALAR tags='checkpoint="no"' physical_time } "keeps track of the physical time, in case time coordinate is reparameterized, a la http://arxiv.org/abs/1404.6523" -CCTK_REAL GlobalCoM type=SCALAR tags='checkpoint="no"' -{ - comx, comy, comz -} "Record CoM to set origin for toroidal/poloidal decompositions" diff --git a/VolumeIntegrals_GRMHDX/param.ccl b/VolumeIntegrals_GRMHDX/param.ccl index ef6bee1..e331c9f 100644 --- a/VolumeIntegrals_GRMHDX/param.ccl +++ b/VolumeIntegrals_GRMHDX/param.ccl @@ -86,8 +86,16 @@ keyword Integration_quantity_keyword[101] "Which quantity to integrate" STEERABL "thermal_energy" :: "Thermal energy" "volume_average_norm_B_ab" :: "Volume average of norm B in density interval (dens_a,dens_b]" "volume_average_norm_B_cd" :: "Volume average of norm B in density interval (dens_c,dens_d]" +"volume_average_norm_B_posz" :: "Volume average of norm B where z > CoM z" +"volume_average_norm_B_negz" :: "Volume average of norm B where z < CoM z" +"volume_average_norm_B_m1_cd" :: "Volume average of m=1 mode of tor., pol. B in (dens_c,dens_d]" +"volume_average_norm_B_m2_cd" :: "Volume average of m=2 mode of tor., pol. B in (dens_c,dens_d]" +"volume_average_norm_B_m3_cd" :: "Volume average of m=3 mode of tor., pol. B in (dens_c,dens_d]" +"volume_average_norm_B_m4_cd" :: "Volume average of m=4 mode of tor., pol. B in (dens_c,dens_d]" "em_energy_ab" :: "Total energy in electro-magnetic sector in density interval (dens_a,dens_b]" "em_energy_cd" :: "Total energy in electro-magnetic sector in density interval (dens_c,dens_d]" +"em_energy_posz" :: "Total energy in electro-magnetic sector where z > CoM z" +"em_energy_negz" :: "Total energy in electro-magnetic sector where z < CoM z" "coordvolume" :: "Coordinate volume of region with rho > dens_atmo" } "nothing" diff --git a/VolumeIntegrals_GRMHDX/schedule.ccl b/VolumeIntegrals_GRMHDX/schedule.ccl index 5ece38b..8d4ec23 100644 --- a/VolumeIntegrals_GRMHDX/schedule.ccl +++ b/VolumeIntegrals_GRMHDX/schedule.ccl @@ -21,11 +21,11 @@ SCHEDULE VI_GRMHDX_file_output_routine_Startup AT CCTK_POST_RECOVER_VARIABLES } ############################# -SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_INITIAL +SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_INITIAL BEFORE AsterX_InitialGroup { LANG: C OPTIONS: GLOBAL - WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) + WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) WRITES: VolIntegral(everywhere) WRITES: volintegral_inside_sphere__center_x(everywhere) WRITES: volintegral_inside_sphere__center_y(everywhere) @@ -39,7 +39,7 @@ SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_POST_RECOVER_VARIABLE { LANG: C OPTIONS: GLOBAL - WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) + WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) WRITES: VolIntegral(everywhere) WRITES: volintegral_inside_sphere__center_x(everywhere) WRITES: volintegral_inside_sphere__center_y(everywhere) @@ -57,7 +57,7 @@ SCHEDULE VI_GRMHDX_InitializeIntegralCounter before VI_GRMHDX_VolumeIntegralGrou WRITES: IntegralCounter(everywhere) } "Initialize IntegralCounter variable" ################## -SCHEDULE GROUP VI_GRMHDX_VolumeIntegralGroup AT poststep BEFORE EstimateError WHILE VolumeIntegrals_GRMHDX::IntegralCounter +SCHEDULE GROUP VI_GRMHDX_VolumeIntegralGroup AT poststep BEFORE (EstimateError Seed_During_Evolution) WHILE VolumeIntegrals_GRMHDX::IntegralCounter { } "Evaluate all volume integrals" @@ -105,7 +105,7 @@ SCHEDULE VI_GRMHDX_DoSum in VI_GRMHDX_VolumeIntegralGroup AFTER VI_GRMHDX_Comput READS: BoxInBox::position_x(everywhere) READS: BoxInBox::position_y(everywhere) READS: BoxInBox::position_z(everywhere) - WRITES: VolIntegral(everywhere) GlobalCoM(everywhere) + WRITES: VolIntegral(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) WRITES: BoxInBox::position_x(everywhere) WRITES: BoxInBox::position_y(everywhere) WRITES: BoxInBox::position_z(everywhere) diff --git a/VolumeIntegrals_GRMHDX/src/Initialize_or_Decrement_IntegralCounter.cxx b/VolumeIntegrals_GRMHDX/src/Initialize_or_Decrement_IntegralCounter.cxx index cdffe56..1797b22 100644 --- a/VolumeIntegrals_GRMHDX/src/Initialize_or_Decrement_IntegralCounter.cxx +++ b/VolumeIntegrals_GRMHDX/src/Initialize_or_Decrement_IntegralCounter.cxx @@ -17,6 +17,12 @@ extern "C" void VI_GRMHDX_InitializeIntegralCounterToZero(CCTK_ARGUMENTS) { *comx = 0; *comy = 0; *comz = 0; + *comx1 = 0; + *comy1 = 0; + *comz1 = 0; + *comx2 = 0; + *comy2 = 0; + *comz2 = 0; // Init Array that Holds Results for (int ii = 0; ii < 101; ii++) { diff --git a/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx b/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx index bb05b20..eb6955e 100644 --- a/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx +++ b/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx @@ -213,6 +213,50 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) { velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y); }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "em_energy_posz")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + double cms_z = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + cms_z = *comz; + } + + magnetic_tot_zsplit(VolIntegrand1, VolIntegrand2, p, + velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz, + gyy, gyz, gzz, cms_x, cms_y, cms_z, 1); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "em_energy_negz")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + double cms_z = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + cms_z = *comz; + } + + magnetic_tot_zsplit(VolIntegrand1, VolIntegrand2, p, + velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz, + gyy, gyz, gzz, cms_x, cms_y, cms_z, -1); + }); } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "volume_average_norm_B_ab")) { grid.loop_all_device<1, 1, 1>( @@ -254,12 +298,128 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) { Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y); }); } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], - "magnetic_energy_comov")) { + "volume_average_norm_B_posz")) { grid.loop_all_device<1, 1, 1>( grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - magnetic_co(VolIntegrand1, p, b2small, w_lorentz, gxx, gxy, gxz, gyy, - gyz, gzz, gamma_lim); + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + double cms_z = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + cms_z = *comz; + } + + volume_norm_B_zsplit(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, cms_z, 1); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_negz")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + double cms_z = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + cms_z = *comz; + } + + volume_norm_B_zsplit(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, cms_z, -1); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m1_cd")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + } + + volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 1); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m2_cd")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + } + + volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 2); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m3_cd")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + } + + volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 3); + }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m4_cd")) { + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // TODO: Get BNS COM from --BNSTrackerGen-- somewhere + // Make sure we have updated values + // before this integral is called! + + double cms_x = 0.0; + double cms_y = 0.0; + if (set_origin_with_VIX) { + cms_x = *comx; + cms_y = *comy; + } + + volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3, + VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, + Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 4); }); } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "magnetic_energy_comov")) { diff --git a/VolumeIntegrals_GRMHDX/src/integrands.cxx b/VolumeIntegrals_GRMHDX/src/integrands.cxx index e437e44..0b3feec 100644 --- a/VolumeIntegrals_GRMHDX/src/integrands.cxx +++ b/VolumeIntegrals_GRMHDX/src/integrands.cxx @@ -307,6 +307,187 @@ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void volume_norm_B_12( } } +/* Volume averaged norm of B in +/- z domain */ +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void volume_norm_B_zsplit( + const Loop::GF3D2 VolIntegrand1, + const Loop::GF3D2 VolIntegrand2, + const Loop::GF3D2 VolIntegrand3, + const Loop::GF3D2 VolIntegrand4, const Loop::PointDesc &p, + const Loop::GF3D2 Bvecx, + const Loop::GF3D2 Bvecy, + const Loop::GF3D2 Bvecz, const Loop::GF3D2 gxx, + const Loop::GF3D2 gxy, const Loop::GF3D2 gxz, + const Loop::GF3D2 gyy, const Loop::GF3D2 gyz, + const Loop::GF3D2 gzz, const double cms_x, + const double cms_y, const double cms_z, const int which_z) { + + if (which_z * (p.z - cms_z) > 0.0) { + const CCTK_REAL gammaDD00 = gxx(p.I); + const CCTK_REAL gammaDD01 = gxy(p.I); + const CCTK_REAL gammaDD02 = gxz(p.I); + const CCTK_REAL gammaDD11 = gyy(p.I); + const CCTK_REAL gammaDD12 = gyz(p.I); + const CCTK_REAL gammaDD22 = gzz(p.I); + + const CCTK_REAL B_contra[3]{Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)}; + + CCTK_REAL posx = p.x - cms_x; + CCTK_REAL posy = p.y - cms_y; + CCTK_REAL posz = p.z; + + // Cylindrical coordinates and some + // quantities for transformation of vectors. + + const double posr2 = posx * posx + posy * posy; + const double posr = sqrt(posr2); + + double xy_over_r2 = posx * posy / posr2; + double x2_over_r2 = posx * posx / posr2; + double y2_over_r2 = posy * posy / posr2; + + if (posr <= 1e-15) { + + xy_over_r2 = 0.0; + x2_over_r2 = 0.5; + y2_over_r2 = 0.5; + } + + const double sqrtgamma = compute_sqrtgamma(p, gxx, gxy, gxz, gyy, gyz, gzz); + + // Covariant B + + const CCTK_REAL B_cov[3]{gammaDD00 * B_contra[0] + gammaDD01 * B_contra[1] + + gammaDD02 * B_contra[2], + gammaDD01 * B_contra[0] + gammaDD11 * B_contra[1] + + gammaDD12 * B_contra[2], + gammaDD02 * B_contra[0] + gammaDD12 * B_contra[1] + + gammaDD22 * B_contra[2]}; + + const CCTK_REAL B_sq = B_cov[0] * B_contra[0] + B_cov[1] * B_contra[1] + + B_cov[2] * B_contra[2]; + + // Finally, compute B_phi*B^phi + + const CCTK_REAL B_2_tor = B_cov[0] * B_contra[0] * y2_over_r2 - + B_cov[0] * B_contra[1] * xy_over_r2 - + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * x2_over_r2; + + // Finally, B_r*B^r + B_z*B^z + + const CCTK_REAL B_2_pol = B_cov[0] * B_contra[0] * x2_over_r2 + + B_cov[0] * B_contra[1] * xy_over_r2 + + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * y2_over_r2 + + B_cov[2] * B_contra[2]; + + VolIntegrand1(p.I) = sqrtgamma * sqrt(abs(B_sq)); + VolIntegrand2(p.I) = sqrtgamma * sqrt(abs(B_2_tor)); + VolIntegrand3(p.I) = sqrtgamma * sqrt(abs(B_2_pol)); + VolIntegrand4(p.I) = sqrtgamma; + + } else { + + VolIntegrand1(p.I) = 0.0; + VolIntegrand2(p.I) = 0.0; + VolIntegrand3(p.I) = 0.0; + VolIntegrand4(p.I) = 0.0; + } +} + +/* Volume averaged mth mode of toroidal, poloidal B in region (dens_1,dens_2]: */ +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void volume_norm_B_m_12( + const Loop::GF3D2 VolIntegrand1, + const Loop::GF3D2 VolIntegrand2, + const Loop::GF3D2 VolIntegrand3, + const Loop::GF3D2 VolIntegrand4, const Loop::PointDesc &p, + const Loop::GF3D2 rho0, const CCTK_REAL dens_1, + const CCTK_REAL dens_2, const Loop::GF3D2 Bvecx, + const Loop::GF3D2 Bvecy, + const Loop::GF3D2 Bvecz, const Loop::GF3D2 gxx, + const Loop::GF3D2 gxy, const Loop::GF3D2 gxz, + const Loop::GF3D2 gyy, const Loop::GF3D2 gyz, + const Loop::GF3D2 gzz, const double cms_x, + const double cms_y, const int mm) { + + const CCTK_REAL my_rho = rho0(p.I); + + CCTK_REAL posx = p.x - cms_x; + CCTK_REAL posy = p.y - cms_y; + CCTK_REAL posz = p.z; + const CCTK_REAL rcyl = posx*posx + posy*posy; // ignore pts close to axis as m>0 modes are discontinuous + + if (my_rho > dens_1 && my_rho <= dens_2 && rcyl > 0.1) { + + const CCTK_REAL gammaDD00 = gxx(p.I); + const CCTK_REAL gammaDD01 = gxy(p.I); + const CCTK_REAL gammaDD02 = gxz(p.I); + const CCTK_REAL gammaDD11 = gyy(p.I); + const CCTK_REAL gammaDD12 = gyz(p.I); + const CCTK_REAL gammaDD22 = gzz(p.I); + + const CCTK_REAL B_contra[3]{Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)}; + + // Cylindrical coordinates and some + // quantities for transformation of vectors. + + const double posr2 = posx * posx + posy * posy; + const double posr = sqrt(posr2); + + double xy_over_r2 = posx * posy / posr2; + double x2_over_r2 = posx * posx / posr2; + double y2_over_r2 = posy * posy / posr2; + + if (posr <= 1e-15) { + + xy_over_r2 = 0.0; + x2_over_r2 = 0.5; + y2_over_r2 = 0.5; + } + + const double sqrtgamma = compute_sqrtgamma(p, gxx, gxy, gxz, gyy, gyz, gzz); + + // Covariant B + + const CCTK_REAL B_cov[3]{gammaDD00 * B_contra[0] + gammaDD01 * B_contra[1] + + gammaDD02 * B_contra[2], + gammaDD01 * B_contra[0] + gammaDD11 * B_contra[1] + + gammaDD12 * B_contra[2], + gammaDD02 * B_contra[0] + gammaDD12 * B_contra[1] + + gammaDD22 * B_contra[2]}; + + const CCTK_REAL B_sq = B_cov[0] * B_contra[0] + B_cov[1] * B_contra[1] + + B_cov[2] * B_contra[2]; + + // Finally, compute B_phi*B^phi + + const CCTK_REAL B_2_tor = B_cov[0] * B_contra[0] * y2_over_r2 - + B_cov[0] * B_contra[1] * xy_over_r2 - + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * x2_over_r2; + + // Finally, B_r*B^r + B_z*B^z + + const CCTK_REAL B_2_pol = B_cov[0] * B_contra[0] * x2_over_r2 + + B_cov[0] * B_contra[1] * xy_over_r2 + + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * y2_over_r2 + + B_cov[2] * B_contra[2]; + + VolIntegrand1(p.I) = sqrtgamma * sqrt(abs(B_2_tor)) * sin(mm * atan2(posy, posx)); + VolIntegrand2(p.I) = sqrtgamma * sqrt(abs(B_2_tor)) * cos(mm * atan2(posy, posx));; + VolIntegrand3(p.I) = sqrtgamma * sqrt(abs(B_2_pol)) * sin(mm * atan2(posy, posx)); + VolIntegrand4(p.I) = sqrtgamma * sqrt(abs(B_2_pol)) * cos(mm * atan2(posy, posx));; + + } else { + + VolIntegrand1(p.I) = 0.0; + VolIntegrand2(p.I) = 0.0; + VolIntegrand3(p.I) = 0.0; + VolIntegrand4(p.I) = 0.0; + } +} + /* Kinetic energy as given by https://arxiv.org/pdf/1206.5911.pdf */ /* See also https://arxiv.org/pdf/2102.01346.pdf for rotational energy */ /* See also https://arxiv.org/abs/astro-ph/0605331 */ @@ -899,6 +1080,139 @@ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void magnetic_tot( 2.0; } +/* Magnetic energy in +/- z domain: */ +CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void magnetic_tot_zsplit( + const Loop::GF3D2 VolIntegrand1, + const Loop::GF3D2 VolIntegrand2, const Loop::PointDesc &p, + const Loop::GF3D2 velx, const Loop::GF3D2 vely, + const Loop::GF3D2 velz, const Loop::GF3D2 Bvecx, + const Loop::GF3D2 Bvecy, + const Loop::GF3D2 Bvecz, const Loop::GF3D2 gxx, + const Loop::GF3D2 gxy, const Loop::GF3D2 gxz, + const Loop::GF3D2 gyy, const Loop::GF3D2 gyz, + const Loop::GF3D2 gzz, const double cms_x, + const double cms_y, const double cms_z, const int which_z) { + + if (which_z * (p.z - cms_z) > 0.0) { + // Notice that CoM_integrand_GAMMA_SPEED_LIMIT is applied to the integrands + // above involving the Lorentz factor w_lorentz. The be consistent we probably + // would have to limit the velocity vel, too. + // This is not done so far. + + const CCTK_REAL gammaDD00 = gxx(p.I); + const CCTK_REAL gammaDD01 = gxy(p.I); + const CCTK_REAL gammaDD02 = gxz(p.I); + const CCTK_REAL gammaDD11 = gyy(p.I); + const CCTK_REAL gammaDD12 = gyz(p.I); + const CCTK_REAL gammaDD22 = gzz(p.I); + + const CCTK_REAL B_contra[3]{Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)}; + + const CCTK_REAL vx = velx(p.I); + const CCTK_REAL vy = vely(p.I); + const CCTK_REAL vz = velz(p.I); + + CCTK_REAL posx = p.x - cms_x; + CCTK_REAL posy = p.y - cms_y; + CCTK_REAL posz = p.z; + + // Cylindrical coordinates and some + // quantities for transformation of vectors. + + const double posr2 = posx * posx + posy * posy; + const double posr = sqrt(posr2); + + double xy_over_r2 = posx * posy / posr2; + double x2_over_r2 = posx * posx / posr2; + double y2_over_r2 = posy * posy / posr2; + + if (posr <= 1e-15) { + + xy_over_r2 = 0.0; + x2_over_r2 = 0.5; + y2_over_r2 = 0.5; + } + + const double sqrtgamma = compute_sqrtgamma(p, gxx, gxy, gxz, gyy, gyz, gzz); + const double gamma = sqrtgamma * sqrtgamma; + + // Covariant quantities + + const CCTK_REAL v_cov[3]{gammaDD00 * vx + gammaDD01 * vy + gammaDD02 * vz, + gammaDD01 * vx + gammaDD11 * vy + gammaDD12 * vz, + gammaDD02 * vx + gammaDD12 * vy + gammaDD22 * vz}; + + const CCTK_REAL B_cov[3]{gammaDD00 * B_contra[0] + gammaDD01 * B_contra[1] + + gammaDD02 * B_contra[2], + gammaDD01 * B_contra[0] + gammaDD11 * B_contra[1] + + gammaDD12 * B_contra[2], + gammaDD02 * B_contra[0] + gammaDD12 * B_contra[1] + + gammaDD22 * B_contra[2]}; + + // Contravariant electric field * sqrtgamma + + const CCTK_REAL sqrtgammaE_contra[3]{ + B_cov[1] * v_cov[2] - B_cov[2] * v_cov[1], + B_cov[2] * v_cov[0] - B_cov[0] * v_cov[2], + B_cov[0] * v_cov[1] - B_cov[1] * v_cov[0]}; + + // Covariant electric field * sqrtgamma + + const CCTK_REAL sqrtgammaE_cov[3]{ + gammaDD00 * sqrtgammaE_contra[0] + gammaDD01 * sqrtgammaE_contra[1] + + gammaDD02 * sqrtgammaE_contra[2], + gammaDD01 * sqrtgammaE_contra[0] + gammaDD11 * sqrtgammaE_contra[1] + + gammaDD12 * sqrtgammaE_contra[2], + gammaDD02 * sqrtgammaE_contra[0] + gammaDD12 * sqrtgammaE_contra[1] + + gammaDD22 * sqrtgammaE_contra[2]}; + + // Finally, compute E_phi*E^phi*sqrtgamma and B_phi*B_phi*sqrtgamma + + CCTK_REAL E_2{0.0}; + CCTK_REAL B_2{0.0}; + + E_2 = sqrtgammaE_cov[0] * sqrtgammaE_contra[0] * y2_over_r2 - + sqrtgammaE_cov[0] * sqrtgammaE_contra[1] * xy_over_r2 - + sqrtgammaE_cov[1] * sqrtgammaE_contra[0] * xy_over_r2 + + sqrtgammaE_cov[1] * sqrtgammaE_contra[1] * x2_over_r2; + E_2 /= sqrtgamma; + + B_2 = B_cov[0] * B_contra[0] * y2_over_r2 - + B_cov[0] * B_contra[1] * xy_over_r2 - + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * x2_over_r2; + B_2 *= sqrtgamma; + + VolIntegrand1(p.I) = (E_2 + B_2) / 2.0; + + // Finally, compute E_r*E^r*sqrtgamma and B_r*B_r*sqrtgamma + + E_2 = 0.0; + B_2 = 0.0; + + E_2 = sqrtgammaE_cov[0] * sqrtgammaE_contra[0] * x2_over_r2 + + sqrtgammaE_cov[0] * sqrtgammaE_contra[1] * xy_over_r2 + + sqrtgammaE_cov[1] * sqrtgammaE_contra[0] * xy_over_r2 + + sqrtgammaE_cov[1] * sqrtgammaE_contra[1] * y2_over_r2; + E_2 /= sqrtgamma; + + B_2 = B_cov[0] * B_contra[0] * x2_over_r2 + + B_cov[0] * B_contra[1] * xy_over_r2 + + B_cov[1] * B_contra[0] * xy_over_r2 + + B_cov[1] * B_contra[1] * y2_over_r2; + B_2 *= sqrtgamma; + + VolIntegrand2(p.I) = + (E_2 + sqrtgammaE_cov[2] * sqrtgammaE_contra[2] / sqrtgamma + B_2 + + B_cov[2] * B_contra[2] * sqrtgamma) / + 2.0; + } else { + + VolIntegrand1(p.I) = 0.0; + VolIntegrand2(p.I) = 0.0; + } +} + /* Magnetic energy split into azimuthal energy component and rest: */ /* Only for specific density region: from dens_1 to dens_2 */ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void magnetic_tot_12( diff --git a/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx b/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx index deac760..00b3a6b 100644 --- a/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx +++ b/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx @@ -33,6 +33,24 @@ VI_GRMHDX_number_of_reductions(int which_integral) { if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "volume_average_norm_B_cd")) return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_posz")) + return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_negz")) + return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m1_cd")) + return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m2_cd")) + return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m3_cd")) + return 4; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_m4_cd")) + return 4; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "kinetic_energy")) return 3; @@ -58,6 +76,10 @@ VI_GRMHDX_number_of_reductions(int which_integral) { return 2; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "em_energy_cd")) return 2; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "em_energy_posz")) + return 2; + if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "em_energy_negz")) + return 2; if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], "coordvolume")) return 1; diff --git a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx index 1946a35..5d4d73a 100644 --- a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx +++ b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx @@ -25,25 +25,19 @@ extern "C" void VI_GRMHDX_DoSum(CCTK_ARGUMENTS) { /* FIXME: Add this symmetry stuff... Should be straightforward. */ CCTK_REAL sym_factor1, sym_factor2, sym_factor3; - /* - if (CCTK_EQUALS(domain,"bitant")){ - sym_factor1 = 2.0e0; - sym_factor2 = 2.0e0; - sym_factor3 = 0.0e0; - } else if (CCTK_EQUALS(domain,"octant")){ - sym_factor1 = 8.0e0; - sym_factor2 = 0.0e0; - sym_factor3 = 0.0e0; - } else { - sym_factor1 = 1.0e0; - sym_factor2 = 1.0e0; - sym_factor3 = 1.0e0; - } - */ - - sym_factor1 = 1.0e0; - sym_factor2 = 1.0e0; - sym_factor3 = 1.0e0; + if (CCTK_EQUALS(VIdomain,"bitant")){ + sym_factor1 = 2.0e0; + sym_factor2 = 2.0e0; + sym_factor3 = 0.0e0; + } else if (CCTK_EQUALS(VIdomain,"octant")){ + sym_factor1 = 8.0e0; + sym_factor2 = 0.0e0; + sym_factor3 = 0.0e0; + } else { + sym_factor1 = 1.0e0; + sym_factor2 = 1.0e0; + sym_factor3 = 1.0e0; + } // CarpetX reductions currently include the cell volume in the sum. Keeping // this code here in case that ever changes. const amrex::Geometry &geom = @@ -182,4 +176,32 @@ extern "C" void VI_GRMHDX_DoSum(CCTK_ARGUMENTS) { which_integral, norm, VolIntegral[4 * (which_integral) + 3]); } } + /* Set CoM tracker for NS 1 */ + if (which_integral == 2 && + CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "centerofmass")) { + const double norm = sym_factor1 * VolIntegral[4 * (which_integral) + 3]; + if (std::isfinite(norm) && std::abs(norm) > 0.0) { + *comx1 = sym_factor2 * VolIntegral[4 * (which_integral) + 0] / norm; + *comy1 = sym_factor2 * VolIntegral[4 * (which_integral) + 1] / norm; + *comz1 = sym_factor3 * VolIntegral[4 * (which_integral) + 2] / norm; + } else if (myproc == 0 && verbose >= 1) { + printf("VolumeIntegrals_GRMHDX: DoSum skipped global COM update: integral=%d norm=%e raw_mass=%e\n", + which_integral, norm, VolIntegral[4 * (which_integral) + 3]); + } + } + /* Set CoM tracker for NS 2 */ + if (which_integral == 3 && + CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "centerofmass")) { + const double norm = sym_factor1 * VolIntegral[4 * (which_integral) + 3]; + if (std::isfinite(norm) && std::abs(norm) > 0.0) { + *comx2 = sym_factor2 * VolIntegral[4 * (which_integral) + 0] / norm; + *comy2 = sym_factor2 * VolIntegral[4 * (which_integral) + 1] / norm; + *comz2 = sym_factor3 * VolIntegral[4 * (which_integral) + 2] / norm; + } else if (myproc == 0 && verbose >= 1) { + printf("VolumeIntegrals_GRMHDX: DoSum skipped global COM update: integral=%d norm=%e raw_mass=%e\n", + which_integral, norm, VolIntegral[4 * (which_integral) + 3]); + } + } }