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 {
      73             :     Epsilon = 0,
      74             :     Pressure,
      75             :     CsSquared,
      76             :     DeltaMu,
      77             :     SpecificEntropy,
      78             :     NumberOfVars
      79             :   };
      80             : 
      81           0 :   Tabulated3D() = default;
      82           0 :   Tabulated3D(const Tabulated3D&) = default;
      83           0 :   Tabulated3D& operator=(const Tabulated3D&) = default;
      84           0 :   Tabulated3D(Tabulated3D&&) = default;
      85           0 :   Tabulated3D& operator=(Tabulated3D&&) = default;
      86           0 :   ~Tabulated3D() override = default;
      87             : 
      88           0 :   explicit Tabulated3D(const std::string& filename,
      89             :                        const std::string& subfilename);
      90             : 
      91           0 :   explicit Tabulated3D(std::vector<double> electron_fraction,
      92             :                        std::vector<double> log_density,
      93             :                        std::vector<double> log_temperature,
      94             :                        std::vector<double> table_data, double energy_shift,
      95             :                        double enthalpy_minimum);
      96             : 
      97           0 :   explicit Tabulated3D(const h5::EosTable& spectre_eos);
      98             : 
      99             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Tabulated3D, 3)
     100             : 
     101             :   template <class DataType>
     102           0 :   void convert_to_table_quantities(
     103             :       const gsl::not_null<Scalar<DataType>*> converted_electron_fraction,
     104             :       const gsl::not_null<Scalar<DataType>*> log_rest_mass_density,
     105             :       const gsl::not_null<Scalar<DataType>*> log_temperature,
     106             :       const Scalar<DataType>& electron_fraction,
     107             :       const Scalar<DataType>& rest_mass_density,
     108             :       const Scalar<DataType>& temperature) const {
     109             :     get(*converted_electron_fraction) = get(electron_fraction);
     110             :     get(*log_rest_mass_density) = get(rest_mass_density);
     111             :     get(*log_temperature) = get(temperature);
     112             : 
     113             :     // Enforce physicality of input
     114             :     // We reuse the same variables here, these are not log yet.
     115             :     enforce_physicality(*converted_electron_fraction, *log_rest_mass_density,
     116             :                         *log_temperature);
     117             : 
     118             :     // Table uses log T and log rho
     119             :     get(*log_rest_mass_density) = log(get(*log_rest_mass_density));
     120             :     get(*log_temperature) = log(get(*log_temperature));
     121             :   }
     122             : 
     123           0 :   std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
     124             :       const override;
     125             : 
     126           0 :   void initialize(std::vector<double> electron_fraction,
     127             :                   std::vector<double> log_density,
     128             :                   std::vector<double> log_temperature,
     129             :                   std::vector<double> table_data, double energy_shift,
     130             :                   double enthalpy_minimum);
     131             : 
     132             : 
     133           0 :   void initialize(const h5::EosTable& spectre_eos);
     134             : 
     135           0 :   bool is_equal(const EquationOfState<IsRelativistic, 3>& rhs) const override;
     136             : 
     137             :   /// \brief Returns `true` if the EOS is barotropic
     138           1 :   bool is_barotropic() const override { return false; }
     139             : 
     140             :   /// \brief Returns `true` if the EOS is in beta-equilibrium
     141           1 :   bool is_equilibrium() const override { return false; }
     142             : 
     143           0 :   bool operator==(const Tabulated3D<IsRelativistic>& rhs) const;
     144             : 
     145           0 :   bool operator!=(const Tabulated3D<IsRelativistic>& rhs) const;
     146             : 
     147             :   template <class DataType>
     148           0 :   Scalar<DataType> equilibrium_electron_fraction_from_density_temperature_impl(
     149             :       const Scalar<DataType>& rest_mass_density,
     150             :       const Scalar<DataType>& temperature) const;
     151             : 
     152             :   /// @{
     153             :   /*!
     154             :    * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
     155             :    * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
     156             :    */
     157           1 :   Scalar<double> equilibrium_electron_fraction_from_density_temperature(
     158             :       const Scalar<double>& rest_mass_density,
     159             :       const Scalar<double>& temperature) const override {
     160             :     return equilibrium_electron_fraction_from_density_temperature_impl<double>(
     161             :         rest_mass_density, temperature);
     162             :   }
     163             : 
     164           1 :   Scalar<DataVector> equilibrium_electron_fraction_from_density_temperature(
     165             :       const Scalar<DataVector>& rest_mass_density,
     166             :       const Scalar<DataVector>& temperature) const override {
     167             :     return equilibrium_electron_fraction_from_density_temperature_impl<
     168             :         DataVector>(rest_mass_density, temperature);
     169             :   }
     170             :   /// @}
     171             :   //
     172             : 
     173             :   template <typename DataType>
     174           0 :   void enforce_physicality(Scalar<DataType>& electron_fraction,
     175             :                            Scalar<DataType>& density,
     176             :                            Scalar<DataType>& temperature) const;
     177             : 
     178           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     179             :       SINGLE_ARG(EquationOfState<IsRelativistic, 3>), Tabulated3D);
     180             : 
     181             :   /// The lower bound of the electron fraction that is valid for this EOS
     182           1 :   double electron_fraction_lower_bound() const override {
     183             :     return table_electron_fraction_.front();
     184             :   }
     185             : 
     186             :   /// The upper bound of the electron fraction that is valid for this EOS
     187           1 :   double electron_fraction_upper_bound() const override {
     188             :     return table_electron_fraction_.back();
     189             :   }
     190             : 
     191             :   /// The lower bound of the rest mass density that is valid for this EOS
     192           1 :   double rest_mass_density_lower_bound() const override {
     193             :     return std::exp((table_log_density_.front()));
     194             :   }
     195             : 
     196             :   /// The upper bound of the rest mass density that is valid for this EOS
     197           1 :   double rest_mass_density_upper_bound() const override {
     198             :     return std::exp((table_log_density_.back()));
     199             :   }
     200             : 
     201             :   /// The lower bound of the temperature that is valid for this EOS
     202           1 :   double temperature_lower_bound() const override {
     203             :     return std::exp((table_log_temperature_.front()));
     204             :   }
     205             : 
     206             :   /// The upper bound of the temperature that is valid for this EOS
     207           1 :   double temperature_upper_bound() const override {
     208             :     return std::exp((table_log_temperature_.back()));
     209             :   }
     210             : 
     211             :   /// The lower bound of the specific internal energy that is valid for this EOS
     212             :   /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$
     213           1 :   double specific_internal_energy_lower_bound(
     214             :       const double rest_mass_density,
     215             :       const double electron_fraction) const override;
     216             : 
     217             :   /// The upper bound of the specific internal energy that is valid for this EOS
     218             :   /// at the given rest mass density \f$\rho\f$
     219           1 :   double specific_internal_energy_upper_bound(
     220             :       const double rest_mass_density,
     221             :       const double electron_fraction) const override;
     222             : 
     223             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     224           1 :   double specific_enthalpy_lower_bound() const override {
     225             :     return enthalpy_minimum_;
     226             :   }
     227             : 
     228             :   /// The baryon mass for this EoS
     229           1 :   double baryon_mass() const override {
     230             :     return hydro::units::geometric::neutron_mass;
     231             :   }
     232             : 
     233             :  private:
     234             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(3)
     235             : 
     236           0 :   void initialize_interpolator();
     237             : 
     238             :   /// Energy shift used to account for negative specific internal energies,
     239             :   /// which are only stored logarithmically
     240           1 :   double energy_shift_ = 0.;
     241             : 
     242             :   /// Enthalpy minium  across the table
     243           1 :   double enthalpy_minimum_ = 1.;
     244             : 
     245             :   /// Main interpolator for the EoS.
     246             :   /// The ordering is  \f$(\log T. \log \rho, Y_e)\f$.
     247             :   /// Assumed to be sorted in ascending order.
     248           1 :   intrp::UniformMultiLinearSpanInterpolation<3, NumberOfVars> interpolator_{};
     249             :   /// Electron fraction
     250           1 :   std::vector<double> table_electron_fraction_{};
     251             :   /// Logarithmic rest-mass denisty
     252           1 :   std::vector<double> table_log_density_{};
     253             :   /// Logarithmic temperature
     254           1 :   std::vector<double> table_log_temperature_{};
     255             :   /// Tabulate data. Entries are stated in the enum
     256           1 :   std::vector<double> table_data_{};
     257             : 
     258             :   /// Tolerance on upper bound for root finding
     259           1 :   static constexpr double upper_bound_tolerance_ = 0.9999;
     260             : };
     261             : 
     262             : /// \cond
     263             : template <bool IsRelativistic>
     264             : PUP::able::PUP_ID EquationsOfState::Tabulated3D<IsRelativistic>::my_PUP_ID = 0;
     265             : /// \endcond
     266             : 
     267             : }  // namespace EquationsOfState