Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <boost/preprocessor/arithmetic/dec.hpp>
       7             : #include <boost/preprocessor/arithmetic/inc.hpp>
       8             : #include <boost/preprocessor/control/expr_iif.hpp>
       9             : #include <boost/preprocessor/list/adt.hpp>
      10             : #include <boost/preprocessor/repetition/for.hpp>
      11             : #include <boost/preprocessor/repetition/repeat.hpp>
      12             : #include <boost/preprocessor/tuple/to_list.hpp>
      13             : #include <limits>
      14             : #include <pup.h>
      15             : 
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "DataStructures/Tensor/TypeAliases.hpp"
      18             : #include "IO/H5/EosTable.hpp"
      19             : #include "IO/H5/File.hpp"
      20             : #include "NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp"
      21             : #include "Options/String.hpp"
      22             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      23             : #include "PointwiseFunctions/Hydro/Units.hpp"
      24             : #include "Utilities/Serialization/CharmPupable.hpp"
      25             : #include "Utilities/TMPL.hpp"
      26             : 
      27             : /// \cond
      28             : class DataVector;
      29             : /// \endcond
      30             : 
      31             : namespace EquationsOfState {
      32             : /*!
      33             :  * \ingroup EquationsOfStateGroup
      34             :  * \brief Nuclear matter equation of state in tabulated form.
      35             :  *
      36             :  * The equation of state takes the form
      37             :  *
      38             :  * \f[
      39             :  * p = p (T, rho, Y_e)
      40             :  * \f]
      41             :  *
      42             :  * where \f$\rho\f$ is the rest mass density, \f$T\f$ is the
      43             :  * temperature, and \f$Y_e\f$ is the electron fraction.
      44             :  * The temperature is given in units of MeV.
      45             :  */
      46             : template <bool IsRelativistic>
      47           1 : class Tabulated3D : public EquationOfState<IsRelativistic, 3> {
      48             :  public:
      49           0 :   static constexpr size_t thermodynamic_dim = 3;
      50           0 :   static constexpr bool is_relativistic = IsRelativistic;
      51             : 
      52           0 :   static constexpr Options::String help = {
      53             :       "A tabulated three-dimensional equation of state.\n"
      54             :       "The energy density, pressure and sound speed "
      55             :       "are tabulated as a function of density, electron_fraction and "
      56             :       "temperature."};
      57             : 
      58           0 :   struct TableFilename {
      59           0 :     using type = std::string;
      60           0 :     static constexpr Options::String help{"File name of the EOS table"};
      61             :   };
      62             : 
      63           0 :   struct TableSubFilename {
      64           0 :     using type = std::string;
      65           0 :     static constexpr Options::String help{
      66             :         "Subfile name of the EOS table, e.g., 'dd2'."};
      67             :   };
      68             : 
      69           0 :   using options = tmpl::list<TableFilename, TableSubFilename>;
      70             : 
      71             :   /// Fields stored in the table
      72           0 :   enum : size_t { Epsilon = 0, Pressure, CsSquared, DeltaMu, NumberOfVars };
      73             : 
      74           0 :   Tabulated3D() = default;
      75           0 :   Tabulated3D(const Tabulated3D&) = default;
      76           0 :   Tabulated3D& operator=(const Tabulated3D&) = default;
      77           0 :   Tabulated3D(Tabulated3D&&) = default;
      78           0 :   Tabulated3D& operator=(Tabulated3D&&) = default;
      79           0 :   ~Tabulated3D() override = default;
      80             : 
      81           0 :   explicit Tabulated3D(const std::string& filename,
      82             :                        const std::string& subfilename);
      83             : 
      84           0 :   explicit Tabulated3D(std::vector<double> electron_fraction,
      85             :                        std::vector<double> log_density,
      86             :                        std::vector<double> log_temperature,
      87             :                        std::vector<double> table_data, double energy_shift,
      88             :                        double enthalpy_minimum);
      89             : 
      90           0 :   explicit Tabulated3D(const h5::EosTable& spectre_eos);
      91             : 
      92             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Tabulated3D, 3)
      93             : 
      94             :   template <class DataType>
      95           0 :   void convert_to_table_quantities(
      96             :       const gsl::not_null<Scalar<DataType>*> converted_electron_fraction,
      97             :       const gsl::not_null<Scalar<DataType>*> log_rest_mass_density,
      98             :       const gsl::not_null<Scalar<DataType>*> log_temperature,
      99             :       const Scalar<DataType>& electron_fraction,
     100             :       const Scalar<DataType>& rest_mass_density,
     101             :       const Scalar<DataType>& temperature) const {
     102             :     get(*converted_electron_fraction) = get(electron_fraction);
     103             :     get(*log_rest_mass_density) = get(rest_mass_density);
     104             :     get(*log_temperature) = get(temperature);
     105             : 
     106             :     // Enforce physicality of input
     107             :     // We reuse the same variables here, these are not log yet.
     108             :     enforce_physicality(*converted_electron_fraction, *log_rest_mass_density,
     109             :                         *log_temperature);
     110             : 
     111             :     // Table uses log T and log rho
     112             :     get(*log_rest_mass_density) = log(get(*log_rest_mass_density));
     113             :     get(*log_temperature) = log(get(*log_temperature));
     114             :   }
     115             : 
     116           0 :   std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
     117             :       const override;
     118             : 
     119           0 :   void initialize(std::vector<double> electron_fraction,
     120             :                   std::vector<double> log_density,
     121             :                   std::vector<double> log_temperature,
     122             :                   std::vector<double> table_data, double energy_shift,
     123             :                   double enthalpy_minimum);
     124             : 
     125             : 
     126           0 :   void initialize(const h5::EosTable& spectre_eos);
     127             : 
     128           0 :   bool is_equal(const EquationOfState<IsRelativistic, 3>& rhs) const override;
     129             : 
     130             :   /// \brief Returns `true` if the EOS is barotropic
     131           1 :   bool is_barotropic() const override { return false; }
     132             : 
     133             :   /// \brief Returns `true` if the EOS is in beta-equilibrium
     134           1 :   bool is_equilibrium() const override { return false; }
     135             : 
     136           0 :   bool operator==(const Tabulated3D<IsRelativistic>& rhs) const;
     137             : 
     138           0 :   bool operator!=(const Tabulated3D<IsRelativistic>& rhs) const;
     139             : 
     140             :   template <class DataType>
     141           0 :   Scalar<DataType> equilibrium_electron_fraction_from_density_temperature_impl(
     142             :       const Scalar<DataType>& rest_mass_density,
     143             :       const Scalar<DataType>& temperature) const;
     144             : 
     145             :   /// @{
     146             :   /*!
     147             :    * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
     148             :    * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
     149             :    */
     150           1 :   Scalar<double> equilibrium_electron_fraction_from_density_temperature(
     151             :       const Scalar<double>& rest_mass_density,
     152             :       const Scalar<double>& temperature) const override {
     153             :     return equilibrium_electron_fraction_from_density_temperature_impl<double>(
     154             :         rest_mass_density, temperature);
     155             :   }
     156             : 
     157           1 :   Scalar<DataVector> equilibrium_electron_fraction_from_density_temperature(
     158             :       const Scalar<DataVector>& rest_mass_density,
     159             :       const Scalar<DataVector>& temperature) const override {
     160             :     return equilibrium_electron_fraction_from_density_temperature_impl<
     161             :         DataVector>(rest_mass_density, temperature);
     162             :   }
     163             :   /// @}
     164             :   //
     165             : 
     166             :   template <typename DataType>
     167           0 :   void enforce_physicality(Scalar<DataType>& electron_fraction,
     168             :                            Scalar<DataType>& density,
     169             :                            Scalar<DataType>& temperature) const;
     170             : 
     171           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     172             :       SINGLE_ARG(EquationOfState<IsRelativistic, 3>), Tabulated3D);
     173             : 
     174             :   /// The lower bound of the electron fraction that is valid for this EOS
     175           1 :   double electron_fraction_lower_bound() const override {
     176             :     return table_electron_fraction_.front();
     177             :   }
     178             : 
     179             :   /// The upper bound of the electron fraction that is valid for this EOS
     180           1 :   double electron_fraction_upper_bound() const override {
     181             :     return table_electron_fraction_.back();
     182             :   }
     183             : 
     184             :   /// The lower bound of the rest mass density that is valid for this EOS
     185           1 :   double rest_mass_density_lower_bound() const override {
     186             :     return std::exp((table_log_density_.front()));
     187             :   }
     188             : 
     189             :   /// The upper bound of the rest mass density that is valid for this EOS
     190           1 :   double rest_mass_density_upper_bound() const override {
     191             :     return std::exp((table_log_density_.back()));
     192             :   }
     193             : 
     194             :   /// The lower bound of the temperature that is valid for this EOS
     195           1 :   double temperature_lower_bound() const override {
     196             :     return std::exp((table_log_temperature_.front()));
     197             :   }
     198             : 
     199             :   /// The upper bound of the temperature that is valid for this EOS
     200           1 :   double temperature_upper_bound() const override {
     201             :     return std::exp((table_log_temperature_.back()));
     202             :   }
     203             : 
     204             :   /// The lower bound of the specific internal energy that is valid for this EOS
     205             :   /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$
     206           1 :   double specific_internal_energy_lower_bound(
     207             :       const double rest_mass_density,
     208             :       const double electron_fraction) const override;
     209             : 
     210             :   /// The upper bound of the specific internal energy that is valid for this EOS
     211             :   /// at the given rest mass density \f$\rho\f$
     212           1 :   double specific_internal_energy_upper_bound(
     213             :       const double rest_mass_density,
     214             :       const double electron_fraction) const override;
     215             : 
     216             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     217           1 :   double specific_enthalpy_lower_bound() const override {
     218             :     return enthalpy_minimum_;
     219             :   }
     220             : 
     221             :   /// The baryon mass for this EoS
     222           1 :   double baryon_mass() const override {
     223             :     return hydro::units::geometric::neutron_mass;
     224             :   }
     225             : 
     226             :  private:
     227             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(3)
     228             : 
     229           0 :   void initialize_interpolator();
     230             : 
     231             :   /// Energy shift used to account for negative specific internal energies,
     232             :   /// which are only stored logarithmically
     233           1 :   double energy_shift_ = 0.;
     234             : 
     235             :   /// Enthalpy minium  across the table
     236           1 :   double enthalpy_minimum_ = 1.;
     237             : 
     238             :   /// Main interpolator for the EoS.
     239             :   /// The ordering is  \f$(\log T. \log \rho, Y_e)\f$.
     240             :   /// Assumed to be sorted in ascending order.
     241           1 :   intrp::UniformMultiLinearSpanInterpolation<3, NumberOfVars> interpolator_{};
     242             :   /// Electron fraction
     243           1 :   std::vector<double> table_electron_fraction_{};
     244             :   /// Logarithmic rest-mass denisty
     245           1 :   std::vector<double> table_log_density_{};
     246             :   /// Logarithmic temperature
     247           1 :   std::vector<double> table_log_temperature_{};
     248             :   /// Tabulate data. Entries are stated in the enum
     249           1 :   std::vector<double> table_data_{};
     250             : 
     251             :   /// Tolerance on upper bound for root finding
     252           1 :   static constexpr double upper_bound_tolerance_ = 0.9999;
     253             : };
     254             : 
     255             : /// \cond
     256             : template <bool IsRelativistic>
     257             : PUP::able::PUP_ID EquationsOfState::Tabulated3D<IsRelativistic>::my_PUP_ID = 0;
     258             : /// \endcond
     259             : 
     260             : }  // namespace EquationsOfState