From 0019f0d4386b4a67aeac18823f78acd98047f054 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Tue, 22 Jul 2025 16:12:24 -0500 Subject: [PATCH 01/11] RNSReader: properly set prims in corotating atmosphere --- RNSReader/param.ccl | 22 +++++++++++----------- RNSReader/schedule.ccl | 4 ++-- RNSReader/src/rnsid_rfr.cxx | 13 ++++++++++++- RNSReader/src/set_vatmo.cxx | 11 +++++++++++ 4 files changed, 36 insertions(+), 14 deletions(-) diff --git a/RNSReader/param.ccl b/RNSReader/param.ccl index 4095c76..53dc7b4 100644 --- a/RNSReader/param.ccl +++ b/RNSReader/param.ccl @@ -13,8 +13,8 @@ USES CCTK_REAL n_rho_atmo private: -# NOTE: These params are used to generate an RNS solution and are thus disabled -# for this stripped-down reader. CarpetX doesn't like treating parameters as non-consts? +# NOTE: These params are used to generate an RNS solution. They are disabled +# for this stripped-down reader and read from the RNS file instead. # KEYWORD rotation_type "Specify type of rotation law" # { @@ -32,16 +32,16 @@ private: # 0: :: "Any positive number" # } 1 -# KEYWORD eos_type "Specify type of equation of state" -# { -# "poly" :: "Polytropic EOS" -# "tab" :: "Tabulated EOS" -# } "poly" +KEYWORD load_eos_type "Specify type of equation of state" +{ + "poly" :: "Polytropic EOS" + "tab" :: "Tabulated EOS" +} "poly" -# STRING eos_file "Equation of state table" -# { -# .* :: "EOS table file" -# } "" +STRING load_eos_file "Equation of state table" +{ + .* :: "EOS table file" +} "" ########## ----------------------------------------------------------------- diff --git a/RNSReader/schedule.ccl b/RNSReader/schedule.ccl index 2bc9e39..9caec9e 100644 --- a/RNSReader/schedule.ccl +++ b/RNSReader/schedule.ccl @@ -60,9 +60,9 @@ if (corot_atmo) { SCHEDULE RNSReader_Init_VelAtmo IN HydroBaseX_InitialData AFTER Hydro_rnsid_init { LANG: C - READS: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) + READS: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) READS: ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere) - WRITES: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) + WRITES: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) } "Set corotating atmosphere from Omega(costh) of RNS surface" } diff --git a/RNSReader/src/rnsid_rfr.cxx b/RNSReader/src/rnsid_rfr.cxx index 1d7a375..4e3bdde 100644 --- a/RNSReader/src/rnsid_rfr.cxx +++ b/RNSReader/src/rnsid_rfr.cxx @@ -99,6 +99,9 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { rotation_type = (char *)malloc(200 * sizeof(char)); double A_diff, axes_ratio; + strcpy(eos_type, load_eos_type); + strcpy(eos_file, load_eos_file); + /* INITIAL DATA VARIABLES */ CCTK_INT m, /* counter */ @@ -287,7 +290,7 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { /* RECOVER FROM 2D FILE */ /* ==================================================== */ message = (char *)malloc(200 * sizeof(char)); - sprintf(message, " Recovering 2D model form file %s", model2D_file); + sprintf(message, " Recovering 2D model from file %s", model2D_file); CCTK_INFO(message); free(message); /* ==================================================== */ @@ -574,6 +577,10 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { /* *************************************** */ /* ATMOSPERE SETTINGs */ /* *************************************** */ + // if (rho(p.I) > 1e-10) { + // printf("Post RNSRead: rho = %e; eps = %e; press = %e;\n", rho(p.I), eps(p.I), press(p.I)); + // } + }); /* END GRID LOOP */ if (strcmp(zero_shift, "yes") == 0) { @@ -610,6 +617,10 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { free(s_gp); free(mu); + free(eos_type); + free(eos_file); + free(rotation_type); + /* Add code for filling time levels here if needed in the future. */ return; diff --git a/RNSReader/src/set_vatmo.cxx b/RNSReader/src/set_vatmo.cxx index b48cafb..8c9f3c6 100644 --- a/RNSReader/src/set_vatmo.cxx +++ b/RNSReader/src/set_vatmo.cxx @@ -56,6 +56,8 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL velxL = velx(p.I); const CCTK_REAL velyL = vely(p.I); const CCTK_REAL rhoL = rho(p.I); + const CCTK_REAL epsL = eps(p.I); + const CCTK_REAL pressL = press(p.I); const CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); // Grading rho @@ -79,6 +81,11 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) : 10 * rho_abs_min; + // TODO: This assumes simple polytrope! properly use EOS calls pls + assert(load_eos_type == "poly"); + CCTK_REAL eps_corot_atm = rho_corot_atm * 100; + CCTK_REAL press_corot_atm = rho_corot_atm * eps_corot_atm; + double alpL = calc_avg_v2c(alp, p); const vec betas_avg( @@ -87,11 +94,15 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { velx(p.I) = (-p.y * omg_atm + betas_avg(0)) / alpL; vely(p.I) = (p.x * omg_atm + betas_avg(1)) / alpL; rho(p.I) = rho_corot_atm; + eps(p.I) = eps_corot_atm; + press(p.I) = press_corot_atm; } else { velx(p.I) = velxL; vely(p.I) = velyL; rho(p.I) = rhoL; + eps(p.I) = epsL; + press(p.I) = pressL; } }); } From 379b105039ade0b1c8553e4bbd3d126b32da1766 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Wed, 10 Sep 2025 09:38:40 -0500 Subject: [PATCH 02/11] VolumeIntegrals_GRMHDX: enable symmetries --- .../src/perform_integration_sum_global.cxx | 32 ++++++++----------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx index c6f4d2e..90e0fbf 100644 --- a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx +++ b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx @@ -20,25 +20,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 = From 2bfc72e71d1bf5c4791101c5eefc4ecd5304d12b Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Tue, 16 Sep 2025 18:18:28 -0500 Subject: [PATCH 03/11] VolumeIntegrals_GRMHDX: add volume averages of m=1-4 modes of B --- VolumeIntegrals_GRMHDX/param.ccl | 4 + .../src/evaluate_integrands_local.cxx | 80 ++++++++++++++++ VolumeIntegrals_GRMHDX/src/integrands.cxx | 92 +++++++++++++++++++ .../src/number_of_reductions.cxx | 12 +++ 4 files changed, 188 insertions(+) diff --git a/VolumeIntegrals_GRMHDX/param.ccl b/VolumeIntegrals_GRMHDX/param.ccl index b5c7079..9bbdf34 100644 --- a/VolumeIntegrals_GRMHDX/param.ccl +++ b/VolumeIntegrals_GRMHDX/param.ccl @@ -86,6 +86,10 @@ 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_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]" "coordvolume" :: "Coordinate volume of region with rho > dens_atmo" diff --git a/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx b/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx index 85cf4f4..51eb4c2 100644 --- a/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx +++ b/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx @@ -249,6 +249,86 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) { VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y); }); + } 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")) { grid.loop_all_device<1, 1, 1>( diff --git a/VolumeIntegrals_GRMHDX/src/integrands.cxx b/VolumeIntegrals_GRMHDX/src/integrands.cxx index a7966c8..c133b0f 100644 --- a/VolumeIntegrals_GRMHDX/src/integrands.cxx +++ b/VolumeIntegrals_GRMHDX/src/integrands.cxx @@ -299,6 +299,98 @@ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void volume_norm_B_12( } } +/* 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); + + if (my_rho > dens_1 && my_rho <= dens_2) { + + 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_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 */ diff --git a/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx b/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx index deac760..2374b2e 100644 --- a/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx +++ b/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx @@ -33,6 +33,18 @@ 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_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; From 57dfa64b2f3afcd279bb7c2a931cd4b2c515dfd3 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Wed, 21 Jan 2026 17:06:11 -0600 Subject: [PATCH 04/11] NRPyPlusTOVID: port to CarpetX --- MHD_Diagnostics/param.ccl | 3 + MHD_Diagnostics/src/master.cxx | 6 + NRPyPlusTOVID/configuration.ccl | 3 + NRPyPlusTOVID/interface.ccl | 18 ++ NRPyPlusTOVID/param.ccl | 66 +++++ NRPyPlusTOVID/schedule.ccl | 44 +++ NRPyPlusTOVID/src/ADMQuantities.h | 67 +++++ NRPyPlusTOVID/src/HydroQuantities.h | 41 +++ NRPyPlusTOVID/src/InitialData.cxx | 191 +++++++++++++ .../convert_TOV_spacetime_vars_to_ADM_vars.h | 40 +++ .../src/interpolate_TOV_solution_to_point.h | 43 +++ NRPyPlusTOVID/src/make.code.defn | 1 + NRPyPlusTOVID/src/tov_interp.h | 254 ++++++++++++++++++ VolumeIntegrals_GRMHDX/src/integrands.cxx | 11 +- 14 files changed, 783 insertions(+), 5 deletions(-) create mode 100644 NRPyPlusTOVID/configuration.ccl create mode 100644 NRPyPlusTOVID/interface.ccl create mode 100644 NRPyPlusTOVID/param.ccl create mode 100644 NRPyPlusTOVID/schedule.ccl create mode 100644 NRPyPlusTOVID/src/ADMQuantities.h create mode 100644 NRPyPlusTOVID/src/HydroQuantities.h create mode 100644 NRPyPlusTOVID/src/InitialData.cxx create mode 100644 NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h create mode 100644 NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h create mode 100644 NRPyPlusTOVID/src/make.code.defn create mode 100644 NRPyPlusTOVID/src/tov_interp.h 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..7e4d1ef 100644 --- a/MHD_Diagnostics/src/master.cxx +++ b/MHD_Diagnostics/src/master.cxx @@ -390,6 +390,11 @@ extern "C" void compute_mhd_derivs(CCTK_ARGUMENTS) { // Use fourth-order ECHO expressions // TODO: check if face index is correct, how to handle densitized B? + divB(p.I) = calc_fd_forward_midpoint(dBx_stag, p) + + calc_fd_forward_midpoint(dBy_stag, p) + + calc_fd_forward_midpoint(dBz_stag, p); + + /* 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. - @@ -410,6 +415,7 @@ extern "C" void compute_mhd_derivs(CCTK_ARGUMENTS) { divB(p.I) = dxi[1] * (Bx_stagg - Bx_stagg_m1) + dxi[2] * (By_stagg - By_stagg_m1) + dxi[3] * (Bz_stagg - Bz_stagg_m1); + */ // =========================================== diff --git a/NRPyPlusTOVID/configuration.ccl b/NRPyPlusTOVID/configuration.ccl new file mode 100644 index 0000000..b56a07c --- /dev/null +++ b/NRPyPlusTOVID/configuration.ccl @@ -0,0 +1,3 @@ +# Configuration definition for thorn NRPyPlusTOVID + +REQUIRES Loop diff --git a/NRPyPlusTOVID/interface.ccl b/NRPyPlusTOVID/interface.ccl new file mode 100644 index 0000000..1403f7c --- /dev/null +++ b/NRPyPlusTOVID/interface.ccl @@ -0,0 +1,18 @@ +implements: NRPyPlusTOVID +inherits: ADMBaseX HydroBaseX + +USES INCLUDE HEADER: loop.hxx +USES INCLUDE HEADER: loop_device.hxx + +PUBLIC: + +CCTK_REAL metric_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { gxx_cc gxy_cc gxz_cc gyy_cc gyz_cc gzz_cc } "cell-centered ADM 3-metric g_ij" + +CCTK_REAL curv_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { kxx_cc kxy_cc kxz_cc kyy_cc kyz_cc kzz_cc } "cell-centered ADM extrinstic curvature K_ij" + +CCTK_REAL lapse_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { alp_cc } "cell-centered ADM lapse function alpha" + +CCTK_REAL shift_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { betax_cc betay_cc betaz_cc} "cell-centered ADM shift function beta^i" + +CCTK_REAL dtlapse_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtalp_cc } "cell-centered Time derivative of ADM lapse function" +CCTK_REAL dtshift_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtbetax_cc dtbetay_cc dtbetaz_cc} "cell-centered Time derivative of ADM shift function" diff --git a/NRPyPlusTOVID/param.ccl b/NRPyPlusTOVID/param.ccl new file mode 100644 index 0000000..eb68e58 --- /dev/null +++ b/NRPyPlusTOVID/param.ccl @@ -0,0 +1,66 @@ +shares: ADMBaseX + +EXTENDS KEYWORD initial_data +{ + "NRPyPlusTOVID" :: "Initial data from NRPyPlusTOVID solution" +} +EXTENDS KEYWORD initial_lapse +{ + "NRPyPlusTOVID" :: "Initial lapse from NRPyPlusTOVID solution" +} +EXTENDS KEYWORD initial_shift +{ + "NRPyPlusTOVID" :: "Initial shift from NRPyPlusTOVID solution" +} +EXTENDS KEYWORD initial_dtlapse +{ + "NRPyPlusTOVID" :: "Initial dtlapse from NRPyPlusTOVID solution" +} +EXTENDS KEYWORD initial_dtshift +{ + "NRPyPlusTOVID" :: "Initial dtshift from NRPyPlusTOVID solution" +} + +shares: HydroBaseX +EXTENDS KEYWORD initial_hydro +{ + "NRPyPlusTOVID" :: "Initial GRHD data from NRPyPlusTOVID solution" +} + +restricted: +CCTK_STRING TOV_filename "Which interpolator should I use" +{ + ".+" :: "Any nonempty string" +} "outputTOVpolytrope.txt" + +restricted: +CCTK_REAL rho_atmosphere "Atmosphere baryonic density" +{ + 0:* :: "Physical values" +-1 :: "forbidden value to make sure it is explicitly set in the parfile" +} -1 + +restricted: +CCTK_REAL K_atmosphere "Polytropic K to be used with the EOS corresponding to rho_atmosphere" +{ + 0:* :: "Physical values" +-1 :: "forbidden value to make sure it is explicitly set in the parfile" +} -1 + +restricted: +CCTK_REAL Gamma_atmosphere "Polytropic Gamma to be used with the EOS corresponding to rho_atmosphere" +{ + 0:* :: "Physical values" +-1 :: "forbidden value to make sure it is explicitly set in the parfile" +} -1 + +CCTK_REAL pert_amp "velocity perturbation amplitude" +{ + 0:* :: "Physical values" +} 0.0 + +restricted: +CCTK_REAL Pressure_depletion_factor "Pressure depletion factor = Pdf: P => (1-Pdf)*P" +{ + 0:* :: "Greater than or equal to zero, where zero is no depletion and default." +} 0.0 diff --git a/NRPyPlusTOVID/schedule.ccl b/NRPyPlusTOVID/schedule.ccl new file mode 100644 index 0000000..5921e24 --- /dev/null +++ b/NRPyPlusTOVID/schedule.ccl @@ -0,0 +1,44 @@ + +schedule NRPyPlusTOVID_ET_InitialData IN HydroBaseX_InitialData +{ + LANG: C + WRITES: metric_cc(everywhere) + WRITES: lapse_cc(everywhere) + WRITES: shift_cc(everywhere) + # WRITES: curv_cc(everywhere) + # WRITES: dtlapse_cc(everywhere) + # WRITES: dtshift_cc(everywhere) + WRITES: HYDROBASEX::eps(everywhere) + WRITES: HYDROBASEX::rho(everywhere) + WRITES: HYDROBASEX::vel(everywhere) + WRITES: HYDROBASEX::press(everywhere) +} "Set up general relativistic hydrodynamic (GRHD) fields for TOV initial data" + +schedule NRPyPlusTOVID_ET_VelPert IN HydroBaseX_InitialData AFTER NRPyPlusTOVID_ET_InitialData +{ + LANG: C + READS: HYDROBASEX::rho(everywhere) + WRITES: HYDROBASEX::vel(everywhere) +} "Set up general relativistic hydrodynamic (GRHD) fields for TOV initial data" + +SCHEDULE NRPyPlusTOVID_Interpolation_C2V IN HydroBaseX_InitialData AFTER NRPyPlusTOVID_ET_InitialData +{ + LANG: C + READS: metric_cc(everywhere) + READS: lapse_cc(everywhere) + READS: shift_cc(everywhere) + # READS: curv_cc(everywhere) + # READS: dtlapse_cc(everywhere) + # READS: dtshift_cc(everywhere) + + WRITES: ADMBaseX::metric(interior) + WRITES: ADMBaseX::lapse(interior) + WRITES: ADMBaseX::shift(interior) + WRITES: ADMBaseX::curv(interior) + WRITES: ADMBaseX::dtlapse(interior) + WRITES: ADMBaseX::dtshift(interior) + + SYNC: ADMBaseX::metric ADMBaseX::lapse ADMBaseX::shift ADMBaseX::curv + SYNC: ADMBaseX::dtlapse ADMBaseX::dtshift +} "Interpolate metric from cell center to vertex" + diff --git a/NRPyPlusTOVID/src/ADMQuantities.h b/NRPyPlusTOVID/src/ADMQuantities.h new file mode 100644 index 0000000..a01b3ec --- /dev/null +++ b/NRPyPlusTOVID/src/ADMQuantities.h @@ -0,0 +1,67 @@ + +namespace NRPyPlusTOVID{ +using namespace std; +using namespace Loop; + +static inline +void ADMQuantities(const PointDesc &p, + const CCTK_REAL IDalpha, + const CCTK_REAL IDgammaDD00,const CCTK_REAL IDgammaDD01, const CCTK_REAL IDgammaDD02, + const CCTK_REAL IDgammaDD11,const CCTK_REAL IDgammaDD12, const CCTK_REAL IDgammaDD22, + + const GF3D2 &alphaGF, + const GF3D2 &betaU0GF, + const GF3D2 &betaU1GF, + const GF3D2 &betaU2GF, + const GF3D2 &gammaDD00GF, + const GF3D2 &gammaDD01GF, + const GF3D2 &gammaDD02GF, + const GF3D2 &gammaDD11GF, + const GF3D2 &gammaDD12GF, + const GF3D2 &gammaDD22GF) { + const CCTK_REAL xx0 = p.x; + const CCTK_REAL xx1 = p.y; + const CCTK_REAL xx2 = p.z; + + /* + * NRPy+ Finite Difference Code Generation, Step 1 of 1: Evaluate SymPy expressions and write to main memory: + */ + const double FDPart3_0 = ((xx1)*(xx1)); + const double FDPart3_1 = ((xx0)*(xx0)); + const double FDPart3_2 = FDPart3_0 + FDPart3_1; + const double FDPart3_3 = IDgammaDD22/((FDPart3_2)*(FDPart3_2)); + const double FDPart3_4 = ((xx2)*(xx2)); + const double FDPart3_5 = FDPart3_2 + FDPart3_4; + const double FDPart3_6 = (1.0/(FDPart3_5)); + const double FDPart3_7 = FDPart3_6*IDgammaDD00; + const double FDPart3_8 = (1.0/(FDPart3_2)); + const double FDPart3_9 = (1.0/sqrt(FDPart3_5)); + const double FDPart3_11 = FDPart3_8*FDPart3_9*IDgammaDD02; + const double FDPart3_12 = xx0*xx1; + const double FDPart3_13 = 2*FDPart3_12; + const double FDPart3_15 = (1.0/((FDPart3_5)*(FDPart3_5))); + const double FDPart3_17 = -FDPart3_4*FDPart3_6 + 1; + const double FDPart3_18 = (1.0/sqrt(FDPart3_17)); + const double FDPart3_19 = FDPart3_18*IDgammaDD01; + const double FDPart3_20 = 2*FDPart3_19*xx2; + const double FDPart3_22 = IDgammaDD11/FDPart3_17; + const double FDPart3_23 = FDPart3_22*FDPart3_4/((FDPart3_5)*(FDPart3_5)*(FDPart3_5)); + const double FDPart3_25 = FDPart3_18*FDPart3_8*IDgammaDD12; + const double FDPart3_26 = pow(FDPart3_5, -3.0/2.0); + const double FDPart3_27 = FDPart3_26*xx2; + const double FDPart3_29 = FDPart3_13*FDPart3_25*FDPart3_27; + const double FDPart3_35 = -FDPart3_26*FDPart3_4 + FDPart3_9; + const double FDPart3_36 = FDPart3_35*xx1; + alphaGF(p.I) = IDalpha; + betaU0GF(p.I) = 0; + betaU1GF(p.I) = 0; + betaU2GF(p.I) = 0; + gammaDD00GF(p.I) = FDPart3_0*FDPart3_3 + FDPart3_1*FDPart3_15*FDPart3_20 + FDPart3_1*FDPart3_23 + FDPart3_1*FDPart3_7 - FDPart3_11*FDPart3_13 - FDPart3_29; + gammaDD01GF(p.I) = -FDPart3_0*FDPart3_11 - FDPart3_0*FDPart3_25*FDPart3_27 + FDPart3_1*FDPart3_18*FDPart3_27*FDPart3_8*IDgammaDD12 + FDPart3_1*FDPart3_8*FDPart3_9*IDgammaDD02 + FDPart3_12*FDPart3_23 - FDPart3_12*FDPart3_3 + FDPart3_12*FDPart3_7 + FDPart3_13*FDPart3_15*FDPart3_19*xx2; + gammaDD02GF(p.I) = -FDPart3_11*xx1*xx2 + FDPart3_15*FDPart3_19*FDPart3_4*xx0 - FDPart3_19*FDPart3_35*FDPart3_9*xx0 - FDPart3_22*FDPart3_27*FDPart3_35*xx0 + FDPart3_25*FDPart3_36 + FDPart3_7*xx0*xx2; + gammaDD11GF(p.I) = FDPart3_0*FDPart3_15*FDPart3_20 + FDPart3_0*FDPart3_23 + FDPart3_0*FDPart3_7 + FDPart3_1*FDPart3_3 + FDPart3_11*FDPart3_13 + FDPart3_29; + gammaDD12GF(p.I) = FDPart3_11*xx0*xx2 + FDPart3_15*FDPart3_19*FDPart3_4*xx1 - FDPart3_19*FDPart3_36*FDPart3_9 - FDPart3_22*FDPart3_27*FDPart3_36 - FDPart3_25*FDPart3_35*xx0 + FDPart3_7*xx1*xx2; + gammaDD22GF(p.I) = -FDPart3_20*FDPart3_35*FDPart3_9 + FDPart3_22*((FDPart3_35)*(FDPart3_35)) + FDPart3_4*FDPart3_6*IDgammaDD00; +} + +} //namespace diff --git a/NRPyPlusTOVID/src/HydroQuantities.h b/NRPyPlusTOVID/src/HydroQuantities.h new file mode 100644 index 0000000..5e8b1a1 --- /dev/null +++ b/NRPyPlusTOVID/src/HydroQuantities.h @@ -0,0 +1,41 @@ + + +namespace NRPyPlusTOVID{ +using namespace std; +using namespace Loop; + +static inline +void HydroQuantities(const PointDesc &p, + const CCTK_REAL IDPressure, const CCTK_REAL IDrho_baryonic, + const CCTK_REAL IDrho__total_energy_density, + const GF3D2 &PressureGF, + const GF3D2 &rho_baryonicGF, + const GF3D2 &epsilonGF, + const GF3D2 &Valencia3velocityU0GF, + const GF3D2 &Valencia3velocityU1GF, + const GF3D2 &Valencia3velocityU2GF) { + DECLARE_CCTK_PARAMETERS; + if(IDrho__total_energy_density <= 0 || IDrho_baryonic <= 0 || IDPressure <= 0) { + rho_baryonicGF(p.I) = rho_atmosphere; + PressureGF(p.I) = K_atmosphere*pow(rho_atmosphere,Gamma_atmosphere); + epsilonGF(p.I) = 0; + Valencia3velocityU0GF(p.I) = 0; + Valencia3velocityU1GF(p.I) = 0; + Valencia3velocityU2GF(p.I) = 0; + } else { + /* + * NRPy+ Finite Difference Code Generation, Step 1 of 1: Evaluate SymPy expressions and write to main memory: + */ + PressureGF(p.I) = IDPressure; + rho_baryonicGF(p.I) = IDrho_baryonic; + epsilonGF(p.I) = IDrho__total_energy_density/IDrho_baryonic - 1; + Valencia3velocityU0GF(p.I) = 0; + Valencia3velocityU1GF(p.I) = 0; + Valencia3velocityU2GF(p.I) = 0; + + // Apply pressure depletion. + PressureGF(p.I) *= (1.0 - Pressure_depletion_factor); + } +} + +} // namespace diff --git a/NRPyPlusTOVID/src/InitialData.cxx b/NRPyPlusTOVID/src/InitialData.cxx new file mode 100644 index 0000000..7b5e84a --- /dev/null +++ b/NRPyPlusTOVID/src/InitialData.cxx @@ -0,0 +1,191 @@ + +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#include +#include + +#include "aster_utils.hxx" + +#include "ADMQuantities.h" +#include "HydroQuantities.h" +#include "interpolate_TOV_solution_to_point.h" +#include "convert_TOV_spacetime_vars_to_ADM_vars.h" + +namespace NRPyPlusTOVID { +using namespace std; +using namespace AsterUtils; +using namespace Loop; + + +// Alias for "vel" vector gridfunction: +// #define velx (&vel[0*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) +// #define vely (&vel[1*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) +// #define velz (&vel[2*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) + +void read_TOV_input_data_from_file(ID_inputs *TOV_in) { + + DECLARE_CCTK_PARAMETERS; + + // Step 1: Set up TOV initial data + // Step 1.a: Read TOV initial data from data file + // Open the data file: + char filename[100]; + sprintf(filename,"%s",TOV_filename); // TOV_filename is a CCTK_PARAMETER + FILE *in1Dpolytrope = fopen(filename, "r"); + if (in1Dpolytrope == NULL) { + fprintf(stderr,"ERROR: could not open file %s\n",filename); + exit(1); + } + // Count the number of lines in the data file: + int numlines_in_file = count_num_lines_in_file(in1Dpolytrope); + // Allocate space for all data arrays: + CCTK_REAL *r_Schw_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *rho_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *rho_baryon_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *P_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *M_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *expnu_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *exp4phi_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + CCTK_REAL *rbar_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); + + // Read from the data file, filling in arrays. + // read_datafile__set_arrays() may be found in TOV/tov_interp.h + if(read_datafile__set_arrays(in1Dpolytrope, r_Schw_arr,rho_arr,rho_baryon_arr,P_arr,M_arr,expnu_arr,exp4phi_arr,rbar_arr) == 1) { + fprintf(stderr,"ERROR WHEN READING FILE %s!\n",filename); + exit(1); + } + fclose(in1Dpolytrope); + REAL Rbar = -100; + int Rbar_idx = -100; + for(int i=1;i0 && rho_arr[i]==0) { Rbar = rbar_arr[i-1]; Rbar_idx = i-1; } + } + if(Rbar<0) { + fprintf(stderr,"Error: could not find rbar=Rbar from data file.\n"); + exit(1); + } + + TOV_in->Rbar = Rbar; + TOV_in->Rbar_idx = Rbar_idx; + + const int interp_stencil_size = 12; + TOV_in->interp_stencil_size = interp_stencil_size; + TOV_in->numlines_in_file = numlines_in_file; + + TOV_in->r_Schw_arr = r_Schw_arr; + TOV_in->rho_arr = rho_arr; + TOV_in->rho_baryon_arr = rho_baryon_arr; + TOV_in->P_arr = P_arr; + TOV_in->M_arr = M_arr; + TOV_in->expnu_arr = expnu_arr; + TOV_in->exp4phi_arr = exp4phi_arr; + TOV_in->rbar_arr = rbar_arr; + /* END TOV INPUT ROUTINE */ +} + +extern "C" void NRPyPlusTOVID_ET_InitialData(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_ET_InitialData; + DECLARE_CCTK_PARAMETERS; + + ID_inputs TOV_in; + read_TOV_input_data_from_file(&TOV_in); + + const array indextype = {1, 1, 1}; + const GF3D2layout CCC_layout(cctkGH, indextype); + + grid.loop_all<1, 1, 1>(grid.nghostzones, + [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_INT idx = CCC_layout.linear(p.I); + CCTK_REAL rr = sqrt(p.x*p.x + p.y*p.y + p.z*p.z); + CCTK_REAL th = acos(p.z/rr); + + CCTK_REAL IDexp_4phi,IDnu,IDPressure,IDrho_baryonic,IDrho__total_energy_density; + interpolate_TOV_solution_to_point(rr, TOV_in, &IDexp_4phi,&IDnu, + &IDPressure,&IDrho_baryonic,&IDrho__total_energy_density); + + CCTK_REAL IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22; + convert_TOV_spacetime_vars_to_ADM_vars(rr, th, IDexp_4phi,IDnu, + &IDalpha,&IDgammaDD00,&IDgammaDD01,&IDgammaDD02,&IDgammaDD11,&IDgammaDD12,&IDgammaDD22); + + HydroQuantities(p, + IDPressure,IDrho_baryonic,IDrho__total_energy_density, + press,rho,eps,velx,vely,velz); + + ADMQuantities(p, + IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22, + alp_cc,betax_cc,betay_cc,betaz_cc, + gxx_cc,gxy_cc,gxz_cc,gyy_cc,gyz_cc,gzz_cc); + }); + + free(TOV_in.r_Schw_arr); + free(TOV_in.rho_arr); + free(TOV_in.rho_baryon_arr); + free(TOV_in.P_arr); + free(TOV_in.M_arr); + free(TOV_in.expnu_arr); + free(TOV_in.exp4phi_arr); + free(TOV_in.rbar_arr); +} + +extern "C" void NRPyPlusTOVID_Interpolation_C2V(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_Interpolation_C2V; + DECLARE_CCTK_PARAMETERS; + + CCTK_INFO("Starting interpolation for ADM variables."); + grid.loop_int<0, 0, 0>(grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + gxx(p.I) = calc_avg_c2v<4>(gxx_cc, p); + gxy(p.I) = calc_avg_c2v<4>(gxy_cc, p); + gxz(p.I) = calc_avg_c2v<4>(gxz_cc, p); + gyy(p.I) = calc_avg_c2v<4>(gyy_cc, p); + gyz(p.I) = calc_avg_c2v<4>(gyz_cc, p); + gzz(p.I) = calc_avg_c2v<4>(gzz_cc, p); + + kxx(p.I) = 0.0; + kxy(p.I) = 0.0; + kxz(p.I) = 0.0; + kyy(p.I) = 0.0; + kyz(p.I) = 0.0; + kzz(p.I) = 0.0; + + alp(p.I) = calc_avg_c2v<4>(alp_cc, p); + betax(p.I) = calc_avg_c2v<4>(betax_cc, p); + betay(p.I) = calc_avg_c2v<4>(betay_cc, p); + betaz(p.I) = calc_avg_c2v<4>(betaz_cc, p); + + dtalp(p.I) = 0.0; + dtbetax(p.I) = 0.0; + dtbetay(p.I) = 0.0; + dtbetaz(p.I) = 0.0; + }); + + CCTK_INFO("Done interpolation for ADM variables."); +} + +extern "C" void NRPyPlusTOVID_ET_VelPert(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_ET_VelPert; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all<1, 1, 1>(grid.nghostzones, + [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + if (rho(p.I) > rho_atmosphere) { + velx(p.I) = pert_amp * p.x; + vely(p.I) = pert_amp * p.y; + velz(p.I) = pert_amp * p.z; + } + else { + velx(p.I) = 0.0; + vely(p.I) = 0.0; + velz(p.I) = 0.0; + } + }); +} + +} // namespace diff --git a/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h b/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h new file mode 100644 index 0000000..ca1d5e4 --- /dev/null +++ b/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h @@ -0,0 +1,40 @@ + +/* This function converts TOV quantities into + * ADM quantities in Spherical coordinates. + */ +void convert_TOV_spacetime_vars_to_ADM_vars( const CCTK_REAL rr, const CCTK_REAL th, + const CCTK_REAL IDexp_4phi, const CCTK_REAL IDexpnu, + CCTK_REAL *IDalpha, + CCTK_REAL *IDgammaDD00, CCTK_REAL *IDgammaDD01, CCTK_REAL *IDgammaDD02, + CCTK_REAL *IDgammaDD11, CCTK_REAL *IDgammaDD12, CCTK_REAL *IDgammaDD22) { + + /*************************************************************** + * Convert TOV quantities to ADM quantities in Spherical basis * + *************************************************************** + * + * First we convert the lapse function: + * .------------------. + * | alpha = e^(nu/2) | + * .------------------. + */ + *IDalpha = sqrt(IDexpnu); + + /* Next we convert the metric function: + * .----------------------------------------. + * | gamma_{00} = e^{4phi} | + * .----------------------------------------. + * | gamma_{11} = e^{4phi} r^2 | + * .----------------------------------------. + * | gamma_{22} = e^{4phi} r^2 sin^2(theta) | + * .----------------------------------------. + * | All other components are zero. | + * .----------------------------------------. + */ + *IDgammaDD00 = IDexp_4phi; + *IDgammaDD11 = IDexp_4phi * rr * rr; + *IDgammaDD22 = IDexp_4phi * rr * rr * sin(th) * sin(th); + *IDgammaDD01 = 0.0; + *IDgammaDD02 = 0.0; + *IDgammaDD12 = 0.0; + +} diff --git a/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h b/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h new file mode 100644 index 0000000..ea98ad6 --- /dev/null +++ b/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h @@ -0,0 +1,43 @@ + +/* Load the TOV_interpolate_1D() function */ +#include "tov_interp.h" + +// Declare initial data input struct: +// stores data from initial data solver, +// so they can be put on the numerical grid. +typedef struct __ID_inputs { + CCTK_REAL Rbar; + int Rbar_idx; + int interp_stencil_size; + int numlines_in_file; + CCTK_REAL *r_Schw_arr,*rho_arr,*rho_baryon_arr,*P_arr,*M_arr,*expnu_arr,*exp4phi_arr,*rbar_arr; +} ID_inputs; + +/* This function returns the TOV quantities at point rr + * by interpolating the data in the TOV initial data file. + */ +void interpolate_TOV_solution_to_point(const CCTK_REAL rr, ID_inputs other_inputs, + CCTK_REAL *exp_4phi, CCTK_REAL *expnu, + CCTK_REAL *Pressure, CCTK_REAL *rho_baryon, CCTK_REAL *rho__total_energy_density) { + + /* The mass valus is not used, but we have to + * store it in this dummy variable because the + * initial data file contains it. + */ + CCTK_REAL M; + + /* Perform the interpolation, returning: + * - rho__total_energy_density + * - rho_baryon + * - Pressure + * - Mass (dummy variable, unused) + * - exp(nu) + * - exp(4phi) + */ + TOV_interpolate_1D(rr,other_inputs.Rbar,other_inputs.Rbar_idx,other_inputs.interp_stencil_size, + other_inputs.numlines_in_file, + other_inputs.r_Schw_arr,other_inputs.rho_arr,other_inputs.rho_baryon_arr,other_inputs.P_arr,other_inputs.M_arr, + other_inputs.expnu_arr,other_inputs.exp4phi_arr,other_inputs.rbar_arr, + rho__total_energy_density,rho_baryon,Pressure,&M,expnu,exp_4phi); + +} diff --git a/NRPyPlusTOVID/src/make.code.defn b/NRPyPlusTOVID/src/make.code.defn new file mode 100644 index 0000000..a8dea9a --- /dev/null +++ b/NRPyPlusTOVID/src/make.code.defn @@ -0,0 +1 @@ +SRCS = InitialData.cxx diff --git a/NRPyPlusTOVID/src/tov_interp.h b/NRPyPlusTOVID/src/tov_interp.h new file mode 100644 index 0000000..e22c0d3 --- /dev/null +++ b/NRPyPlusTOVID/src/tov_interp.h @@ -0,0 +1,254 @@ +// This C header file reads a TOV solution from data file and performs +// 1D interpolation of the solution to a desired radius. + +// Author: Zachariah B. Etienne +// zachetie **at** gmail **dot* com + +#include "stdio.h" +#include "stdlib.h" +#include "math.h" +#include "string.h" + +#define REAL double + +//#define STANDALONE_UNIT_TEST + +int count_num_lines_in_file(FILE *in1Dpolytrope) { + int numlines_in_file = 0; + char * line = NULL; + + size_t len = 0; + ssize_t read; + while ((read = getline(&line, &len, in1Dpolytrope)) != -1) { + numlines_in_file++; + } + rewind(in1Dpolytrope); + + free(line); + return numlines_in_file; +} + +int read_datafile__set_arrays(FILE *in1Dpolytrope, REAL *restrict r_Schw_arr,REAL *restrict rho_arr,REAL *restrict rho_baryon_arr,REAL *restrict P_arr, + REAL *restrict M_arr,REAL *restrict expnu_arr,REAL *restrict exp4phi_arr,REAL *restrict rbar_arr) { + char * line = NULL; + + size_t len = 0; + ssize_t read; + + int which_line = 0; + while ((read = getline(&line, &len, in1Dpolytrope)) != -1) { + // Define the line delimiters (i.e., the stuff that goes between the data on a given + // line of data. Here, we define both spaces " " and tabs "\t" as data delimiters. + const char delimiters[] = " \t"; + + //Now we define "token", a pointer to the first column of data + char *token; + + //Each successive time we call strtok(NULL,blah), we read in a new column of data from + // the originally defined character array, as pointed to by token. + + token=strtok(line, delimiters); if(token==NULL) { printf("BADDDD\n"); return 1; } + r_Schw_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + rho_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + rho_baryon_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + P_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + M_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + expnu_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + exp4phi_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); + rbar_arr[which_line] = strtod(token, NULL); + + which_line++; + } + free(line); + return 0; +} + +// Find interpolation index using Bisection root-finding algorithm: +static inline int bisection_idx_finder(const REAL rrbar, const int numlines_in_file, const REAL *restrict rbar_arr) { + int x1 = 0; + int x2 = numlines_in_file-1; + REAL y1 = rrbar-rbar_arr[x1]; + REAL y2 = rrbar-rbar_arr[x2]; + if(y1*y2 >= 0) { + fprintf(stderr,"INTERPOLATION BRACKETING ERROR %e | %e %e\n",rrbar,y1,y2); + exit(1); + } + for(int i=0;i (B)) ? (A) : (B) ) + + int idxmin = MAX(0,idx-interp_stencil_size/2-1); + +#ifdef MIN +#undef MIN +#endif +#define MIN(A, B) ( ((A) < (B)) ? (A) : (B) ) + + // -= Do not allow the interpolation stencil to cross the star's surface =- + // max index is when idxmin + (interp_stencil_size-1) = Rbar_idx + // -> idxmin at most can be Rbar_idx - interp_stencil_size + 1 + if(rrbar < Rbar) { + idxmin = MIN(idxmin,Rbar_idx - interp_stencil_size + 1); + } else { + idxmin = MAX(idxmin,Rbar_idx+1); + idxmin = MIN(idxmin,numlines_in_file - interp_stencil_size + 1); + } + // Now perform the Lagrange polynomial interpolation: + + // First set the interpolation coefficients: + REAL rbar_sample[interp_stencil_size]; + for(int i=idxmin;i Rbar) { + *rho = 0; + *rho_baryon = 0; + *P = 0; + *M = M_arr[Rbar_idx+1]; + *expnu = 1. - 2.*(*M) / r_Schw; + *exp4phi = pow(r_Schw / rrbar,2.0); + } +} + +// To compile, copy this file to tov_interp.c, and then run: +// gcc -Ofast tov_interp.c -o tov_interp -DSTANDALONE_UNIT_TEST +#ifdef STANDALONE_UNIT_TEST +int main() { + + // Open the data file: + char filename[100]; + sprintf(filename,"../outputTOVpolytrope.txt"); + FILE *in1Dpolytrope = fopen(filename, "r"); + if (in1Dpolytrope == NULL) { + fprintf(stderr,"ERROR: could not open file %s\n",filename); + exit(1); + } + + // Count the number of lines in the data file: + int numlines_in_file = count_num_lines_in_file(in1Dpolytrope); + + // Allocate space for all data arrays: + REAL *r_Schw_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *rho_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *rho_baryon_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *P_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *M_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *expnu_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *exp4phi_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + REAL *rbar_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); + + // Read from the data file, filling in arrays + if(read_datafile__set_arrays(in1Dpolytrope, r_Schw_arr,rho_arr,rho_baryon_arr,P_arr,M_arr,expnu_arr,exp4phi_arr,rbar_arr) == 1) { + fprintf(stderr,"ERROR WHEN READING FILE %s!\n",filename); + exit(1); + } + fclose(in1Dpolytrope); + + REAL Rbar = -100; + int Rbar_idx = -100; + for(int i=1;i0 && rho_arr[i]==0) { Rbar = rbar_arr[i-1]; Rbar_idx = i-1; } + } + if(Rbar<0) { + fprintf(stderr,"Error: could not find r=R from data file.\n"); + exit(1); + } + + // Next, interpolate! + // Create trial radius array: + int num_r_pts = 100000; + //REAL *r_out_arr = (REAL *)malloc(sizeof(REAL)*num_r_pts); + struct drand48_data randBuffer; + srand48_r(1313, &randBuffer); +#pragma omp parallel for + for(int i=0;i dens_1 && my_rho <= dens_2) { + 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); @@ -327,10 +332,6 @@ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void volume_norm_B_m_1 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. From 97776524ec3e3a4b639f64ebad5e30a0bb4a5ca8 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Thu, 29 Jan 2026 23:06:26 -0600 Subject: [PATCH 05/11] RNSReader: hardcode tabulated EOS requirement --- MHD_Diagnostics/src/master.cxx | 6 +++--- RNSReader/param.ccl | 3 +++ RNSReader/schedule.ccl | 7 ++++--- RNSReader/src/set_vatmo.cxx | 18 ++++++++++++++---- 4 files changed, 24 insertions(+), 10 deletions(-) diff --git a/MHD_Diagnostics/src/master.cxx b/MHD_Diagnostics/src/master.cxx index 7e4d1ef..f7cd5aa 100644 --- a/MHD_Diagnostics/src/master.cxx +++ b/MHD_Diagnostics/src/master.cxx @@ -390,9 +390,9 @@ extern "C" void compute_mhd_derivs(CCTK_ARGUMENTS) { // Use fourth-order ECHO expressions // TODO: check if face index is correct, how to handle densitized B? - divB(p.I) = calc_fd_forward_midpoint(dBx_stag, p) + - calc_fd_forward_midpoint(dBy_stag, p) + - calc_fd_forward_midpoint(dBz_stag, p); + 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); /* double Bx_stagg = 13. * dBx_stag(p.I + p.DI[0]) / 12. - diff --git a/RNSReader/param.ccl b/RNSReader/param.ccl index 53dc7b4..03390ec 100644 --- a/RNSReader/param.ccl +++ b/RNSReader/param.ccl @@ -5,6 +5,9 @@ # #USES CCTK_INT timelevels +SHARES:EOSX +USES KEYWORD evolution_eos + SHARES: Con2PrimFactory USES CCTK_REAL rho_abs_min USES CCTK_REAL atmo_tol diff --git a/RNSReader/schedule.ccl b/RNSReader/schedule.ccl index 9caec9e..65856cf 100644 --- a/RNSReader/schedule.ccl +++ b/RNSReader/schedule.ccl @@ -57,12 +57,13 @@ if (CCTK_Equals(initial_data,"RNSReader")) } if (corot_atmo) { - SCHEDULE RNSReader_Init_VelAtmo IN HydroBaseX_InitialData AFTER Hydro_rnsid_init + SCHEDULE RNSReader_Init_VelAtmo IN HydroBaseX_PostInitial AFTER ID_TabEOS_init { LANG: C - READS: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) + READS: HydroBaseX::velx(everywhere) HydroBaseX::vely(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) + READS: HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) READS: ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere) - WRITES: HydroBaseX::vel(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) + WRITES: HydroBaseX::velx(everywhere) HydroBaseX::vely(everywhere) HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) } "Set corotating atmosphere from Omega(costh) of RNS surface" } diff --git a/RNSReader/src/set_vatmo.cxx b/RNSReader/src/set_vatmo.cxx index 8c9f3c6..78a3d57 100644 --- a/RNSReader/src/set_vatmo.cxx +++ b/RNSReader/src/set_vatmo.cxx @@ -5,6 +5,7 @@ #include "rnsreader_utils.hxx" #include "consts.h" #include "aster_utils.hxx" +#include "setup_eos.hxx" #include #include @@ -15,11 +16,19 @@ using namespace Loop; using namespace amrex; using namespace std; using namespace AsterUtils; +using namespace EOSX; extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_RNSReader_Init_VelAtmo; DECLARE_CCTK_PARAMETERS; + auto eos_3p_tab3d = global_eos_3p_tab3d; + if (not CCTK_EQUALS(evolution_eos, "Tabulated3d")) { + CCTK_VERROR("Invalid evolution EOS type '%s'. Please, set " + "EOSX::evolution_eos = \"Tabulated3d\" in your parameter file.", + evolution_eos); + } + Omega_th_reader* vel_th_reader = nullptr; // Open the Omega(th) file FILE *vel_th_file = fopen(vel_th_filename, "r"); @@ -56,6 +65,8 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL velxL = velx(p.I); const CCTK_REAL velyL = vely(p.I); const CCTK_REAL rhoL = rho(p.I); + const CCTK_REAL tempL = temperature(p.I); + const CCTK_REAL YeL = Ye(p.I); const CCTK_REAL epsL = eps(p.I); const CCTK_REAL pressL = press(p.I); const CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); @@ -77,14 +88,13 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { ? (omg_surf * pow((r_corot_min / radial_distance), n_omg_atmo)) : omg_surf; + // Pad density and recalculate primitives CCTK_REAL rho_corot_atm = (radial_distance > r_atmo) ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) : 10 * rho_abs_min; - // TODO: This assumes simple polytrope! properly use EOS calls pls - assert(load_eos_type == "poly"); - CCTK_REAL eps_corot_atm = rho_corot_atm * 100; - CCTK_REAL press_corot_atm = rho_corot_atm * eps_corot_atm; + CCTK_REAL eps_corot_atm = eos_3p_tab3d->eps_from_valid_rho_temp_ye(rho_corot_atm, tempL, YeL); + CCTK_REAL press_corot_atm = eos_3p_tab3d->press_from_valid_rho_temp_ye(rho_corot_atm, tempL, YeL); double alpL = calc_avg_v2c(alp, p); From 655ef2db337687d0f35f1305bf504c2b0a5c86e4 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Sat, 28 Mar 2026 18:40:27 -0500 Subject: [PATCH 06/11] RNSReader: set corotating atmosphere to have cylindrical symmetry --- RNSReader/src/equil_util.cxx | 4 ++-- RNSReader/src/include/consts.h | 1 - RNSReader/src/rnsid_rfr.cxx | 5 +++++ RNSReader/src/rnsreader_utils.hxx | 15 +++++++++++++++ RNSReader/src/set_vatmo.cxx | 28 +++++++++++++++------------- 5 files changed, 37 insertions(+), 16 deletions(-) diff --git a/RNSReader/src/equil_util.cxx b/RNSReader/src/equil_util.cxx index 1bb0ad4..808db45 100644 --- a/RNSReader/src/equil_util.cxx +++ b/RNSReader/src/equil_util.cxx @@ -479,7 +479,7 @@ double zbrent_diff(double (*func)(double, double, double, double, double, fb = fc; fc = fa; } - tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol; + tol1 = 2.0 * DBL_EPSILON * fabs(b) + 0.5 * tol; xm = 0.5 * (c - b); if (fabs(xm) <= tol1 || fb == 0.0) { return_value = b; @@ -566,7 +566,7 @@ double zbrent_rot(double (*func)(double, double, double, double, double, double, fb = fc; fc = fa; } - tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol; + tol1 = 2.0 * DBL_EPSILON * fabs(b) + 0.5 * tol; xm = 0.5 * (c - b); if (fabs(xm) <= tol1 || fb == 0.0) { return_value = b; diff --git a/RNSReader/src/include/consts.h b/RNSReader/src/include/consts.h index fa2e6c1..34999a6 100644 --- a/RNSReader/src/include/consts.h +++ b/RNSReader/src/include/consts.h @@ -91,4 +91,3 @@ #define MAX_NTAB 16001 /* MAX ELEMENTS IN EOS TAB */ #define MAX_ITERATION 2000 /* max number of iteration in iterate (euil.c) */ #define ITMAX 2000 /* MAX ITERATION IN BRENT (equil_util.c) */ -#define EPS DBL_EPSILON diff --git a/RNSReader/src/rnsid_rfr.cxx b/RNSReader/src/rnsid_rfr.cxx index 4e3bdde..5086444 100644 --- a/RNSReader/src/rnsid_rfr.cxx +++ b/RNSReader/src/rnsid_rfr.cxx @@ -431,6 +431,11 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { alp_cc(p.I) = exp_nu_ijk; + dtalp_cc(p.I) = 0.0; + dtbetax_cc(p.I) = 0.0; + dtbetay_cc(p.I) = 0.0; + dtbetaz_cc(p.I) = 0.0; + if (x_i == 0.0 && y_j == 0.0) { gxx_cc(p.I) = SQ(exp_alpha_ijk); diff --git a/RNSReader/src/rnsreader_utils.hxx b/RNSReader/src/rnsreader_utils.hxx index 05d730d..7ee7b54 100644 --- a/RNSReader/src/rnsreader_utils.hxx +++ b/RNSReader/src/rnsreader_utils.hxx @@ -74,6 +74,16 @@ public: // First find the central interpolation stencil index: int idx = bisection_idx_finder(th, numlines_in_file, th_arr); + if (idx == 0) + { + *f_of_th = v_th_arr[0]; + return; + } + if (idx == numlines_in_file - 1) + { + *f_of_th = v_th_arr[numlines_in_file - 1]; + return; + } #ifdef MAX #undef MAX @@ -122,6 +132,11 @@ public: int x2 = numlines_in_file - 1; CCTK_REAL y1 = rrbar - rbar_arr[x1]; CCTK_REAL y2 = rrbar - rbar_arr[x2]; + if (rrbar <= rbar_arr[0]) + return 0; + if (rrbar >= rbar_arr[numlines_in_file - 1]) + return numlines_in_file - 1; + if (y1 * y2 >= 0) { // Cannot print on GPU // fprintf(stderr,"INTERPOLATION BRACKETING ERROR %e | %e diff --git a/RNSReader/src/set_vatmo.cxx b/RNSReader/src/set_vatmo.cxx index 78a3d57..056e649 100644 --- a/RNSReader/src/set_vatmo.cxx +++ b/RNSReader/src/set_vatmo.cxx @@ -29,8 +29,8 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { evolution_eos); } + // Open the Omega(r_cyl) file. Everything still says vel(theta) because I am a lazy asshole Omega_th_reader* vel_th_reader = nullptr; - // Open the Omega(th) file FILE *vel_th_file = fopen(vel_th_filename, "r"); // Check if everything is OK with the file @@ -53,11 +53,11 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const vec, dim> gf_beta{betax, betay, betaz}; // Calculate grading of corot. atmo. to match numerical atmo. at r_corot_max - const CCTK_REAL rho_atm_r = (r_corot_max > r_atmo) - ? (rho_abs_min * pow((r_atmo / r_corot_max), n_rho_atmo)) - : rho_abs_min; - const CCTK_REAL n_corot = (log(rho_atm_r) - log(10 * rho_abs_min)) / (log(r_atmo) - log(r_corot_max)); - CCTK_VINFO("Calculated n_corot to be %e.", n_corot); + // const CCTK_REAL rho_atm_r = (r_corot_max > r_atmo) + // ? (rho_abs_min * pow((r_atmo / r_corot_max), n_rho_atmo)) + // : rho_abs_min; + // const CCTK_REAL n_corot = (log(rho_atm_r) - log(10 * rho_abs_min)) / (log(r_atmo) - log(r_corot_max)); + // CCTK_VINFO("Calculated n_corot to be %e.", n_corot); grid.loop_all<1, 1, 1>(grid.nghostzones, [=] CCTK_HOST(const Loop::PointDesc &p) @@ -69,6 +69,7 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL YeL = Ye(p.I); const CCTK_REAL epsL = eps(p.I); const CCTK_REAL pressL = press(p.I); + const CCTK_REAL cyl_distance = sqrt(p.x * p.x + p.y * p.y); const CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); // Grading rho @@ -78,10 +79,10 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol); if (rhoL <= rho_atmo_cut && radial_distance < r_corot_max) { - CCTK_REAL costh = std::abs(p.z) / radial_distance; // cos(theta) on RNS grid runs from 0 to 1 + // CCTK_REAL costh = std::abs(p.z) / radial_distance; // cos(theta) on RNS grid runs from 0 to 1 CCTK_REAL omg_surf; vel_th_reader->interpolate_1d_quantity_as_function_of_th( - MDIV - 1, costh, &omg_surf); // Interp omega to cos(theta) from file + MDIV - 1, cyl_distance, &omg_surf); // This used to be tabulated in cos(th), now I tabulate in r_cyl // Grading Omega CCTK_REAL omg_atm = (radial_distance > r_corot_min) @@ -89,12 +90,13 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { : omg_surf; // Pad density and recalculate primitives - CCTK_REAL rho_corot_atm = (radial_distance > r_atmo) - ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) - : 10 * rho_abs_min; + // CCTK_REAL rho_corot_atm = (radial_distance > r_atmo) + // ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) + // : 10 * rho_abs_min; + CCTK_REAL rho_corot_atm = 10.0 * rho_atm; - CCTK_REAL eps_corot_atm = eos_3p_tab3d->eps_from_valid_rho_temp_ye(rho_corot_atm, tempL, YeL); - CCTK_REAL press_corot_atm = eos_3p_tab3d->press_from_valid_rho_temp_ye(rho_corot_atm, tempL, YeL); + CCTK_REAL eps_corot_atm = eos_3p_tab3d->eps_from_rho_temp_ye(rho_corot_atm, tempL, YeL); + CCTK_REAL press_corot_atm = eos_3p_tab3d->press_from_rho_temp_ye(rho_corot_atm, tempL, YeL); double alpL = calc_avg_v2c(alp, p); From 7b54130a3d1b81367b08be44e5df3d3ef25f2bfa Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Mon, 30 Mar 2026 14:29:34 -0500 Subject: [PATCH 07/11] MeudonBNSX: remove redeclaration of Lorene namespace. Clarify EOS compatibility comments --- MeudonBNSX/src/Bin_NS.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) 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.); From c9baadd74845cb0277250cddcdabc239a0100a0a Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Tue, 31 Mar 2026 12:37:27 -0500 Subject: [PATCH 08/11] VolumeIntegrals_GRMHDX: WIP add tracking data for individual NSs --- VolumeIntegrals_GRMHDX/interface.ccl | 10 +++++++ .../src/perform_integration_sum_global.cxx | 28 +++++++++++++++++++ 2 files changed, 38 insertions(+) diff --git a/VolumeIntegrals_GRMHDX/interface.ccl b/VolumeIntegrals_GRMHDX/interface.ccl index c54058d..30e8a64 100644 --- a/VolumeIntegrals_GRMHDX/interface.ccl +++ b/VolumeIntegrals_GRMHDX/interface.ccl @@ -37,3 +37,13 @@ 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" diff --git a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx index 0a21f09..5d4d73a 100644 --- a/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx +++ b/VolumeIntegrals_GRMHDX/src/perform_integration_sum_global.cxx @@ -176,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]); + } + } } From 90a04d0c9c02fd36ac41ede6dd63e5cd36957ea0 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Wed, 1 Apr 2026 13:24:22 -0500 Subject: [PATCH 09/11] VolumeIntegrals_GRMHDX: store tracked NS CoMs in array --- VolumeIntegrals_GRMHDX/interface.ccl | 30 ++++++++++--------- VolumeIntegrals_GRMHDX/schedule.ccl | 10 +++---- ...nitialize_or_Decrement_IntegralCounter.cxx | 6 ++++ 3 files changed, 27 insertions(+), 19 deletions(-) diff --git a/VolumeIntegrals_GRMHDX/interface.ccl b/VolumeIntegrals_GRMHDX/interface.ccl index 30e8a64..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,17 +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" - -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" 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++) { From 7a3528306894bd01153611b9f8bea276cf28af61 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Sun, 5 Apr 2026 16:44:27 -0500 Subject: [PATCH 10/11] MHD_Diagnostics: clean up comments. General changes to make merging with main easier --- MHD_Diagnostics/src/master.cxx | 26 +- NRPyPlusTOVID/configuration.ccl | 3 - NRPyPlusTOVID/interface.ccl | 18 -- NRPyPlusTOVID/param.ccl | 66 ----- NRPyPlusTOVID/schedule.ccl | 44 --- NRPyPlusTOVID/src/ADMQuantities.h | 67 ----- NRPyPlusTOVID/src/HydroQuantities.h | 41 --- NRPyPlusTOVID/src/InitialData.cxx | 191 ------------- .../convert_TOV_spacetime_vars_to_ADM_vars.h | 40 --- .../src/interpolate_TOV_solution_to_point.h | 43 --- NRPyPlusTOVID/src/make.code.defn | 1 - NRPyPlusTOVID/src/tov_interp.h | 254 ------------------ RNSReader/param.ccl | 25 +- RNSReader/schedule.ccl | 7 +- RNSReader/src/equil_util.cxx | 4 +- RNSReader/src/include/consts.h | 1 + RNSReader/src/rnsid_rfr.cxx | 18 +- RNSReader/src/rnsreader_utils.hxx | 15 -- RNSReader/src/set_vatmo.cxx | 45 +--- 19 files changed, 30 insertions(+), 879 deletions(-) delete mode 100644 NRPyPlusTOVID/configuration.ccl delete mode 100644 NRPyPlusTOVID/interface.ccl delete mode 100644 NRPyPlusTOVID/param.ccl delete mode 100644 NRPyPlusTOVID/schedule.ccl delete mode 100644 NRPyPlusTOVID/src/ADMQuantities.h delete mode 100644 NRPyPlusTOVID/src/HydroQuantities.h delete mode 100644 NRPyPlusTOVID/src/InitialData.cxx delete mode 100644 NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h delete mode 100644 NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h delete mode 100644 NRPyPlusTOVID/src/make.code.defn delete mode 100644 NRPyPlusTOVID/src/tov_interp.h diff --git a/MHD_Diagnostics/src/master.cxx b/MHD_Diagnostics/src/master.cxx index f7cd5aa..befda0a 100644 --- a/MHD_Diagnostics/src/master.cxx +++ b/MHD_Diagnostics/src/master.cxx @@ -387,36 +387,12 @@ 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? + // 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); - /* - 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); - */ - // =========================================== double my_rho = rho(p.I); diff --git a/NRPyPlusTOVID/configuration.ccl b/NRPyPlusTOVID/configuration.ccl deleted file mode 100644 index b56a07c..0000000 --- a/NRPyPlusTOVID/configuration.ccl +++ /dev/null @@ -1,3 +0,0 @@ -# Configuration definition for thorn NRPyPlusTOVID - -REQUIRES Loop diff --git a/NRPyPlusTOVID/interface.ccl b/NRPyPlusTOVID/interface.ccl deleted file mode 100644 index 1403f7c..0000000 --- a/NRPyPlusTOVID/interface.ccl +++ /dev/null @@ -1,18 +0,0 @@ -implements: NRPyPlusTOVID -inherits: ADMBaseX HydroBaseX - -USES INCLUDE HEADER: loop.hxx -USES INCLUDE HEADER: loop_device.hxx - -PUBLIC: - -CCTK_REAL metric_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { gxx_cc gxy_cc gxz_cc gyy_cc gyz_cc gzz_cc } "cell-centered ADM 3-metric g_ij" - -CCTK_REAL curv_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { kxx_cc kxy_cc kxz_cc kyy_cc kyz_cc kzz_cc } "cell-centered ADM extrinstic curvature K_ij" - -CCTK_REAL lapse_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { alp_cc } "cell-centered ADM lapse function alpha" - -CCTK_REAL shift_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { betax_cc betay_cc betaz_cc} "cell-centered ADM shift function beta^i" - -CCTK_REAL dtlapse_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtalp_cc } "cell-centered Time derivative of ADM lapse function" -CCTK_REAL dtshift_cc TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtbetax_cc dtbetay_cc dtbetaz_cc} "cell-centered Time derivative of ADM shift function" diff --git a/NRPyPlusTOVID/param.ccl b/NRPyPlusTOVID/param.ccl deleted file mode 100644 index eb68e58..0000000 --- a/NRPyPlusTOVID/param.ccl +++ /dev/null @@ -1,66 +0,0 @@ -shares: ADMBaseX - -EXTENDS KEYWORD initial_data -{ - "NRPyPlusTOVID" :: "Initial data from NRPyPlusTOVID solution" -} -EXTENDS KEYWORD initial_lapse -{ - "NRPyPlusTOVID" :: "Initial lapse from NRPyPlusTOVID solution" -} -EXTENDS KEYWORD initial_shift -{ - "NRPyPlusTOVID" :: "Initial shift from NRPyPlusTOVID solution" -} -EXTENDS KEYWORD initial_dtlapse -{ - "NRPyPlusTOVID" :: "Initial dtlapse from NRPyPlusTOVID solution" -} -EXTENDS KEYWORD initial_dtshift -{ - "NRPyPlusTOVID" :: "Initial dtshift from NRPyPlusTOVID solution" -} - -shares: HydroBaseX -EXTENDS KEYWORD initial_hydro -{ - "NRPyPlusTOVID" :: "Initial GRHD data from NRPyPlusTOVID solution" -} - -restricted: -CCTK_STRING TOV_filename "Which interpolator should I use" -{ - ".+" :: "Any nonempty string" -} "outputTOVpolytrope.txt" - -restricted: -CCTK_REAL rho_atmosphere "Atmosphere baryonic density" -{ - 0:* :: "Physical values" --1 :: "forbidden value to make sure it is explicitly set in the parfile" -} -1 - -restricted: -CCTK_REAL K_atmosphere "Polytropic K to be used with the EOS corresponding to rho_atmosphere" -{ - 0:* :: "Physical values" --1 :: "forbidden value to make sure it is explicitly set in the parfile" -} -1 - -restricted: -CCTK_REAL Gamma_atmosphere "Polytropic Gamma to be used with the EOS corresponding to rho_atmosphere" -{ - 0:* :: "Physical values" --1 :: "forbidden value to make sure it is explicitly set in the parfile" -} -1 - -CCTK_REAL pert_amp "velocity perturbation amplitude" -{ - 0:* :: "Physical values" -} 0.0 - -restricted: -CCTK_REAL Pressure_depletion_factor "Pressure depletion factor = Pdf: P => (1-Pdf)*P" -{ - 0:* :: "Greater than or equal to zero, where zero is no depletion and default." -} 0.0 diff --git a/NRPyPlusTOVID/schedule.ccl b/NRPyPlusTOVID/schedule.ccl deleted file mode 100644 index 5921e24..0000000 --- a/NRPyPlusTOVID/schedule.ccl +++ /dev/null @@ -1,44 +0,0 @@ - -schedule NRPyPlusTOVID_ET_InitialData IN HydroBaseX_InitialData -{ - LANG: C - WRITES: metric_cc(everywhere) - WRITES: lapse_cc(everywhere) - WRITES: shift_cc(everywhere) - # WRITES: curv_cc(everywhere) - # WRITES: dtlapse_cc(everywhere) - # WRITES: dtshift_cc(everywhere) - WRITES: HYDROBASEX::eps(everywhere) - WRITES: HYDROBASEX::rho(everywhere) - WRITES: HYDROBASEX::vel(everywhere) - WRITES: HYDROBASEX::press(everywhere) -} "Set up general relativistic hydrodynamic (GRHD) fields for TOV initial data" - -schedule NRPyPlusTOVID_ET_VelPert IN HydroBaseX_InitialData AFTER NRPyPlusTOVID_ET_InitialData -{ - LANG: C - READS: HYDROBASEX::rho(everywhere) - WRITES: HYDROBASEX::vel(everywhere) -} "Set up general relativistic hydrodynamic (GRHD) fields for TOV initial data" - -SCHEDULE NRPyPlusTOVID_Interpolation_C2V IN HydroBaseX_InitialData AFTER NRPyPlusTOVID_ET_InitialData -{ - LANG: C - READS: metric_cc(everywhere) - READS: lapse_cc(everywhere) - READS: shift_cc(everywhere) - # READS: curv_cc(everywhere) - # READS: dtlapse_cc(everywhere) - # READS: dtshift_cc(everywhere) - - WRITES: ADMBaseX::metric(interior) - WRITES: ADMBaseX::lapse(interior) - WRITES: ADMBaseX::shift(interior) - WRITES: ADMBaseX::curv(interior) - WRITES: ADMBaseX::dtlapse(interior) - WRITES: ADMBaseX::dtshift(interior) - - SYNC: ADMBaseX::metric ADMBaseX::lapse ADMBaseX::shift ADMBaseX::curv - SYNC: ADMBaseX::dtlapse ADMBaseX::dtshift -} "Interpolate metric from cell center to vertex" - diff --git a/NRPyPlusTOVID/src/ADMQuantities.h b/NRPyPlusTOVID/src/ADMQuantities.h deleted file mode 100644 index a01b3ec..0000000 --- a/NRPyPlusTOVID/src/ADMQuantities.h +++ /dev/null @@ -1,67 +0,0 @@ - -namespace NRPyPlusTOVID{ -using namespace std; -using namespace Loop; - -static inline -void ADMQuantities(const PointDesc &p, - const CCTK_REAL IDalpha, - const CCTK_REAL IDgammaDD00,const CCTK_REAL IDgammaDD01, const CCTK_REAL IDgammaDD02, - const CCTK_REAL IDgammaDD11,const CCTK_REAL IDgammaDD12, const CCTK_REAL IDgammaDD22, - - const GF3D2 &alphaGF, - const GF3D2 &betaU0GF, - const GF3D2 &betaU1GF, - const GF3D2 &betaU2GF, - const GF3D2 &gammaDD00GF, - const GF3D2 &gammaDD01GF, - const GF3D2 &gammaDD02GF, - const GF3D2 &gammaDD11GF, - const GF3D2 &gammaDD12GF, - const GF3D2 &gammaDD22GF) { - const CCTK_REAL xx0 = p.x; - const CCTK_REAL xx1 = p.y; - const CCTK_REAL xx2 = p.z; - - /* - * NRPy+ Finite Difference Code Generation, Step 1 of 1: Evaluate SymPy expressions and write to main memory: - */ - const double FDPart3_0 = ((xx1)*(xx1)); - const double FDPart3_1 = ((xx0)*(xx0)); - const double FDPart3_2 = FDPart3_0 + FDPart3_1; - const double FDPart3_3 = IDgammaDD22/((FDPart3_2)*(FDPart3_2)); - const double FDPart3_4 = ((xx2)*(xx2)); - const double FDPart3_5 = FDPart3_2 + FDPart3_4; - const double FDPart3_6 = (1.0/(FDPart3_5)); - const double FDPart3_7 = FDPart3_6*IDgammaDD00; - const double FDPart3_8 = (1.0/(FDPart3_2)); - const double FDPart3_9 = (1.0/sqrt(FDPart3_5)); - const double FDPart3_11 = FDPart3_8*FDPart3_9*IDgammaDD02; - const double FDPart3_12 = xx0*xx1; - const double FDPart3_13 = 2*FDPart3_12; - const double FDPart3_15 = (1.0/((FDPart3_5)*(FDPart3_5))); - const double FDPart3_17 = -FDPart3_4*FDPart3_6 + 1; - const double FDPart3_18 = (1.0/sqrt(FDPart3_17)); - const double FDPart3_19 = FDPart3_18*IDgammaDD01; - const double FDPart3_20 = 2*FDPart3_19*xx2; - const double FDPart3_22 = IDgammaDD11/FDPart3_17; - const double FDPart3_23 = FDPart3_22*FDPart3_4/((FDPart3_5)*(FDPart3_5)*(FDPart3_5)); - const double FDPart3_25 = FDPart3_18*FDPart3_8*IDgammaDD12; - const double FDPart3_26 = pow(FDPart3_5, -3.0/2.0); - const double FDPart3_27 = FDPart3_26*xx2; - const double FDPart3_29 = FDPart3_13*FDPart3_25*FDPart3_27; - const double FDPart3_35 = -FDPart3_26*FDPart3_4 + FDPart3_9; - const double FDPart3_36 = FDPart3_35*xx1; - alphaGF(p.I) = IDalpha; - betaU0GF(p.I) = 0; - betaU1GF(p.I) = 0; - betaU2GF(p.I) = 0; - gammaDD00GF(p.I) = FDPart3_0*FDPart3_3 + FDPart3_1*FDPart3_15*FDPart3_20 + FDPart3_1*FDPart3_23 + FDPart3_1*FDPart3_7 - FDPart3_11*FDPart3_13 - FDPart3_29; - gammaDD01GF(p.I) = -FDPart3_0*FDPart3_11 - FDPart3_0*FDPart3_25*FDPart3_27 + FDPart3_1*FDPart3_18*FDPart3_27*FDPart3_8*IDgammaDD12 + FDPart3_1*FDPart3_8*FDPart3_9*IDgammaDD02 + FDPart3_12*FDPart3_23 - FDPart3_12*FDPart3_3 + FDPart3_12*FDPart3_7 + FDPart3_13*FDPart3_15*FDPart3_19*xx2; - gammaDD02GF(p.I) = -FDPart3_11*xx1*xx2 + FDPart3_15*FDPart3_19*FDPart3_4*xx0 - FDPart3_19*FDPart3_35*FDPart3_9*xx0 - FDPart3_22*FDPart3_27*FDPart3_35*xx0 + FDPart3_25*FDPart3_36 + FDPart3_7*xx0*xx2; - gammaDD11GF(p.I) = FDPart3_0*FDPart3_15*FDPart3_20 + FDPart3_0*FDPart3_23 + FDPart3_0*FDPart3_7 + FDPart3_1*FDPart3_3 + FDPart3_11*FDPart3_13 + FDPart3_29; - gammaDD12GF(p.I) = FDPart3_11*xx0*xx2 + FDPart3_15*FDPart3_19*FDPart3_4*xx1 - FDPart3_19*FDPart3_36*FDPart3_9 - FDPart3_22*FDPart3_27*FDPart3_36 - FDPart3_25*FDPart3_35*xx0 + FDPart3_7*xx1*xx2; - gammaDD22GF(p.I) = -FDPart3_20*FDPart3_35*FDPart3_9 + FDPart3_22*((FDPart3_35)*(FDPart3_35)) + FDPart3_4*FDPart3_6*IDgammaDD00; -} - -} //namespace diff --git a/NRPyPlusTOVID/src/HydroQuantities.h b/NRPyPlusTOVID/src/HydroQuantities.h deleted file mode 100644 index 5e8b1a1..0000000 --- a/NRPyPlusTOVID/src/HydroQuantities.h +++ /dev/null @@ -1,41 +0,0 @@ - - -namespace NRPyPlusTOVID{ -using namespace std; -using namespace Loop; - -static inline -void HydroQuantities(const PointDesc &p, - const CCTK_REAL IDPressure, const CCTK_REAL IDrho_baryonic, - const CCTK_REAL IDrho__total_energy_density, - const GF3D2 &PressureGF, - const GF3D2 &rho_baryonicGF, - const GF3D2 &epsilonGF, - const GF3D2 &Valencia3velocityU0GF, - const GF3D2 &Valencia3velocityU1GF, - const GF3D2 &Valencia3velocityU2GF) { - DECLARE_CCTK_PARAMETERS; - if(IDrho__total_energy_density <= 0 || IDrho_baryonic <= 0 || IDPressure <= 0) { - rho_baryonicGF(p.I) = rho_atmosphere; - PressureGF(p.I) = K_atmosphere*pow(rho_atmosphere,Gamma_atmosphere); - epsilonGF(p.I) = 0; - Valencia3velocityU0GF(p.I) = 0; - Valencia3velocityU1GF(p.I) = 0; - Valencia3velocityU2GF(p.I) = 0; - } else { - /* - * NRPy+ Finite Difference Code Generation, Step 1 of 1: Evaluate SymPy expressions and write to main memory: - */ - PressureGF(p.I) = IDPressure; - rho_baryonicGF(p.I) = IDrho_baryonic; - epsilonGF(p.I) = IDrho__total_energy_density/IDrho_baryonic - 1; - Valencia3velocityU0GF(p.I) = 0; - Valencia3velocityU1GF(p.I) = 0; - Valencia3velocityU2GF(p.I) = 0; - - // Apply pressure depletion. - PressureGF(p.I) *= (1.0 - Pressure_depletion_factor); - } -} - -} // namespace diff --git a/NRPyPlusTOVID/src/InitialData.cxx b/NRPyPlusTOVID/src/InitialData.cxx deleted file mode 100644 index 7b5e84a..0000000 --- a/NRPyPlusTOVID/src/InitialData.cxx +++ /dev/null @@ -1,191 +0,0 @@ - -#include -#include -#include - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - -#include -#include - -#include "aster_utils.hxx" - -#include "ADMQuantities.h" -#include "HydroQuantities.h" -#include "interpolate_TOV_solution_to_point.h" -#include "convert_TOV_spacetime_vars_to_ADM_vars.h" - -namespace NRPyPlusTOVID { -using namespace std; -using namespace AsterUtils; -using namespace Loop; - - -// Alias for "vel" vector gridfunction: -// #define velx (&vel[0*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) -// #define vely (&vel[1*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) -// #define velz (&vel[2*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]]) - -void read_TOV_input_data_from_file(ID_inputs *TOV_in) { - - DECLARE_CCTK_PARAMETERS; - - // Step 1: Set up TOV initial data - // Step 1.a: Read TOV initial data from data file - // Open the data file: - char filename[100]; - sprintf(filename,"%s",TOV_filename); // TOV_filename is a CCTK_PARAMETER - FILE *in1Dpolytrope = fopen(filename, "r"); - if (in1Dpolytrope == NULL) { - fprintf(stderr,"ERROR: could not open file %s\n",filename); - exit(1); - } - // Count the number of lines in the data file: - int numlines_in_file = count_num_lines_in_file(in1Dpolytrope); - // Allocate space for all data arrays: - CCTK_REAL *r_Schw_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *rho_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *rho_baryon_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *P_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *M_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *expnu_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *exp4phi_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - CCTK_REAL *rbar_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file); - - // Read from the data file, filling in arrays. - // read_datafile__set_arrays() may be found in TOV/tov_interp.h - if(read_datafile__set_arrays(in1Dpolytrope, r_Schw_arr,rho_arr,rho_baryon_arr,P_arr,M_arr,expnu_arr,exp4phi_arr,rbar_arr) == 1) { - fprintf(stderr,"ERROR WHEN READING FILE %s!\n",filename); - exit(1); - } - fclose(in1Dpolytrope); - REAL Rbar = -100; - int Rbar_idx = -100; - for(int i=1;i0 && rho_arr[i]==0) { Rbar = rbar_arr[i-1]; Rbar_idx = i-1; } - } - if(Rbar<0) { - fprintf(stderr,"Error: could not find rbar=Rbar from data file.\n"); - exit(1); - } - - TOV_in->Rbar = Rbar; - TOV_in->Rbar_idx = Rbar_idx; - - const int interp_stencil_size = 12; - TOV_in->interp_stencil_size = interp_stencil_size; - TOV_in->numlines_in_file = numlines_in_file; - - TOV_in->r_Schw_arr = r_Schw_arr; - TOV_in->rho_arr = rho_arr; - TOV_in->rho_baryon_arr = rho_baryon_arr; - TOV_in->P_arr = P_arr; - TOV_in->M_arr = M_arr; - TOV_in->expnu_arr = expnu_arr; - TOV_in->exp4phi_arr = exp4phi_arr; - TOV_in->rbar_arr = rbar_arr; - /* END TOV INPUT ROUTINE */ -} - -extern "C" void NRPyPlusTOVID_ET_InitialData(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_ET_InitialData; - DECLARE_CCTK_PARAMETERS; - - ID_inputs TOV_in; - read_TOV_input_data_from_file(&TOV_in); - - const array indextype = {1, 1, 1}; - const GF3D2layout CCC_layout(cctkGH, indextype); - - grid.loop_all<1, 1, 1>(grid.nghostzones, - [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - CCTK_INT idx = CCC_layout.linear(p.I); - CCTK_REAL rr = sqrt(p.x*p.x + p.y*p.y + p.z*p.z); - CCTK_REAL th = acos(p.z/rr); - - CCTK_REAL IDexp_4phi,IDnu,IDPressure,IDrho_baryonic,IDrho__total_energy_density; - interpolate_TOV_solution_to_point(rr, TOV_in, &IDexp_4phi,&IDnu, - &IDPressure,&IDrho_baryonic,&IDrho__total_energy_density); - - CCTK_REAL IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22; - convert_TOV_spacetime_vars_to_ADM_vars(rr, th, IDexp_4phi,IDnu, - &IDalpha,&IDgammaDD00,&IDgammaDD01,&IDgammaDD02,&IDgammaDD11,&IDgammaDD12,&IDgammaDD22); - - HydroQuantities(p, - IDPressure,IDrho_baryonic,IDrho__total_energy_density, - press,rho,eps,velx,vely,velz); - - ADMQuantities(p, - IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22, - alp_cc,betax_cc,betay_cc,betaz_cc, - gxx_cc,gxy_cc,gxz_cc,gyy_cc,gyz_cc,gzz_cc); - }); - - free(TOV_in.r_Schw_arr); - free(TOV_in.rho_arr); - free(TOV_in.rho_baryon_arr); - free(TOV_in.P_arr); - free(TOV_in.M_arr); - free(TOV_in.expnu_arr); - free(TOV_in.exp4phi_arr); - free(TOV_in.rbar_arr); -} - -extern "C" void NRPyPlusTOVID_Interpolation_C2V(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_Interpolation_C2V; - DECLARE_CCTK_PARAMETERS; - - CCTK_INFO("Starting interpolation for ADM variables."); - grid.loop_int<0, 0, 0>(grid.nghostzones, - [=] CCTK_HOST(const Loop::PointDesc &p) - CCTK_ATTRIBUTE_ALWAYS_INLINE { - gxx(p.I) = calc_avg_c2v<4>(gxx_cc, p); - gxy(p.I) = calc_avg_c2v<4>(gxy_cc, p); - gxz(p.I) = calc_avg_c2v<4>(gxz_cc, p); - gyy(p.I) = calc_avg_c2v<4>(gyy_cc, p); - gyz(p.I) = calc_avg_c2v<4>(gyz_cc, p); - gzz(p.I) = calc_avg_c2v<4>(gzz_cc, p); - - kxx(p.I) = 0.0; - kxy(p.I) = 0.0; - kxz(p.I) = 0.0; - kyy(p.I) = 0.0; - kyz(p.I) = 0.0; - kzz(p.I) = 0.0; - - alp(p.I) = calc_avg_c2v<4>(alp_cc, p); - betax(p.I) = calc_avg_c2v<4>(betax_cc, p); - betay(p.I) = calc_avg_c2v<4>(betay_cc, p); - betaz(p.I) = calc_avg_c2v<4>(betaz_cc, p); - - dtalp(p.I) = 0.0; - dtbetax(p.I) = 0.0; - dtbetay(p.I) = 0.0; - dtbetaz(p.I) = 0.0; - }); - - CCTK_INFO("Done interpolation for ADM variables."); -} - -extern "C" void NRPyPlusTOVID_ET_VelPert(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_NRPyPlusTOVID_ET_VelPert; - DECLARE_CCTK_PARAMETERS; - - grid.loop_all<1, 1, 1>(grid.nghostzones, - [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - if (rho(p.I) > rho_atmosphere) { - velx(p.I) = pert_amp * p.x; - vely(p.I) = pert_amp * p.y; - velz(p.I) = pert_amp * p.z; - } - else { - velx(p.I) = 0.0; - vely(p.I) = 0.0; - velz(p.I) = 0.0; - } - }); -} - -} // namespace diff --git a/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h b/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h deleted file mode 100644 index ca1d5e4..0000000 --- a/NRPyPlusTOVID/src/convert_TOV_spacetime_vars_to_ADM_vars.h +++ /dev/null @@ -1,40 +0,0 @@ - -/* This function converts TOV quantities into - * ADM quantities in Spherical coordinates. - */ -void convert_TOV_spacetime_vars_to_ADM_vars( const CCTK_REAL rr, const CCTK_REAL th, - const CCTK_REAL IDexp_4phi, const CCTK_REAL IDexpnu, - CCTK_REAL *IDalpha, - CCTK_REAL *IDgammaDD00, CCTK_REAL *IDgammaDD01, CCTK_REAL *IDgammaDD02, - CCTK_REAL *IDgammaDD11, CCTK_REAL *IDgammaDD12, CCTK_REAL *IDgammaDD22) { - - /*************************************************************** - * Convert TOV quantities to ADM quantities in Spherical basis * - *************************************************************** - * - * First we convert the lapse function: - * .------------------. - * | alpha = e^(nu/2) | - * .------------------. - */ - *IDalpha = sqrt(IDexpnu); - - /* Next we convert the metric function: - * .----------------------------------------. - * | gamma_{00} = e^{4phi} | - * .----------------------------------------. - * | gamma_{11} = e^{4phi} r^2 | - * .----------------------------------------. - * | gamma_{22} = e^{4phi} r^2 sin^2(theta) | - * .----------------------------------------. - * | All other components are zero. | - * .----------------------------------------. - */ - *IDgammaDD00 = IDexp_4phi; - *IDgammaDD11 = IDexp_4phi * rr * rr; - *IDgammaDD22 = IDexp_4phi * rr * rr * sin(th) * sin(th); - *IDgammaDD01 = 0.0; - *IDgammaDD02 = 0.0; - *IDgammaDD12 = 0.0; - -} diff --git a/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h b/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h deleted file mode 100644 index ea98ad6..0000000 --- a/NRPyPlusTOVID/src/interpolate_TOV_solution_to_point.h +++ /dev/null @@ -1,43 +0,0 @@ - -/* Load the TOV_interpolate_1D() function */ -#include "tov_interp.h" - -// Declare initial data input struct: -// stores data from initial data solver, -// so they can be put on the numerical grid. -typedef struct __ID_inputs { - CCTK_REAL Rbar; - int Rbar_idx; - int interp_stencil_size; - int numlines_in_file; - CCTK_REAL *r_Schw_arr,*rho_arr,*rho_baryon_arr,*P_arr,*M_arr,*expnu_arr,*exp4phi_arr,*rbar_arr; -} ID_inputs; - -/* This function returns the TOV quantities at point rr - * by interpolating the data in the TOV initial data file. - */ -void interpolate_TOV_solution_to_point(const CCTK_REAL rr, ID_inputs other_inputs, - CCTK_REAL *exp_4phi, CCTK_REAL *expnu, - CCTK_REAL *Pressure, CCTK_REAL *rho_baryon, CCTK_REAL *rho__total_energy_density) { - - /* The mass valus is not used, but we have to - * store it in this dummy variable because the - * initial data file contains it. - */ - CCTK_REAL M; - - /* Perform the interpolation, returning: - * - rho__total_energy_density - * - rho_baryon - * - Pressure - * - Mass (dummy variable, unused) - * - exp(nu) - * - exp(4phi) - */ - TOV_interpolate_1D(rr,other_inputs.Rbar,other_inputs.Rbar_idx,other_inputs.interp_stencil_size, - other_inputs.numlines_in_file, - other_inputs.r_Schw_arr,other_inputs.rho_arr,other_inputs.rho_baryon_arr,other_inputs.P_arr,other_inputs.M_arr, - other_inputs.expnu_arr,other_inputs.exp4phi_arr,other_inputs.rbar_arr, - rho__total_energy_density,rho_baryon,Pressure,&M,expnu,exp_4phi); - -} diff --git a/NRPyPlusTOVID/src/make.code.defn b/NRPyPlusTOVID/src/make.code.defn deleted file mode 100644 index a8dea9a..0000000 --- a/NRPyPlusTOVID/src/make.code.defn +++ /dev/null @@ -1 +0,0 @@ -SRCS = InitialData.cxx diff --git a/NRPyPlusTOVID/src/tov_interp.h b/NRPyPlusTOVID/src/tov_interp.h deleted file mode 100644 index e22c0d3..0000000 --- a/NRPyPlusTOVID/src/tov_interp.h +++ /dev/null @@ -1,254 +0,0 @@ -// This C header file reads a TOV solution from data file and performs -// 1D interpolation of the solution to a desired radius. - -// Author: Zachariah B. Etienne -// zachetie **at** gmail **dot* com - -#include "stdio.h" -#include "stdlib.h" -#include "math.h" -#include "string.h" - -#define REAL double - -//#define STANDALONE_UNIT_TEST - -int count_num_lines_in_file(FILE *in1Dpolytrope) { - int numlines_in_file = 0; - char * line = NULL; - - size_t len = 0; - ssize_t read; - while ((read = getline(&line, &len, in1Dpolytrope)) != -1) { - numlines_in_file++; - } - rewind(in1Dpolytrope); - - free(line); - return numlines_in_file; -} - -int read_datafile__set_arrays(FILE *in1Dpolytrope, REAL *restrict r_Schw_arr,REAL *restrict rho_arr,REAL *restrict rho_baryon_arr,REAL *restrict P_arr, - REAL *restrict M_arr,REAL *restrict expnu_arr,REAL *restrict exp4phi_arr,REAL *restrict rbar_arr) { - char * line = NULL; - - size_t len = 0; - ssize_t read; - - int which_line = 0; - while ((read = getline(&line, &len, in1Dpolytrope)) != -1) { - // Define the line delimiters (i.e., the stuff that goes between the data on a given - // line of data. Here, we define both spaces " " and tabs "\t" as data delimiters. - const char delimiters[] = " \t"; - - //Now we define "token", a pointer to the first column of data - char *token; - - //Each successive time we call strtok(NULL,blah), we read in a new column of data from - // the originally defined character array, as pointed to by token. - - token=strtok(line, delimiters); if(token==NULL) { printf("BADDDD\n"); return 1; } - r_Schw_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - rho_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - rho_baryon_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - P_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - M_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - expnu_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - exp4phi_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters ); - rbar_arr[which_line] = strtod(token, NULL); - - which_line++; - } - free(line); - return 0; -} - -// Find interpolation index using Bisection root-finding algorithm: -static inline int bisection_idx_finder(const REAL rrbar, const int numlines_in_file, const REAL *restrict rbar_arr) { - int x1 = 0; - int x2 = numlines_in_file-1; - REAL y1 = rrbar-rbar_arr[x1]; - REAL y2 = rrbar-rbar_arr[x2]; - if(y1*y2 >= 0) { - fprintf(stderr,"INTERPOLATION BRACKETING ERROR %e | %e %e\n",rrbar,y1,y2); - exit(1); - } - for(int i=0;i (B)) ? (A) : (B) ) - - int idxmin = MAX(0,idx-interp_stencil_size/2-1); - -#ifdef MIN -#undef MIN -#endif -#define MIN(A, B) ( ((A) < (B)) ? (A) : (B) ) - - // -= Do not allow the interpolation stencil to cross the star's surface =- - // max index is when idxmin + (interp_stencil_size-1) = Rbar_idx - // -> idxmin at most can be Rbar_idx - interp_stencil_size + 1 - if(rrbar < Rbar) { - idxmin = MIN(idxmin,Rbar_idx - interp_stencil_size + 1); - } else { - idxmin = MAX(idxmin,Rbar_idx+1); - idxmin = MIN(idxmin,numlines_in_file - interp_stencil_size + 1); - } - // Now perform the Lagrange polynomial interpolation: - - // First set the interpolation coefficients: - REAL rbar_sample[interp_stencil_size]; - for(int i=idxmin;i Rbar) { - *rho = 0; - *rho_baryon = 0; - *P = 0; - *M = M_arr[Rbar_idx+1]; - *expnu = 1. - 2.*(*M) / r_Schw; - *exp4phi = pow(r_Schw / rrbar,2.0); - } -} - -// To compile, copy this file to tov_interp.c, and then run: -// gcc -Ofast tov_interp.c -o tov_interp -DSTANDALONE_UNIT_TEST -#ifdef STANDALONE_UNIT_TEST -int main() { - - // Open the data file: - char filename[100]; - sprintf(filename,"../outputTOVpolytrope.txt"); - FILE *in1Dpolytrope = fopen(filename, "r"); - if (in1Dpolytrope == NULL) { - fprintf(stderr,"ERROR: could not open file %s\n",filename); - exit(1); - } - - // Count the number of lines in the data file: - int numlines_in_file = count_num_lines_in_file(in1Dpolytrope); - - // Allocate space for all data arrays: - REAL *r_Schw_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *rho_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *rho_baryon_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *P_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *M_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *expnu_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *exp4phi_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - REAL *rbar_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file); - - // Read from the data file, filling in arrays - if(read_datafile__set_arrays(in1Dpolytrope, r_Schw_arr,rho_arr,rho_baryon_arr,P_arr,M_arr,expnu_arr,exp4phi_arr,rbar_arr) == 1) { - fprintf(stderr,"ERROR WHEN READING FILE %s!\n",filename); - exit(1); - } - fclose(in1Dpolytrope); - - REAL Rbar = -100; - int Rbar_idx = -100; - for(int i=1;i0 && rho_arr[i]==0) { Rbar = rbar_arr[i-1]; Rbar_idx = i-1; } - } - if(Rbar<0) { - fprintf(stderr,"Error: could not find r=R from data file.\n"); - exit(1); - } - - // Next, interpolate! - // Create trial radius array: - int num_r_pts = 100000; - //REAL *r_out_arr = (REAL *)malloc(sizeof(REAL)*num_r_pts); - struct drand48_data randBuffer; - srand48_r(1313, &randBuffer); -#pragma omp parallel for - for(int i=0;i 1e-10) { - // printf("Post RNSRead: rho = %e; eps = %e; press = %e;\n", rho(p.I), eps(p.I), press(p.I)); - // } - }); /* END GRID LOOP */ if (strcmp(zero_shift, "yes") == 0) { @@ -622,10 +610,6 @@ extern "C" void Hydro_rnsid_init(CCTK_ARGUMENTS) { free(s_gp); free(mu); - free(eos_type); - free(eos_file); - free(rotation_type); - /* Add code for filling time levels here if needed in the future. */ return; diff --git a/RNSReader/src/rnsreader_utils.hxx b/RNSReader/src/rnsreader_utils.hxx index 7ee7b54..05d730d 100644 --- a/RNSReader/src/rnsreader_utils.hxx +++ b/RNSReader/src/rnsreader_utils.hxx @@ -74,16 +74,6 @@ public: // First find the central interpolation stencil index: int idx = bisection_idx_finder(th, numlines_in_file, th_arr); - if (idx == 0) - { - *f_of_th = v_th_arr[0]; - return; - } - if (idx == numlines_in_file - 1) - { - *f_of_th = v_th_arr[numlines_in_file - 1]; - return; - } #ifdef MAX #undef MAX @@ -132,11 +122,6 @@ public: int x2 = numlines_in_file - 1; CCTK_REAL y1 = rrbar - rbar_arr[x1]; CCTK_REAL y2 = rrbar - rbar_arr[x2]; - if (rrbar <= rbar_arr[0]) - return 0; - if (rrbar >= rbar_arr[numlines_in_file - 1]) - return numlines_in_file - 1; - if (y1 * y2 >= 0) { // Cannot print on GPU // fprintf(stderr,"INTERPOLATION BRACKETING ERROR %e | %e diff --git a/RNSReader/src/set_vatmo.cxx b/RNSReader/src/set_vatmo.cxx index 056e649..b48cafb 100644 --- a/RNSReader/src/set_vatmo.cxx +++ b/RNSReader/src/set_vatmo.cxx @@ -5,7 +5,6 @@ #include "rnsreader_utils.hxx" #include "consts.h" #include "aster_utils.hxx" -#include "setup_eos.hxx" #include #include @@ -16,21 +15,13 @@ using namespace Loop; using namespace amrex; using namespace std; using namespace AsterUtils; -using namespace EOSX; extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_RNSReader_Init_VelAtmo; DECLARE_CCTK_PARAMETERS; - auto eos_3p_tab3d = global_eos_3p_tab3d; - if (not CCTK_EQUALS(evolution_eos, "Tabulated3d")) { - CCTK_VERROR("Invalid evolution EOS type '%s'. Please, set " - "EOSX::evolution_eos = \"Tabulated3d\" in your parameter file.", - evolution_eos); - } - - // Open the Omega(r_cyl) file. Everything still says vel(theta) because I am a lazy asshole Omega_th_reader* vel_th_reader = nullptr; + // Open the Omega(th) file FILE *vel_th_file = fopen(vel_th_filename, "r"); // Check if everything is OK with the file @@ -53,11 +44,11 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const vec, dim> gf_beta{betax, betay, betaz}; // Calculate grading of corot. atmo. to match numerical atmo. at r_corot_max - // const CCTK_REAL rho_atm_r = (r_corot_max > r_atmo) - // ? (rho_abs_min * pow((r_atmo / r_corot_max), n_rho_atmo)) - // : rho_abs_min; - // const CCTK_REAL n_corot = (log(rho_atm_r) - log(10 * rho_abs_min)) / (log(r_atmo) - log(r_corot_max)); - // CCTK_VINFO("Calculated n_corot to be %e.", n_corot); + const CCTK_REAL rho_atm_r = (r_corot_max > r_atmo) + ? (rho_abs_min * pow((r_atmo / r_corot_max), n_rho_atmo)) + : rho_abs_min; + const CCTK_REAL n_corot = (log(rho_atm_r) - log(10 * rho_abs_min)) / (log(r_atmo) - log(r_corot_max)); + CCTK_VINFO("Calculated n_corot to be %e.", n_corot); grid.loop_all<1, 1, 1>(grid.nghostzones, [=] CCTK_HOST(const Loop::PointDesc &p) @@ -65,11 +56,6 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL velxL = velx(p.I); const CCTK_REAL velyL = vely(p.I); const CCTK_REAL rhoL = rho(p.I); - const CCTK_REAL tempL = temperature(p.I); - const CCTK_REAL YeL = Ye(p.I); - const CCTK_REAL epsL = eps(p.I); - const CCTK_REAL pressL = press(p.I); - const CCTK_REAL cyl_distance = sqrt(p.x * p.x + p.y * p.y); const CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); // Grading rho @@ -79,24 +65,19 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol); if (rhoL <= rho_atmo_cut && radial_distance < r_corot_max) { - // CCTK_REAL costh = std::abs(p.z) / radial_distance; // cos(theta) on RNS grid runs from 0 to 1 + CCTK_REAL costh = std::abs(p.z) / radial_distance; // cos(theta) on RNS grid runs from 0 to 1 CCTK_REAL omg_surf; vel_th_reader->interpolate_1d_quantity_as_function_of_th( - MDIV - 1, cyl_distance, &omg_surf); // This used to be tabulated in cos(th), now I tabulate in r_cyl + MDIV - 1, costh, &omg_surf); // Interp omega to cos(theta) from file // Grading Omega CCTK_REAL omg_atm = (radial_distance > r_corot_min) ? (omg_surf * pow((r_corot_min / radial_distance), n_omg_atmo)) : omg_surf; - // Pad density and recalculate primitives - // CCTK_REAL rho_corot_atm = (radial_distance > r_atmo) - // ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) - // : 10 * rho_abs_min; - CCTK_REAL rho_corot_atm = 10.0 * rho_atm; - - CCTK_REAL eps_corot_atm = eos_3p_tab3d->eps_from_rho_temp_ye(rho_corot_atm, tempL, YeL); - CCTK_REAL press_corot_atm = eos_3p_tab3d->press_from_rho_temp_ye(rho_corot_atm, tempL, YeL); + CCTK_REAL rho_corot_atm = (radial_distance > r_atmo) + ? (10 * rho_abs_min * pow((r_atmo / radial_distance), n_corot)) + : 10 * rho_abs_min; double alpL = calc_avg_v2c(alp, p); @@ -106,15 +87,11 @@ extern "C" void RNSReader_Init_VelAtmo(CCTK_ARGUMENTS) { velx(p.I) = (-p.y * omg_atm + betas_avg(0)) / alpL; vely(p.I) = (p.x * omg_atm + betas_avg(1)) / alpL; rho(p.I) = rho_corot_atm; - eps(p.I) = eps_corot_atm; - press(p.I) = press_corot_atm; } else { velx(p.I) = velxL; vely(p.I) = velyL; rho(p.I) = rhoL; - eps(p.I) = epsL; - press(p.I) = pressL; } }); } From 19824347409f334f2a0e3bdc9e3d61dfb36744e8 Mon Sep 17 00:00:00 2001 From: Allen Wen Date: Wed, 8 Apr 2026 00:03:16 -0500 Subject: [PATCH 11/11] VolumeIntegrals_GRMHDX: add integrals for computing +/- asymmetry in magnetic fields --- VolumeIntegrals_GRMHDX/param.ccl | 4 + .../src/evaluate_integrands_local.cxx | 96 +++++++- VolumeIntegrals_GRMHDX/src/integrands.cxx | 221 ++++++++++++++++++ .../src/number_of_reductions.cxx | 10 + 4 files changed, 323 insertions(+), 8 deletions(-) diff --git a/VolumeIntegrals_GRMHDX/param.ccl b/VolumeIntegrals_GRMHDX/param.ccl index 884010f..e331c9f 100644 --- a/VolumeIntegrals_GRMHDX/param.ccl +++ b/VolumeIntegrals_GRMHDX/param.ccl @@ -86,12 +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/src/evaluate_integrands_local.cxx b/VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx index 238c0a5..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>( @@ -253,6 +297,50 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) { VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y); }); + } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], + "volume_average_norm_B_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; + } + + 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>( @@ -341,14 +429,6 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) { magnetic_co(VolIntegrand1, p, b2small, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, gamma_lim); }); - } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral], - "magnetic_energy_comov")) { - 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); - }); } else { /* Print a warning if no integrand is computed because * Integration_quantity_keyword unrecognized. */ diff --git a/VolumeIntegrals_GRMHDX/src/integrands.cxx b/VolumeIntegrals_GRMHDX/src/integrands.cxx index 91c7d71..0b3feec 100644 --- a/VolumeIntegrals_GRMHDX/src/integrands.cxx +++ b/VolumeIntegrals_GRMHDX/src/integrands.cxx @@ -307,6 +307,94 @@ 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, @@ -992,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 2374b2e..00b3a6b 100644 --- a/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx +++ b/VolumeIntegrals_GRMHDX/src/number_of_reductions.cxx @@ -33,6 +33,12 @@ 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; @@ -70,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;