SpECTRE Documentation Coverage Report
Current view: top level - PointwiseFunctions/Hydro/EquationsOfState - Tabulated3d.hpp Hit Total Coverage
Commit: 9a905b0737f373631c1b8e8389b8f26e67fa5313 Lines: 22 54 40.7 %
Date: 2024-03-28 09:03:18
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14