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