Release version 1.6

* Breaking change in spline EOS file format
  EOS of type eos_barotr_spline are now interpolating the soundspeed as
  function of density instead pseudo-enthalpy, because otherwise
  phase transitions cannot be represented. This required a change in
  the file format for this EOS type as well. Old files of type
  eos_barotr_spline cannot be read anymore. Files of all other EOS types
  are unaffected.
* The Python API can now be build also for Python 3.11

EOS/create/make_read_pp.py

@@ -64,7 +64,7 @@ def make_read_eos(name, p1_cgs, ga, rho_max_si=1e20):
 
     path     = f"{name}.eos.h5"
 
-    uc = pyr.units.geom_solar
+    uc = pyr.units.geom_solar()
     rho_poly = sly_ds_si[0] / uc.density
     rho_b    = rho_b_si / uc.density 
     rho_max = rho_max_si / uc.density

bindings/python/pyproject.toml

@@ -1,18 +1,18 @@
 [build-system]
-requires = ["setuptools>=42", "wheel", "pybind11~=2.6.1"]
+requires = ["setuptools>=42", "wheel", "pybind11>=2.6.1"]
 build-backend = "setuptools.build_meta"
 
 
 [project]
   name = "RePrimAnd"
-  version = "1.5"
+  version = "1.6"
   authors = [
     { name = "Wolfgang Kastaun", email = "[email protected]" }
   ]  
   description = "Library for handling equations of state for supranuclear matter, computing neutron star properties, and utilities for numerical relativity"
   readme = "README.md" 
   license = { file = "LICENSE" }
-  requires-python = '>=3.7, <3.11'
+  requires-python = '>=3.7'
   dependencies = [
     "numpy"
   ]
@@ -26,7 +26,6 @@ build-backend = "setuptools.build_meta"
 [tool.cibuildwheel]
 
   archs = ["x86_64"]
-  skip="cp311-manylinux_x86_64 cp311-musllinux_x86_64"
 
   before-all = "/project/docker/cibw_prepare.sh"
 

docsrc/sphinx/eos_barotr_available.rst

@@ -126,24 +126,24 @@ This EOS implements all functions using monotonic cubic spline
 interpolation. The EOS is therefore differentiable in principle. 
 Of course, there can still be steep gradients.
 
-Internally, most properties are represented in terms of the 
-pseudo enthalpy. This has some advantages for computing hydrostatic
-equilibrium models and also with regard to phase 
-transitions (although the density has a jump at those).
+Internally, most properties are internally represented as functions 
+of the pseudo enthalpy. This has some advantages for computing hydrostatic
+equilibrium models and with regard to phase transitions. 
 When calling the EOS using density as independent variable, another 
 interpolation spline is used to first compute the pseudo enthalpy
 (in presence of phase transitions, it has a plateau as function of 
 density). The desired quantity is then computed from the pseudo enthalpy
 using the same interpolation splines used for evaluating the EOS as 
-function of pseudo enthalpy.
-The monotonic interpolation ensures that the EOS does not produce 
-unphysical overshoots.
+function of pseudo enthalpy. On exception is the soundspeed, which is 
+internally represented as function of density (this is because at 
+phase transition it has infinitly sharp features as function of 
+pseudo enthalpy). The monotonic interpolation ensures that the EOS 
+does not produce unphysical overshoots.
 
 The spline sample points are spaced regularly with respect to 
 logarithm of pseudo-enthalpy-minus-one  :math:`\log(g-1)` or 
-mass density :math:`\rho`. This allows 
-efficient computation with cost nearly independent of the sample 
-resolution.
+mass density :math:`\rho`. The regular sampling allows efficient 
+computation with cost nearly independent of the sample resolution.
 In order to use the number of sample points efficiently, the spline 
 interpolation covers a user-specified range of magnitudes. Below,
 a generalized polytrope (meaning an additional offset in specific
@@ -154,7 +154,12 @@ If the temperature is not provided when creating the EOS, it is assumed
 to be a zero-temperature EOS.
 
 To set up this type of EOS, one provides individual functions which are then
-sampled to create the interpolation splines. In addition, one has to specify 
+sampled to create the interpolation splines. There are different options
+which quantities need to be provided. One always needs density and pressure.
+Providing the pseudo-enthalpy is optional, it can be recomputed from density, 
+energy, and pressure. Providing the specific energy is optional for 
+isentropic EOS, where it can be recomputed from pressure and density.
+In addition, one has to specify 
 the maximum validity range, the matching point to the polytrope, and its 
 exponent. Note that this polytrope is completely determined by the density,
 energy, and pressure at the matching point. Since the pseudo-enthalpy is an integral

library/EOS_Barotropic/eos_barotr_spline.cc

@@ -15,16 +15,15 @@ using namespace EOS_Toolkit::implementations;
 
 using func_t  = std::function<real_t(real_t)>;
 
-auto eos_barotr_spline::get_rggm1(const lglgspl_t& gm1_rho, 
-  const lgspl_t& eps_gm1, const lglgspl_t& p_gm1, 
-  const lgspl_t& hm1_gm1, const lglgspl_t& rho_gm1, 
-  const lgspl_t& csnd_gm1, const opt_t& temp_gm1, 
+auto eos_barotr_spline::get_rggm1(const lgspl_t& eps_gm1, 
+  const lglgspl_t& p_gm1, const lgspl_t& hm1_gm1, 
+  const lglgspl_t& rho_gm1, const opt_t& temp_gm1, 
   const opt_t& efrac_gm1)
 -> range
 {
   range rg{ 
     intersect(eps_gm1.range_x(), p_gm1.range_x(), hm1_gm1.range_x(),
-              rho_gm1.range_x(), csnd_gm1.range_x())
+              rho_gm1.range_x())
   };
 
   if (temp_gm1) 
@@ -43,7 +42,7 @@ auto eos_barotr_spline::get_rggm1(const lglgspl_t& gm1_rho,
 eos_barotr_spline::eos_barotr_spline(
     lglgspl_t gm1_rho_, lglgspl_t rho_gm1_,
     lgspl_t eps_gm1_, lglgspl_t p_gm1_, 
-    lgspl_t hm1_gm1_, lgspl_t csnd_gm1_, 
+    lgspl_t hm1_gm1_, lgspl_t csnd_rho_, 
     opt_t temp_gm1_, opt_t efrac_gm1_, 
     bool isentropic_,
     const eos_barotr_gpoly& poly_)
@@ -53,12 +52,12 @@ eos_barotr_spline::eos_barotr_spline(
   p_gm1{ std::move(p_gm1_) }, 
   hm1_gm1{ std::move(hm1_gm1_) },
   rho_gm1{ std::move(rho_gm1_) }, 
-  csnd_gm1{ std::move(csnd_gm1_) }, 
+  csnd_rho{ std::move(csnd_rho_) }, 
   temp_gm1{ std::move(temp_gm1_) }, 
   efrac_gm1{ std::move(efrac_gm1_) }, 
   poly{ poly_ },  
-  rggm1{ get_rggm1(gm1_rho, eps_gm1, p_gm1, hm1_gm1, rho_gm1,
-                   csnd_gm1,temp_gm1, efrac_gm1) },
+  rggm1{ get_rggm1(eps_gm1, p_gm1, hm1_gm1, rho_gm1,
+                   temp_gm1, efrac_gm1) },
   rgrho{ 0, rho_gm1(rggm1.max()) },
   gm1_low{ poly.range_gm1().max() },
   rho_low{ poly.range_rho().max() },
@@ -69,7 +68,6 @@ eos_barotr_spline::eos_barotr_spline(
   if (!(eps_gm1.contains(gm1_low) &&
         p_gm1.contains(gm1_low) &&
         rho_gm1.contains(gm1_low) &&
-        csnd_gm1.contains(gm1_low) &&
         hm1_gm1.contains(gm1_low) &&
         (!temp_gm1 || temp_gm1->contains(gm1_low)) &&
         (!efrac_gm1 || efrac_gm1->contains(gm1_low)) ))
@@ -79,7 +77,7 @@ eos_barotr_spline::eos_barotr_spline(
   }
 
 
-  if (!gm1_rho.contains(rho_low))
+  if (!(gm1_rho.contains(rho_low) && csnd_rho.contains(rho_low)))
   {
     throw runtime_error("eos_barotr_spline: matching polytrope "
                         "outside sampled range for rho");     
@@ -89,10 +87,10 @@ eos_barotr_spline::eos_barotr_spline(
     throw runtime_error("eos_barotr_spline: negative mass density "
                         "in rho(gm1)");
   }
-  if (csnd_gm1.range_y().max() >= 1.0) {
+  if (csnd_rho.range_y().max() >= 1.0) {
     throw runtime_error("eos_barotr_spline: sound speed >= 1");
   }
-  if (csnd_gm1.range_y().min() < 0.0) {
+  if (csnd_rho.range_y().min() < 0.0) {
     throw runtime_error("eos_barotr_spline: sound speed < 0");
   }
   if (p_gm1.range_y().min() < 0.0) {
@@ -150,7 +148,12 @@ real_t eos_barotr_spline::hm1(real_t gm1) const
 
 real_t eos_barotr_spline::csnd(real_t gm1) const
 {
-  return (gm1 >= gm1_low) ? csnd_gm1(gm1) : poly.csnd(gm1);
+  return (gm1 >= gm1_low) ? csnd_rho(rho_gm1(gm1)) : poly.csnd(gm1);
+}
+
+real_t eos_barotr_spline::csnd_from_rho_gm1(real_t rho, real_t gm1) const 
+{
+  return (rho >= rho_low) ? csnd_rho(rho) : poly.csnd(gm1);
 }
 
 real_t eos_barotr_spline::temp(real_t gm1) const
@@ -205,7 +208,7 @@ auto eos_barotr_spline::descr_str() const -> std::string
 
 auto EOS_Toolkit::make_eos_barotr_spline(
   func_t gm1_rho, func_t rho_gm1, func_t eps_gm1, func_t press_gm1, 
-  func_t csnd_gm1, func_t temp_gm1, func_t efrac_gm1, 
+  func_t csnd_rho, func_t temp_gm1, func_t efrac_gm1, 
   bool isentropic, interval<real_t> rg_rho, real_t n_poly,
   units u, std::size_t pts_per_mag)
 -> eos_barotr
@@ -287,9 +290,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(
   };
 
   auto scsnd{ 
-    interpol_logspl_impl::from_function(
-      [&] (real_t gm1) -> real_t {return csnd_gm1(gm1_old(gm1));},
-      rg_gm1, npts_gm1)
+    interpol_logspl_impl::from_function(csnd_rho, rg_rho, npts_rho)
   };
 
   boost::optional<interpol_logspl_impl> stemp;
@@ -340,7 +341,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(const eos_barotr& eos,
             [&eos](real_t gm1) {return eos.at_gm1(gm1).rho();},  
             [&eos](real_t gm1) {return eos.at_gm1(gm1).eps();},  
             [&eos](real_t gm1) {return eos.at_gm1(gm1).press();},  
-            [&eos](real_t gm1) {return eos.at_gm1(gm1).csnd();},  
+            [&eos](real_t rho) {return eos.at_rho(rho).csnd();},  
             temp_gm1, efrac_gm1, 
             eos.is_isentropic(), rg_rho, n_poly,
             eos.units_to_SI(), pts_per_mag
@@ -442,7 +443,6 @@ auto EOS_Toolkit::make_eos_barotr_spline(
 
   auto eps_gm1 = [&] (real_t gm1) {return eps_rho(rho_gm1(gm1));};
   auto press_gm1 = [&] (real_t gm1) {return press_rho(rho_gm1(gm1));};
-  auto csnd_gm1 = [&] (real_t gm1) {return csnd_rho(rho_gm1(gm1));};
 
   auto temp_gm1 = [&]() -> func_t {
     if (temp.empty()) return nullptr;
@@ -463,7 +463,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(
   }
 
   return make_eos_barotr_spline(gm1_rho, rho_gm1, eps_gm1, 
-                      press_gm1, csnd_gm1, temp_gm1, efrac_gm1, 
+                      press_gm1, csnd_rho, temp_gm1, efrac_gm1, 
                       isentropic, rg_rho, n_poly, uc, pts_per_mag);
 }
 
@@ -522,7 +522,6 @@ auto EOS_Toolkit::make_eos_barotr_spline(
 
   auto eps_gm1 = [&] (real_t gm1) {return eps_rho(rho_gm1(gm1));};
   auto press_gm1 = [&] (real_t gm1) {return press_rho(rho_gm1(gm1));};
-  auto csnd_gm1 = [&] (real_t gm1) {return csnd_rho(rho_gm1(gm1));};
 
   auto temp_gm1 = [&]() -> func_t {
     if (temp.empty()) return nullptr;
@@ -543,7 +542,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(
   }
 
   return make_eos_barotr_spline(gm1_rho, rho_gm1, eps_gm1, 
-                      press_gm1, csnd_gm1, temp_gm1, efrac_gm1, 
+                      press_gm1, csnd_rho, temp_gm1, efrac_gm1, 
                       true, rg_rho, n_poly, uc, pts_per_mag);
 
 }
@@ -562,7 +561,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(
   auto rho_gm1{ make_interpol_pchip_spline(gm1, rho) };
   auto eps_gm1{ make_interpol_pchip_spline(gm1, eps) };
   auto press_gm1{ make_interpol_pchip_spline(gm1, press) };
-  auto csnd_gm1{ make_interpol_pchip_spline(gm1, csnd) };
+  auto csnd_rho{ make_interpol_pchip_spline(rho, csnd) };
 
   func_t temp_gm1{nullptr};
   if (!temp.empty()) 
@@ -583,7 +582,7 @@ auto EOS_Toolkit::make_eos_barotr_spline(
   }
 
   return make_eos_barotr_spline(gm1_rho, rho_gm1, eps_gm1, 
-                      press_gm1, csnd_gm1, temp_gm1, efrac_gm1, 
+                      press_gm1, csnd_rho, temp_gm1, efrac_gm1, 
                       isentropic, rg_rho, n_poly, uc, pts_per_mag);
 
 }

library/EOS_Barotropic/eos_barotr_spline_file.cc

@@ -7,7 +7,7 @@
 namespace EOS_Toolkit {
 namespace implementations {
 
-const std::string eos_barotr_spline::datastore_id{ "barotr_spline" };
+const std::string eos_barotr_spline::datastore_id{ "barotr_spline_v2" };
 
 struct reader_eos_barotr_spline : reader_eos_barotr 
 {
@@ -42,15 +42,15 @@ eos_barotr reader_eos_barotr_spline::load(const datasource g,
   lgspl_t seps        = g["eps_from_gm1"];
   lgspl_t shm1        = g["hm1_from_gm1"];
   lglgspl_t spress_si = g["press_from_gm1"];
-  lgspl_t scsnd_si    = g["csnd_from_gm1"];
+  lgspl_t scsnd_si    = g["csnd_from_rho"];
 
   opt_t stemp_mev     = g["temp_from_gm1"];
   opt_t sefrac        = g["efrac_from_gm1"];
 
   lglgspl_t sgm1   = sgm1_si.rescale_x(1./u.density());
   lglgspl_t srho   = srho_si / u.density();
   lglgspl_t spress = spress_si / u.pressure();
-  lgspl_t scsnd    = scsnd_si / u.velocity();
+  lgspl_t scsnd    = scsnd_si.rescale_x(1./u.density()) / u.velocity();
 
   return eos_barotr{ 
     std::make_shared<eos_barotr_spline>(sgm1, srho, seps, spress, 
@@ -71,7 +71,7 @@ void eos_barotr_spline::save(datasink g) const
   g["eps_from_gm1"] = eps_gm1;
   g["hm1_from_gm1"] = hm1_gm1;
   g["press_from_gm1"] = p_gm1 * u.pressure();
-  g["csnd_from_gm1"] = csnd_gm1 * u.velocity();
+  g["csnd_from_rho"] = csnd_rho.rescale_x(u.density()) * u.velocity();
 
 
   if (!zerotemp) 

library/EOS_Barotropic/eos_barotr_spline_impl.h

@@ -104,6 +104,14 @@ class eos_barotr_spline : public eos_barotr_impl {
     real_t gm1      ///< \f$ g-1 \f$
   ) const final;
 
+  ///Compute adiabatic soundspeed \f$ c_s \f$
+  /**Assumes input is in the valid range, no checks are performed.*/
+  real_t csnd_from_rho_gm1(
+    real_t rho,     ///<Rest mass density  \f$ \rho \f$
+    real_t gm1      ///< \f$ g-1 \f$
+  ) const final;
+
+
   ///Returns temperature \f$ T \f$
   /**Assumes input is in the valid range, no checks are performed.*/
   real_t temp(
@@ -124,9 +132,8 @@ class eos_barotr_spline : public eos_barotr_impl {
 
   private:
 
-  static auto get_rggm1(const lglgspl_t&, const lgspl_t&, 
-                        const lglgspl_t&, const lgspl_t&, 
-                        const lglgspl_t&, const lgspl_t&, 
+  static auto get_rggm1(const lgspl_t&, const lglgspl_t&, 
+                        const lgspl_t&, const lglgspl_t&, 
                         const opt_t&, const opt_t&)  
   -> range;
 
@@ -135,7 +142,7 @@ class eos_barotr_spline : public eos_barotr_impl {
   lglgspl_t p_gm1;
   lgspl_t hm1_gm1;
   lglgspl_t rho_gm1;
-  lgspl_t csnd_gm1;
+  lgspl_t csnd_rho;
   opt_t temp_gm1;
   opt_t efrac_gm1;
 

library/EOS_Barotropic/eos_barotropic.cc

@@ -71,13 +71,13 @@ auto eos_barotr::is_gm1_valid(real_t gm1) const -> bool
 auto eos_barotr::at_rho(real_t rho) const -> state
 {
   if (!is_rho_valid(rho)) return {};
-  return {impl(), impl().gm1_from_rho(rho)};
+  return {impl(), impl().gm1_from_rho(rho), rho};
 }
 
 auto eos_barotr::at_gm1(real_t gm1) const -> state
 {
   if (!is_gm1_valid(gm1)) return {};
-  return {impl(), gm1};
+  return {impl(), gm1, impl().rho(gm1)};
 }
 
 
@@ -99,14 +99,14 @@ auto eos_barotr::state::press() const -> real_t
 
 auto eos_barotr::state::rho() const -> real_t 
 {
-  real_t rho = impl().rho(gm1_);  
-  assert(rho >= 0);
-  return rho;
+  if (!am_ok()) throw(eos_barotr_invalid::invalid());
+  assert(rho_ >= 0);
+  return rho_;
 }
 
 auto eos_barotr::state::csnd() const -> real_t 
 {
-  real_t cs = impl().csnd(gm1_);  
+  real_t cs = impl().csnd_from_rho_gm1(rho_, gm1_);  
   assert(cs < 1.0);
   assert(cs >= 0);
   return cs;
@@ -238,6 +238,11 @@ auto eos_barotr::descr_str() const -> std::string
   return impl().descr_str();
 }
 
+auto eos_barotr_impl::csnd_from_rho_gm1(real_t rho, real_t gm1) const 
+-> real_t
+{
+  return csnd(gm1);
+}
 
 void eos_barotr_impl::save(datasink s) const
 {