diff --git a/AsterSeeds/param.ccl b/AsterSeeds/param.ccl index 5dfdf4df..32be9f9e 100644 --- a/AsterSeeds/param.ccl +++ b/AsterSeeds/param.ccl @@ -42,6 +42,28 @@ KEYWORD test_case "Name of the testcase" STEERABLE=never private: +# parameters for perturbing neutron-star initial data + +KEYWORD ns_perturbation_mode "Neutron-star perturbation applied after hydro initial data are read" STEERABLE=never +{ + "none" :: "Do not perturb the neutron-star initial data" + "m0" :: "Deplete pressure by ns_perturbation_strength inside the star" + "m1" :: "rho -> rho * (1 + B * (x + y) / Req)" + "m2" :: "rho -> rho * (1 + B * (x^2 - y^2) / Req^2)" +} "none" + +CCTK_REAL ns_perturbation_strength "Fractional perturbation strength B" +{ + *:* :: "" +} 0.0 + +CCTK_REAL ns_perturbation_Req "Equatorial radius Req used for m=1 and m=2 perturbations" +{ + (0.0:* :: "Positive" +} 1.0 + +# --------------------------------------- + # parameters for rotating intial data REAL rotate_angle_x "Rotation angle about x-axis (in units of pi)" STEERABLE=never { diff --git a/AsterSeeds/schedule.ccl b/AsterSeeds/schedule.ccl index 9744427e..324aa691 100644 --- a/AsterSeeds/schedule.ccl +++ b/AsterSeeds/schedule.ccl @@ -252,7 +252,7 @@ if (CCTK_Equals(test_type, "TabEOSTest")) { #Initial conditions for temperature, entropy and Ye -if (set_Ye_postinitial) { +if (set_Ye_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) { SCHEDULE SetYe IN HydroBaseX_PostInitial { LANG: C @@ -261,22 +261,89 @@ if (set_Ye_postinitial) { } "Set initial Ye" } -if (set_temperature_postinitial) { - SCHEDULE SetTemp IN HydroBaseX_PostInitial AFTER SetYe - { - LANG: C - READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) - WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere) - - } "Set initial temperature" +if (set_temperature_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) { + if (set_Ye_postinitial) { + SCHEDULE SetTemp IN HydroBaseX_PostInitial AFTER SetYe + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere) + + } "Set initial temperature" + } else { + SCHEDULE SetTemp IN HydroBaseX_PostInitial + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere) + + } "Set initial temperature" + } } -if (set_entropy_postinitial) { - SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetTemp - { - LANG: C - READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) - WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere) - - } "Set initial entropy" +if (set_entropy_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) { + if (set_temperature_postinitial) { + SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetTemp + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere) + + } "Set initial entropy" + } else if (set_Ye_postinitial) { + SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetYe + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere) + + } "Set initial entropy" + } else { + SCHEDULE SetEntropy IN HydroBaseX_PostInitial + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere) + + } "Set initial entropy" + } +} + +if (!CCTK_Equals(ns_perturbation_mode, "none")) { + if (CCTK_Equals(evolution_eos, "Tabulated3d")) { + SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER ID_TabEOS_HydroQuantities_recompute_HydroBase_variables + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere) + } "Apply neutron-star pressure or density perturbation after tabulated-EOS hydro quantities" + } else if (set_entropy_postinitial) { + SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetEntropy + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere) + } "Apply neutron-star pressure or density perturbation after hydro initial data" + } else if (set_temperature_postinitial) { + SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetTemp + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere) + } "Apply neutron-star pressure or density perturbation after hydro initial data" + } else if (set_Ye_postinitial) { + SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetYe + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere) + } "Apply neutron-star pressure or density perturbation after hydro initial data" + } else { + SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial + { + LANG: C + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere) + } "Apply neutron-star pressure or density perturbation after hydro initial data" + } } diff --git a/AsterSeeds/src/SetTempEntropyYe.cxx b/AsterSeeds/src/SetTempEntropyYe.cxx index 2c97d0fa..57b637c8 100644 --- a/AsterSeeds/src/SetTempEntropyYe.cxx +++ b/AsterSeeds/src/SetTempEntropyYe.cxx @@ -1,44 +1,48 @@ -#include #include #include #include #include -#include -#include -#include +#include #include +namespace AsterSeeds { +using namespace Loop; using namespace EOSX; + enum class eos_3param { IdealGas, Hybrid, Tabulated }; template -void SetTemp_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) { +void AsterSeeds_SetTemp_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) { DECLARE_CCTK_ARGUMENTSX_SetTemp; DECLARE_CCTK_PARAMETERS; grid.loop_all_device<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL epsL = eps(p.I); temperature(p.I) = - eos_3p->temp_from_rho_eps_ye(rho(p.I), eps(p.I), Ye(p.I)); + eos_3p->temp_from_rho_eps_ye(rho(p.I), epsL, Ye(p.I)); + eps(p.I) = epsL; }); } template -void SetEntropy_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) { +void AsterSeeds_SetEntropy_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) { DECLARE_CCTK_ARGUMENTSX_SetEntropy; DECLARE_CCTK_PARAMETERS; grid.loop_all_device<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL epsL = eps(p.I); entropy(p.I) = - eos_3p->kappa_from_rho_eps_ye(rho(p.I), eps(p.I), Ye(p.I)); + eos_3p->kappa_from_rho_eps_ye(rho(p.I), epsL, Ye(p.I)); + eps(p.I) = epsL; }); } @@ -74,16 +78,16 @@ extern "C" void SetTemp(CCTK_ARGUMENTS) { switch (eos_3p_type) { case eos_3param::IdealGas: { auto eos_3p_ig = global_eos_3p_ig; - SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_ig); + AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_ig); break; } case eos_3param::Hybrid: { - if (global_eos_3p_hyb_poly) { - auto eos_3p_hyb = global_eos_3p_hyb_poly; - SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); - } else if (global_eos_3p_hyb_pwpoly) { + if (global_eos_3p_hyb_pwpoly) { auto eos_3p_hyb = global_eos_3p_hyb_pwpoly; - SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + } else if (global_eos_3p_hyb_poly) { + auto eos_3p_hyb = global_eos_3p_hyb_poly; + AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); } else { CCTK_ERROR( "Hybrid EOS selected but no hybrid EOS object was initialized"); @@ -92,7 +96,7 @@ extern "C" void SetTemp(CCTK_ARGUMENTS) { } case eos_3param::Tabulated: { auto eos_3p_tab3d = global_eos_3p_tab3d; - SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d); + AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d); break; } default: @@ -121,16 +125,16 @@ extern "C" void SetEntropy(CCTK_ARGUMENTS) { switch (eos_3p_type) { case eos_3param::IdealGas: { auto eos_3p_ig = global_eos_3p_ig; - SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_ig); + AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_ig); break; } case eos_3param::Hybrid: { - if (global_eos_3p_hyb_poly) { - auto eos_3p_hyb = global_eos_3p_hyb_poly; - SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); - } else if (global_eos_3p_hyb_pwpoly) { + if (global_eos_3p_hyb_pwpoly) { auto eos_3p_hyb = global_eos_3p_hyb_pwpoly; - SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + } else if (global_eos_3p_hyb_poly) { + auto eos_3p_hyb = global_eos_3p_hyb_poly; + AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); } else { CCTK_ERROR( "Hybrid EOS selected but no hybrid EOS object was initialized"); @@ -139,10 +143,12 @@ extern "C" void SetEntropy(CCTK_ARGUMENTS) { } case eos_3param::Tabulated: { auto eos_3p_tab3d = global_eos_3p_tab3d; - SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d); + AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d); break; } default: assert(0); } } + +} // namespace AsterSeeds diff --git a/AsterSeeds/src/magTOV.cxx b/AsterSeeds/src/magTOV.cxx index 16dfa014..75ac3b2f 100644 --- a/AsterSeeds/src/magTOV.cxx +++ b/AsterSeeds/src/magTOV.cxx @@ -52,6 +52,34 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_TOV(CCTK_ARGUMENTS) { Avec_z_cent(p.I) = 0.0; }); + } else if (CCTK_EQUALS(Afield_config, "external dipole UIUC")) { + + /* computing cell centered vector potential components */ + grid.loop_all<1, 1, 1>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + const CCTK_REAL pi = 2 * acos(0.0); + + CCTK_REAL x_local = p.x - dipole_x[0]; + CCTK_REAL y_local = p.y - dipole_y[0]; + CCTK_REAL z_local = p.z - dipole_z[0]; + CCTK_REAL cylrad2 = x_local * x_local + y_local * y_local; + CCTK_REAL sphrad2 = + x_local * x_local + y_local * y_local + z_local * z_local; + CCTK_REAL r02 = r0 * r0; + + // See e.g. Ruiz+ 2020, Eq. 1. Here, B0 is I0 from Eq. 1, and we + // have preemptively canceled out the factor of cylrad2. + CCTK_REAL Aphi_local = + pi * B0 * r02 / pow(r02 + sphrad2, 1.5) * + (1.0 + (15.0 * r02 * (r02 + cylrad2) / + (8.0 * pow(r02 + sphrad2, 2.0)))); + + Avec_x_cent(p.I) = -y_local * Aphi_local; + Avec_y_cent(p.I) = x_local * Aphi_local; + Avec_z_cent(p.I) = 0.0; + }); + } else { CCTK_ERROR("Vector potential configuration not defined"); } diff --git a/AsterSeeds/src/make.code.defn b/AsterSeeds/src/make.code.defn index 96819da3..a0f10e36 100644 --- a/AsterSeeds/src/make.code.defn +++ b/AsterSeeds/src/make.code.defn @@ -9,6 +9,7 @@ SRCS = \ atmosphere.cxx \ bbhcloud.cxx \ initial_betafloor.cxx \ + perturbNS.cxx \ magTOV.cxx \ magBNS.cxx \ SetTempEntropyYe.cxx diff --git a/AsterSeeds/src/perturbNS.cxx b/AsterSeeds/src/perturbNS.cxx new file mode 100644 index 00000000..e38e5d96 --- /dev/null +++ b/AsterSeeds/src/perturbNS.cxx @@ -0,0 +1,177 @@ +#include + +#include +#include +#include + +#include +#include +#include + +#include "setup_eos.hxx" + +namespace AsterSeeds { +using namespace std; +using namespace Loop; +using namespace EOSX; + +enum class eos_3param { IdealGas, Hybrid, Tabulated }; +enum class ns_perturbation_t { m0, m1, m2 }; + +template +void AsterSeeds_PerturbNS_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_PerturbNS; + DECLARE_CCTK_PARAMETERS; + + ns_perturbation_t perturbation_mode; + + if (CCTK_EQUALS(ns_perturbation_mode, "m0")) { + perturbation_mode = ns_perturbation_t::m0; + if (ns_perturbation_strength < 0.0 || + ns_perturbation_strength >= 1.0) { + CCTK_VERROR( + "AsterSeeds::ns_perturbation_strength must be in [0, 1) for m0, " + "but is %.17g", + static_cast(ns_perturbation_strength)); + } + } else if (CCTK_EQUALS(ns_perturbation_mode, "m1")) { + perturbation_mode = ns_perturbation_t::m1; + } else if (CCTK_EQUALS(ns_perturbation_mode, "m2")) { + perturbation_mode = ns_perturbation_t::m2; + } else { + CCTK_VERROR("Unknown AsterSeeds::ns_perturbation_mode '%s'", + ns_perturbation_mode); + } + + if ((perturbation_mode == ns_perturbation_t::m1 || + perturbation_mode == ns_perturbation_t::m2) && + ns_perturbation_Req <= 0.0) { + CCTK_VERROR("AsterSeeds::ns_perturbation_Req must be positive, but is %.17g", + static_cast(ns_perturbation_Req)); + } + + grid.loop_all<1, 1, 1>( + grid.nghostzones, + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); + CCTK_REAL rho_atm = + (radial_distance > r_atmo) + ? (rho_abs_min * pow((r_atmo / radial_distance), n_rho_atmo)) + : rho_abs_min; + rho_atm = max(eos_3p->rgrho.min, rho_atm); + const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol); + + CCTK_REAL rhoL = rho(p.I); + if (rhoL <= rho_atmo_cut) { + return; + } + + CCTK_REAL pressL = press(p.I); + CCTK_REAL epsL = eps(p.I); + CCTK_REAL tempL = temperature(p.I); + CCTK_REAL entropyL = entropy(p.I); + const CCTK_REAL YeL = Ye(p.I); + + // Perturbations from Tsokaros et al., + // "Dynamically stable ergostars exist: General relativistic models and + // simulations", Phys. Rev. Lett. 123, 231103 (2019), + // arXiv:1907.03765. + switch (perturbation_mode) { + case ns_perturbation_t::m0: { + pressL *= (1.0 - ns_perturbation_strength); + tempL = eos_3p->temp_from_rho_press_ye(rhoL, pressL, YeL); + epsL = eos_3p->eps_from_rho_temp_ye(rhoL, tempL, YeL); + entropyL = eos_3p->kappa_from_rho_eps_ye(rhoL, epsL, YeL); + break; + } + case ns_perturbation_t::m1: { + const CCTK_REAL inv_Req = 1.0 / ns_perturbation_Req; + const CCTK_REAL factor = + 1.0 + ns_perturbation_strength * (p.x + p.y) * inv_Req; + rhoL = max(rhoL * factor, rho_atmo_cut); + pressL = eos_3p->press_from_rho_temp_ye(rhoL, tempL, YeL); + epsL = eos_3p->eps_from_rho_temp_ye(rhoL, tempL, YeL); + entropyL = eos_3p->kappa_from_rho_eps_ye(rhoL, epsL, YeL); + break; + } + case ns_perturbation_t::m2: { + const CCTK_REAL inv_Req = 1.0 / ns_perturbation_Req; + const CCTK_REAL factor = + 1.0 + ns_perturbation_strength * (p.x * p.x - p.y * p.y) * + inv_Req * inv_Req; + rhoL = max(rhoL * factor, rho_atmo_cut); + pressL = eos_3p->press_from_rho_temp_ye(rhoL, tempL, YeL); + epsL = eos_3p->eps_from_rho_temp_ye(rhoL, tempL, YeL); + entropyL = eos_3p->kappa_from_rho_eps_ye(rhoL, epsL, YeL); + break; + } + default: + assert(0); + } + + rho(p.I) = rhoL; + press(p.I) = pressL; + eps(p.I) = epsL; + temperature(p.I) = tempL; + entropy(p.I) = entropyL; + }); +} + +extern "C" void AsterSeeds_PerturbNS(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_PerturbNS; + DECLARE_CCTK_PARAMETERS; + + // defining EOS objects + eos_3param eos_3p_type; + + if (CCTK_EQUALS(evolution_eos, "IdealGas")) { + eos_3p_type = eos_3param::IdealGas; + } else if (CCTK_EQUALS(evolution_eos, "Hybrid")) { + eos_3p_type = eos_3param::Hybrid; + } else if (CCTK_EQUALS(evolution_eos, "Tabulated3d")) { + eos_3p_type = eos_3param::Tabulated; + } else { + CCTK_ERROR("Unknown value for parameter \"evolution_eos\""); + } + + switch (eos_3p_type) { + case eos_3param::IdealGas: { + // Get local eos object + auto eos_3p_ig = global_eos_3p_ig; + + AsterSeeds_PerturbNS_typeEoS(CCTK_PASS_CTOC, eos_3p_ig); + break; + } + case eos_3param::Hybrid: { + if (global_eos_3p_hyb_pwpoly) { + // pwpoly cold + pwpoly-hybrid + auto eos_3p_hyb = global_eos_3p_hyb_pwpoly; + + AsterSeeds_PerturbNS_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + + } else if (global_eos_3p_hyb_poly) { + // poly cold + poly-hybrid + auto eos_3p_hyb = global_eos_3p_hyb_poly; + + AsterSeeds_PerturbNS_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb); + + } else { + CCTK_ERROR( + "Hybrid EOS selected but no hybrid EOS object was initialized"); + } + + break; + } + case eos_3param::Tabulated: { + // Get local eos object + auto eos_3p_tab3d = global_eos_3p_tab3d; + + AsterSeeds_PerturbNS_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d); + break; + } + default: + assert(0); + } +} + +} // namespace AsterSeeds diff --git a/EOSX/scripts/CompOSE/README.md b/EOSX/scripts/CompOSE/README.md new file mode 100644 index 00000000..b7a5df08 --- /dev/null +++ b/EOSX/scripts/CompOSE/README.md @@ -0,0 +1,1019 @@ +# Generating CompOSE HS(DD2) EOS Tables + +This guide gives a complete, practical workflow for generating CompOSE EOS tables for the **HS(DD2) with electrons** EOS. + +It covers two outputs: + +1. a full 3D HDF5 table in `(T, n_b, Y_q)`, and +2. a cold beta-equilibrium slice. + +The same workflow can be adapted to other CompOSE EOS tables by changing the EOS download link and using the ranges from that EOS’s `eos.t`, `eos.nb`, and `eos.yq` files. + +--- + +## 1. What you will produce + +After the full 3D workflow: + +```text +eos.table.HSDD2.full3d.h5 +``` + +After the beta-equilibrium workflow: + +```text +eos.table.HSDD2.beta_T0p1MeV.dat +``` + +Optional one-column beta-equilibrium charge/electron fraction file: + +```text +Ye_beq_HSDD2.txt +``` + +--- + +## 2. Download the CompOSE code + +Create a working directory: + +```bash +mkdir -p ~/compose_work +cd ~/compose_work +``` + +Clone the CompOSE source code: + +```bash +git clone https://gitlab.obspm.fr/data_and_software_compose/code-compose.git +cd code-compose +``` + +If cloning with `git` is not available, download the source archive instead: + +```bash +mkdir -p ~/compose_work +cd ~/compose_work + +curl -L -o code-compose-master.tar.gz \ + https://gitlab.obspm.fr/data_and_software_compose/code-compose/-/archive/master/code-compose-master.tar.gz + +tar -xzf code-compose-master.tar.gz +cd code-compose-master +``` + +In the commands below, replace `~/compose_work/code-compose` with your actual CompOSE directory if it is different. + +--- + +## 3. Download the HS(DD2) EOS files + +The HS(DD2) CompOSE page is: + +```text +https://compose.obspm.fr/eos/18 +``` + +Download the EOS archive into the CompOSE source directory: + +```bash +cd ~/compose_work/code-compose + +wget -O hsdd2_eos.zip \ + 'https://compose.obspm.fr/download//3D/Hempel_SchaffnerBielich/hs_dd2_compose/with_electrons/eos.zip' +``` + +Unpack it: + +```bash +unzip hsdd2_eos.zip +``` + +You should now have files like: + +```text +eos.thermo +eos.compo +eos.micro +eos.t +eos.nb +eos.yq +eos.pdf +``` + +Check: + +```bash +ls -lh eos.* +``` + +--- + +## 4. Check the EOS grid ranges + +The files `eos.t`, `eos.nb`, and `eos.yq` contain the table grid. Use this command to print the number of points and the min/max values: + +```bash +for f in eos.t eos.nb eos.yq; do + awk 'NR==2{n=$1} NR==3{min=$1} NR>2{max=$1} + END{print FILENAME, "npoints =", n, "min =", min, "max =", max}' "$f" +done +``` + +For HS(DD2), you should get: + +```text +eos.t npoints = 81 min = 1.0000000E-01 max = 1.5848932E+02 +eos.nb npoints = 326 min = 1.0000000E-12 max = 1.0000000E+01 +eos.yq npoints = 60 min = 1.0000000E-02 max = 6.0000000E-01 +``` + +The first two lines of these files are metadata. The second line is the number of grid points, not a physical value. + +--- + +## 5. Build CompOSE + +### 5.1 Build with HDF5 support + +The full 3D table should be generated as HDF5. On a cluster, load an HDF5 module first. The module name depends on the machine. For example: + +```bash +module load cray-hdf5 +``` + +Then check that the HDF5 Fortran compiler wrapper is available: + +```bash +which h5fc +``` + +Build CompOSE with HDF5 enabled: + +```bash +cd ~/compose_work/code-compose + +unset USE_HDF5 +make clean +rm -f *.o *.mod compose + +make USE_HDF5=1 FC=h5fc compose 2>&1 | tee build_hdf5.log +``` + +Check that HDF5 was actually used: + +```bash +grep -E "h5fc|have_hdf5|modhdf5|hdf5compose" build_hdf5.log +``` + +You want to see lines containing: + +```text +h5fc +-Dhave_hdf5 +modhdf5.f90 +hdf5compose.f90 +``` + +If the build uses only `gfortran` and does not mention `modhdf5.f90` or `hdf5compose.f90`, then HDF5 was not enabled. + +### 5.2 Build without HDF5 + +For ASCII-only beta-equilibrium tables, a simple build can also work: + +```bash +make clean +make compose +``` + +But for the full 3D table used by an HDF5 reader, use the HDF5 build above. + +--- + +## 6. Generate the full 3D HDF5 table + +The full 3D table depends on: + +```text +T, n_b, Y_q +``` + +This is the table used by a 3D tabulated EOS reader. + +--- + +### 6.1 Generate `eos.quantities` + +Run: + +```bash +./compose +``` + +Choose task: + +```text +1 +``` + +If it asks whether to generate a new `eos.quantities`, enter: + +```text +1 +``` + +CompOSE expects numbers at prompts. Do not type `y` or `yes`. + +#### Regular thermodynamic quantities + +When asked: + +```text +How many regular thermodynamic quantities do you want to select for the file eos.table? +``` + +enter: + +```text +7 +``` + +Then enter these indices, one at a time: + +```text +1 +2 +3 +4 +5 +7 +12 +``` + +They mean: + +| Index | Quantity | +|---:|---| +| 1 | pressure `p` | +| 2 | entropy per baryon `S` | +| 3 | shifted baryon chemical potential `mu_b - m_n` | +| 4 | charge chemical potential `mu_q` | +| 5 | lepton chemical potential `mu_l` | +| 7 | scaled internal energy per baryon `E/m_n - 1` | +| 12 | sound speed squared `c_s^2` | + +#### Free-energy derivatives + +Enter: + +```text +0 +``` + +#### Composition particles + +CompOSE lists particle numbers and particle indices. For HS(DD2), the relevant entries are: + +| Particle number | Particle | Particle index | +|---:|---|---:| +| 2 | neutron | 10 | +| 3 | proton | 11 | +| 4 | alpha | 4002 | + +When asked: + +```text +How many particles do you want to select for the file eos.table? +``` + +enter: + +```text +3 +``` + +Then CompOSE asks for particle numbers. Enter: + +```text +2 +3 +4 +``` + +Do not enter `10`, `11`, and `4002` here. Those are particle indices, not particle numbers. + +The electron particle can be selected if you need it for diagnostics, but it is not needed for the usual reader setup because `Y_q` is already a table coordinate. + +#### Average nucleus set + +When asked: + +```text +How many sets do you want to select for the file eos.table? +``` + +enter: + +```text +1 +``` + +If CompOSE asks for the set number, enter: + +```text +1 +``` + +This selects the average nucleus set with index `999`. + +#### Microscopic quantities + +Enter: + +```text +0 +``` + +#### Error estimates + +Enter: + +```text +0 +``` + +#### Output format + +If CompOSE asks: + +```text +1: ASCII, else: HDF5 +``` + +enter: + +```text +2 +``` + +This selects HDF5 output. + +--- + +### 6.2 Generate `eos.parameters` + +Run: + +```bash +./compose +``` + +Choose task: + +```text +2 +``` + +#### Interpolation order + +When asked for the interpolation order in `T`, `n_b`, and `Y_q`, enter: + +```text +1 +1 +1 +``` + +#### Beta equilibrium + +When asked: + +```text +Please select if you want to calculate the EoS of matter in beta-equilibrium. +1: yes, else: no +``` + +enter: + +```text +0 +``` + +#### Fixed entropy + +When asked: + +```text +Please select if you want to calculate the EoS for given entropy per baryon. +1: yes, else: no +``` + +enter: + +```text +0 +``` + +#### Tabulation scheme + +When asked: + +```text +Please select the tabulation scheme for the parameters from +0: explicit listing of parameter values +1: loop form of parameter values +``` + +enter: + +```text +1 +``` + +#### Temperature grid + +Enter: + +```text +1.0000000E-01 1.5848932E+02 +81 +1 +``` + +This means: + +```text +T_min T_max = 1.0000000E-01 1.5848932E+02 +number of T points = 81 +logarithmic scaling = 1 +``` + +#### Baryon density grid + +Enter: + +```text +1.0000000E-12 1.0000000E+01 +326 +1 +``` + +This means: + +```text +n_b,min n_b,max = 1.0000000E-12 1.0000000E+01 +number of n_b points = 326 +logarithmic scaling = 1 +``` + +#### Charge fraction grid + +Enter: + +```text +1.0000000E-02 6.0000000E-01 +60 +0 +``` + +This means: + +```text +Y_q,min Y_q,max = 1.0000000E-02 6.0000000E-01 +number of Y_q points = 60 +linear scaling = 0 +``` + +CompOSE should now write: + +```text +eos.parameters +``` + +--- + +### 6.3 Generate the table + +Run: + +```bash +./compose +``` + +Choose task: + +```text +3 +``` + +A successful HDF5 run should include messages like: + +```text +format of output table = 2 (HDF5) +call writing HDF5 table +writing 7 thermodynamic quantities into file +writing 3 pairs into file +writing 1 quadruples into file +``` + +For HDF5 output, CompOSE usually writes the actual HDF5 file as: + +```text +eoscompose.h5 +``` + +Copy it to a clearer name: + +```bash +cp eoscompose.h5 eos.table.HSDD2.full3d.h5 +``` + +--- + +### 6.4 Check the HDF5 file + +Check that it is really HDF5: + +```bash +file eos.table.HSDD2.full3d.h5 +``` + +Expected: + +```text +Hierarchical Data Format (version 5) data +``` + +List the groups: + +```bash +h5ls eos.table.HSDD2.full3d.h5 +``` + +Expected groups include: + +```text +Composition_pairs +Composition_quadrupels +Parameters +Thermo_qty +metadata +``` + +Check the selected quantities: + +```bash +h5dump -d /Thermo_qty/index_thermo eos.table.HSDD2.full3d.h5 +h5dump -d /Composition_pairs/index_yi eos.table.HSDD2.full3d.h5 +h5dump -d /Composition_quadrupels/index_av eos.table.HSDD2.full3d.h5 +``` + +Expected values: + +```text +/Thermo_qty/index_thermo = 1, 2, 3, 4, 5, 7, 12 +/Composition_pairs/index_yi = 10, 11, 4002 +/Composition_quadrupels/index_av = 999 +``` + +--- + +## 7. Generate a cold beta-equilibrium slice + +The beta-equilibrium slice is a 1D table along `n_b`. It is useful for cold EOS checks and TOV-style preprocessing. + +For HS(DD2), use the lowest table temperature: + +```text +T = 1.0000000E-01 MeV +``` + +Before changing the CompOSE input files, save the full 3D setup: + +```bash +cp eos.quantities eos.quantities.full3d +cp eos.parameters eos.parameters.full3d +``` + +--- + +### 7.1 Generate beta-slice `eos.quantities` + +Run: + +```bash +./compose +``` + +Choose task: + +```text +1 +``` + +Regenerate `eos.quantities`: + +```text +1 +``` + +For a minimal beta-equilibrium table, select three regular thermodynamic quantities: + +```text +3 +21 +1 +3 +``` + +These are: + +| Index | Quantity | +|---:|---| +| 21 | internal energy per baryon `E` [MeV] | +| 1 | pressure `p` | +| 3 | shifted baryon chemical potential `mu_b - m_n` | + +For the remaining groups, enter: + +```text +0 +0 +0 +0 +0 +``` + +That means: + +```text +free-energy derivatives: 0 +composition particles: 0 +average sets: 0 +microscopic quantities: 0 +error estimates: 0 +``` + +For a simple beta-equilibrium text table, choose ASCII output: + +```text +1 +``` + +If you want extra diagnostics, choose six regular quantities instead: + +```text +6 +21 +1 +3 +12 +15 +2 +``` + +This gives: + +```text +E, p, mu_b - m_n, c_s^2, Gamma, entropy +``` + +--- + +### 7.2 Generate beta-slice `eos.parameters` + +Run: + +```bash +./compose +``` + +Choose task: + +```text +2 +``` + +Use interpolation order 1 in all variables: + +```text +1 +1 +1 +``` + +Turn on beta equilibrium: + +```text +1 +``` + +Turn off fixed entropy: + +```text +0 +``` + +Use loop form: + +```text +1 +``` + +Use one temperature point at the lowest table temperature: + +```text +1.0000000E-01 1.0000000E-01 +1 +0 +``` + +Use the full baryon density range: + +```text +1.0000000E-12 1.0000000E+01 +326 +1 +``` + +Use the full charge-fraction range as the search range for the beta-equilibrium solve: + +```text +1.0000000E-02 6.0000000E-01 +``` + +If CompOSE asks for a number of `Y_q` points and scaling, enter: + +```text +60 +0 +``` + +If it does not ask for those two values in beta-equilibrium mode, that is fine. + +--- + +### 7.3 Generate the beta-equilibrium table + +Run: + +```bash +./compose +``` + +Choose task: + +```text +3 +``` + +Rename the output immediately: + +```bash +mv eos.table eos.table.HSDD2.beta_T0p1MeV.dat +cp eos.beta eos.beta.HSDD2_T0p1MeV.dat 2>/dev/null || true +cp eos.report eos.report.HSDD2.beta_T0p1MeV.txt 2>/dev/null || true +``` + +Check: + +```bash +head eos.table.HSDD2.beta_T0p1MeV.dat +wc -l eos.table.HSDD2.beta_T0p1MeV.dat +``` + +For the minimal quantity selection, the columns are: + +```text +1 T [MeV] +2 n_b [fm^-3] +3 Y_q +4 E [MeV] +5 p [MeV fm^-3] +6 mu_b - m_n [MeV] +``` + +--- + +## 8. Create a one-column beta-equilibrium `Y_e` file + +For charge-neutral matter in this table, the beta-equilibrium `Y_e`/charge fraction is the third column of the beta table. + +Create a one-column file: + +```bash +awk '!/^#/ {printf "%.18e\n", $3}' \ + eos.table.HSDD2.beta_T0p1MeV.dat > Ye_beq_HSDD2.txt +``` + +Check it: + +```bash +head Ye_beq_HSDD2.txt +wc -l Ye_beq_HSDD2.txt +``` + +For the full HS(DD2) `n_b` grid, the file should have: + +```text +326 +``` + +lines. + +--- + +## 9. Optional: extract a beta slice from the full 3D HDF5 table + +If you already have the full 3D HDF5 table, you can also post-process it by solving: + +```text +mu_l(T, n_b, Y_q) = 0 +``` + +along the `Y_q` direction. + +If you have the helper script `generate_beta_eq_slice.py`, run: + +```bash +python3 generate_beta_eq_slice.py eos.table.HSDD2.full3d.h5 \ + --cold \ + --out-prefix HSDD2_beta_cold +``` + +Then make the one-column file: + +```bash +awk '!/^#/ {printf "%.18e\n", $3}' HSDD2_beta_cold.dat > Ye_beq_HSDD2.txt +``` + +This method is useful for quick post-processing. The built-in CompOSE beta-equilibrium mode is the preferred method for production beta slices. + +--- + +## 10. Restore the full 3D setup + +If you saved the full 3D input files and want to regenerate the full 3D table later: + +```bash +cp eos.quantities.full3d eos.quantities +cp eos.parameters.full3d eos.parameters + +./compose +# choose task 3 +``` + +--- + +## 11. Notes for EOS readers + +A typical HDF5 reader should use: + +```text +/Parameters +/Thermo_qty +/Composition_pairs +/Composition_quadrupels +``` + +The full 3D table generated above contains: + +```text +index_thermo = 1, 2, 3, 4, 5, 7, 12 +index_yi = 10, 11, 4002 +index_av = 999 +``` + +For chemical potentials: + +```text +3 = mu_b - m_n +4 = mu_q +5 = mu_l +``` + +A common mapping is: + +```text +mu_n = mu_b - m_n +mu_p = mu_n + mu_q +mu_e = mu_l - mu_q +``` + +If your reader has a variable for `mu_l` or `MUNU`, fill it from CompOSE index `5`. + +--- + +## 12. Quick troubleshooting + +### The generated `.h5` file is not HDF5 + +Check: + +```bash +file eos.table.HSDD2.full3d.h5 +``` + +If it says `ASCII text`, the file is ASCII even if it has a `.h5` name. Regenerate with HDF5 output selected. + +A real HDF5 file should say: + +```text +Hierarchical Data Format (version 5) data +``` + +### HDF5 output says it wrote `eos.table`, but there is no `eos.table` + +For HDF5 output, the actual file is usually: + +```text +eoscompose.h5 +``` + +Use: + +```bash +ls -lh eoscompose.h5 +cp eoscompose.h5 eos.table.HSDD2.full3d.h5 +``` + +### The build did not enable HDF5 + +Check the build log: + +```bash +grep -E "h5fc|have_hdf5|modhdf5|hdf5compose" build_hdf5.log +``` + +If those strings are missing, rebuild with: + +```bash +make clean +rm -f *.o *.mod compose +make USE_HDF5=1 FC=h5fc compose 2>&1 | tee build_hdf5.log +``` + +### CompOSE crashes after typing `yes` + +CompOSE prompts expect numbers, not words. Use: + +```text +1 +``` + +instead of: + +```text +yes +``` + +### Files get overwritten + +CompOSE reuses filenames such as: + +```text +eos.quantities +eos.parameters +eos.table +eoscompose.h5 +eos.report +``` + +Save important files before switching workflows: + +```bash +cp eos.quantities eos.quantities.full3d +cp eos.parameters eos.parameters.full3d +cp eoscompose.h5 eos.table.HSDD2.full3d.h5 +``` + +--- + +## 13. Final checklist + +For the full 3D table: + +```bash +file eos.table.HSDD2.full3d.h5 +h5ls eos.table.HSDD2.full3d.h5 +h5dump -d /Thermo_qty/index_thermo eos.table.HSDD2.full3d.h5 +h5dump -d /Composition_pairs/index_yi eos.table.HSDD2.full3d.h5 +h5dump -d /Composition_quadrupels/index_av eos.table.HSDD2.full3d.h5 +``` + +Expected: + +```text +index_thermo = 1, 2, 3, 4, 5, 7, 12 +index_yi = 10, 11, 4002 +index_av = 999 +``` + +For the beta-equilibrium slice: + +```bash +head eos.table.HSDD2.beta_T0p1MeV.dat +awk '!/^#/ {printf "%.18e\n", $3}' \ + eos.table.HSDD2.beta_T0p1MeV.dat > Ye_beq_HSDD2.txt +wc -l Ye_beq_HSDD2.txt +``` diff --git a/EOSX/scripts/CompOSE/generate_beta_eq_slice.py b/EOSX/scripts/CompOSE/generate_beta_eq_slice.py new file mode 100644 index 00000000..798db18d --- /dev/null +++ b/EOSX/scripts/CompOSE/generate_beta_eq_slice.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python3 +""" +make_beta_slice_from_compose_h5.py + +Extract a beta-equilibrium slice from an already-created full 3D CompOSE HDF5 +table by solving + + mu_l(T, n_b, Y_q) = 0 + +along the Y_q direction at each (T, n_b). + +For the HS(DD2) HDF5 table generated with the selection + + index_thermo = 1, 2, 3, 4, 5, 7, 12 + +CompOSE thermodynamic index 5 is the lepton chemical potential mu_l. For +neutrino-less beta equilibrium, mu_l = 0 is the condition used here. + +Example usage: + + python make_beta_slice_from_compose_h5.py eos.table.HSDD2.h5 \ + --out-prefix HSDD2_beta_allT + +Lowest-temperature/cold slice only: + + python make_beta_slice_from_compose_h5.py eos.table.HSDD2.h5 \ + --cold --out-prefix HSDD2_beta_cold + +A specific temperature index: + + python make_beta_slice_from_compose_h5.py eos.table.HSDD2.h5 \ + --temperature-index 0 --out-prefix HSDD2_beta_T0 + +Outputs: + + .dat + ASCII table with columns: + T, nb, Yq_beta, status, + all selected thermodynamic quantities, + optional composition pair quantities, + optional average-nucleus quantities. + + .h5 + Compact postprocessed HDF5 file with /Beta_slice/table. + This is NOT the same 3D CompOSE layout as the input file. + +status column: + 0 = mu_l=0 crossing was bracketed and linearly interpolated. + 1 = no crossing was bracketed; nearest |mu_l| point was used. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path +from typing import Dict, List, Tuple + +import h5py +import numpy as np + + +THERMO_NAMES = { + 1: "press_MeV_fm^-3", + 2: "entropy", + 3: "mu_b_minus_m_n_MeV", + 4: "mu_q_MeV", + 5: "mu_l_MeV", + 7: "eps_scaled_E_over_mn_minus_1", + 12: "cs2", + 15: "Gamma", + 21: "E_MeV", +} + + +def map_axes_by_unique_size(shape: Tuple[int, ...], sizes: Dict[str, int]) -> Dict[str, int]: + """ + Map raw HDF5 dataset axes to nb/t/yq/q by matching dimension sizes. + + For HS(DD2), the dimensions are usually unique: + nb=326, t=81, yq=60, q=7 or 3 or 1. + """ + axes: Dict[str, int] = {} + used = set() + + for name, size in sizes.items(): + matches = [i for i, s in enumerate(shape) if s == size and i not in used] + if len(matches) != 1: + raise RuntimeError( + f"Could not uniquely identify axis {name} of size {size}. " + f"Dataset shape is {shape}; matches={matches}. " + "If two dimensions have the same size, edit the transpose logic manually." + ) + axes[name] = matches[0] + used.add(matches[0]) + + return axes + + +def to_nb_t_yq_q(raw: np.ndarray, nrho: int, ntemp: int, nye: int, nq: int) -> np.ndarray: + axes = map_axes_by_unique_size(raw.shape, {"nb": nrho, "t": ntemp, "yq": nye, "q": nq}) + return np.transpose(raw, (axes["nb"], axes["t"], axes["yq"], axes["q"])) + + +def interp_along_yq_to_mu_target( + yq: np.ndarray, + values: np.ndarray, + mu_l: np.ndarray, + target_mu_l: float, +) -> Tuple[np.ndarray, float, int]: + """ + Interpolate all values to mu_l = target_mu_l along Y_q. + + values has shape (nye, ncols). + """ + f = mu_l - target_mu_l + finite = np.isfinite(yq) & np.isfinite(f) & np.all(np.isfinite(values), axis=1) + + if np.count_nonzero(finite) < 2: + return np.full(values.shape[1], np.nan), np.nan, 1 + + y = yq[finite] + ff = f[finite] + vv = values[finite, :] + + # Exact/nearly exact point. + iz = int(np.argmin(np.abs(ff))) + if abs(ff[iz]) < 1.0e-12: + return vv[iz].copy(), float(y[iz]), 0 + + # Locate sign changes. + changes = np.where(ff[:-1] * ff[1:] <= 0.0)[0] + + if len(changes) == 0: + # No bracket; use nearest |mu_l|. + return vv[iz].copy(), float(y[iz]), 1 + + # If multiple crossings exist, choose the local bracket closest to zero. + i0 = int(min(changes, key=lambda i: min(abs(ff[i]), abs(ff[i + 1])))) + i1 = i0 + 1 + + f0 = ff[i0] + f1 = ff[i1] + if f1 == f0: + w = 0.0 + else: + w = -f0 / (f1 - f0) + + w = float(np.clip(w, 0.0, 1.0)) + y_beta = y[i0] + w * (y[i1] - y[i0]) + out = vv[i0, :] + w * (vv[i1, :] - vv[i0, :]) + + return out, float(y_beta), 0 + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("input_h5", help="Full 3D CompOSE HDF5 table") + parser.add_argument("--out-prefix", default="beta_slice", help="Output prefix") + parser.add_argument("--target-mul", type=float, default=0.0, help="Target mu_l in MeV") + parser.add_argument("--cold", action="store_true", help="Only extract the lowest-temperature slice") + parser.add_argument("--temperature-index", type=int, default=None, help="Extract one temperature index") + args = parser.parse_args() + + input_h5 = Path(args.input_h5) + out_prefix = Path(args.out_prefix) + + with h5py.File(input_h5, "r") as f: + nb = np.asarray(f["/Parameters/nb"][:], dtype=float) + temp = np.asarray(f["/Parameters/t"][:], dtype=float) + yq = np.asarray(f["/Parameters/yq"][:], dtype=float) + + index_thermo = np.asarray(f["/Thermo_qty/index_thermo"][:], dtype=int) + thermo_raw = np.asarray(f["/Thermo_qty/thermo"][:], dtype=float) + + nrho = len(nb) + ntemp = len(temp) + nye = len(yq) + nthermo = len(index_thermo) + + thermo = to_nb_t_yq_q(thermo_raw, nrho, ntemp, nye, nthermo) + + if 5 not in index_thermo: + raise RuntimeError("Thermodynamic index 5, mu_l, is not present.") + i_mul = int(np.where(index_thermo == 5)[0][0]) + + comp_pair_indices = np.array([], dtype=int) + comp_pairs = None + if "/Composition_pairs" in f: + comp_pair_indices = np.asarray(f["/Composition_pairs/index_yi"][:], dtype=int) + yi_raw = np.asarray(f["/Composition_pairs/yi"][:], dtype=float) + comp_pairs = to_nb_t_yq_q(yi_raw, nrho, ntemp, nye, len(comp_pair_indices)) + + av_indices = np.array([], dtype=int) + av_data = None + av_group_name = None + if "/Composition_quadrupels" in f: + av_group_name = "/Composition_quadrupels" + elif "/Composition_quadruples" in f: + av_group_name = "/Composition_quadruples" + + if av_group_name is not None: + g = f[av_group_name] + av_indices = np.asarray(g["index_av"][:], dtype=int) + navsets = len(av_indices) + + yav = to_nb_t_yq_q(np.asarray(g["yav"][:], dtype=float), nrho, ntemp, nye, navsets) + aav = to_nb_t_yq_q(np.asarray(g["aav"][:], dtype=float), nrho, ntemp, nye, navsets) + zav = to_nb_t_yq_q(np.asarray(g["zav"][:], dtype=float), nrho, ntemp, nye, navsets) + nav = to_nb_t_yq_q(np.asarray(g["nav"][:], dtype=float), nrho, ntemp, nye, navsets) + + # Interleave by set: Yav, Aav, Zav, Nav for each set. + av_parts = [] + for iset in range(navsets): + av_parts.extend( + [ + yav[:, :, :, iset:iset+1], + aav[:, :, :, iset:iset+1], + zav[:, :, :, iset:iset+1], + nav[:, :, :, iset:iset+1], + ] + ) + av_data = np.concatenate(av_parts, axis=3) + + if args.temperature_index is not None: + t_indices = [args.temperature_index] + elif args.cold: + t_indices = [0] + else: + t_indices = list(range(ntemp)) + + for jt in t_indices: + if jt < 0 or jt >= ntemp: + raise RuntimeError(f"Temperature index {jt} is outside [0, {ntemp - 1}]") + + thermo_col_names = [THERMO_NAMES.get(int(idx), f"thermo_{int(idx)}") for idx in index_thermo] + comp_col_names: List[str] = [f"Y_particle_{int(idx)}" for idx in comp_pair_indices] + + av_col_names: List[str] = [] + for idx in av_indices: + av_col_names.extend( + [ + f"Yav_set_{int(idx)}", + f"Aav_set_{int(idx)}", + f"Zav_set_{int(idx)}", + f"Nav_set_{int(idx)}", + ] + ) + + columns = ["T_MeV", "nb_fm^-3", "Yq_beta", "status"] + columns += thermo_col_names + comp_col_names + av_col_names + + rows = [] + n_cross = 0 + n_fallback = 0 + + for jt in t_indices: + for ir in range(nrho): + pieces = [thermo[ir, jt, :, :]] + + if comp_pairs is not None: + pieces.append(comp_pairs[ir, jt, :, :]) + + if av_data is not None: + pieces.append(av_data[ir, jt, :, :]) + + values = np.concatenate(pieces, axis=1) + mu_l = thermo[ir, jt, :, i_mul] + + out_values, yq_beta, status = interp_along_yq_to_mu_target( + yq, values, mu_l, args.target_mul + ) + + if status == 0: + n_cross += 1 + else: + n_fallback += 1 + + row = [temp[jt], nb[ir], yq_beta, float(status)] + row.extend(out_values.tolist()) + rows.append(row) + + table = np.asarray(rows, dtype=float) + + dat_path = out_prefix.with_suffix(".dat") + h5_path = out_prefix.with_suffix(".h5") + + np.savetxt(dat_path, table, header=" ".join(columns), fmt="%.16e") + + with h5py.File(h5_path, "w") as fout: + fout.attrs["source_file"] = str(input_h5) + fout.attrs["beta_condition"] = "mu_l = target_mu_l" + fout.attrs["target_mu_l_MeV"] = args.target_mul + fout.attrs["status_meaning"] = "0=crossing found; 1=no crossing, nearest |mu_l| used" + + gp = fout.create_group("Parameters") + gp.create_dataset("nb", data=nb) + gp.create_dataset("t", data=temp[t_indices]) + gp.create_dataset("columns", data=np.asarray(columns, dtype="S")) + + gb = fout.create_group("Beta_slice") + gb.create_dataset("table", data=table) + gb.create_dataset("index_thermo", data=index_thermo) + if comp_pair_indices.size: + gb.create_dataset("index_yi", data=comp_pair_indices) + if av_indices.size: + gb.create_dataset("index_av", data=av_indices) + + print(f"Wrote {dat_path}") + print(f"Wrote {h5_path}") + print(f"Rows: {len(rows)}") + print(f"Crossings found: {n_cross}") + print(f"Fallback/min-|mu_l| points: {n_fallback}") + if n_fallback: + print("WARNING: Some points did not bracket mu_l=0 over the available Yq range.") + + +if __name__ == "__main__": + main() diff --git a/EOSX/scripts/README.md b/EOSX/scripts/StellarCollapse/README.md similarity index 100% rename from EOSX/scripts/README.md rename to EOSX/scripts/StellarCollapse/README.md diff --git a/EOSX/scripts/adjust_table_mineps0.py b/EOSX/scripts/StellarCollapse/adjust_table_mineps0.py similarity index 100% rename from EOSX/scripts/adjust_table_mineps0.py rename to EOSX/scripts/StellarCollapse/adjust_table_mineps0.py diff --git a/EOSX/scripts/extend_table.py b/EOSX/scripts/StellarCollapse/extend_table.py similarity index 100% rename from EOSX/scripts/extend_table.py rename to EOSX/scripts/StellarCollapse/extend_table.py diff --git a/EOSX/scripts/generate_Ye_of_rho.py b/EOSX/scripts/StellarCollapse/generate_Ye_of_rho.py similarity index 100% rename from EOSX/scripts/generate_Ye_of_rho.py rename to EOSX/scripts/StellarCollapse/generate_Ye_of_rho.py diff --git a/EOSX/scripts/soften_table.py b/EOSX/scripts/StellarCollapse/soften_table.py similarity index 100% rename from EOSX/scripts/soften_table.py rename to EOSX/scripts/StellarCollapse/soften_table.py diff --git a/EOSX/src/eos_3p_tabulated3d/eos_readtable_compose.hxx b/EOSX/src/eos_3p_tabulated3d/eos_readtable_compose.hxx index 4cca9dbf..c32ca37f 100644 --- a/EOSX/src/eos_3p_tabulated3d/eos_readtable_compose.hxx +++ b/EOSX/src/eos_3p_tabulated3d/eos_readtable_compose.hxx @@ -336,10 +336,10 @@ CCTK_HOST inline void eos_readtable_compose(const std::string &filename, } int nav = 0; - if (hdf5_link_exists(file, "/Composition_quadruples")) { + if (hdf5_link_exists(file, "/Composition_quadrupels")) { hid_t av_id; HDF5_ERROR(av_id = - H5Gopen2(file, "/Composition_quadruples", H5P_DEFAULT)); + H5Gopen2(file, "/Composition_quadrupels", H5P_DEFAULT)); READ_ATTR_HDF5_COMPOSE(av_id, "pointsav", &nav, H5T_NATIVE_INT); if (nav > 0) {