SpECTRE Documentation Coverage Report
 Current view: top level - PointwiseFunctions/Hydro/EquationsOfState - Enthalpy.hpp Hit Total Coverage Commit: cd74d65bdc718fd7e344eaec61dc6334dd4d366b Lines: 6 86 7.0 % Date: 2022-08-12 23:56:47 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 7 : #include 8 : #include 9 : #include 10 : #include 11 : #include 12 : #include 13 : #include 14 : #include 15 : #include 16 : 17 : #include "DataStructures/Tensor/TypeAliases.hpp" 18 : #include "Options/Options.hpp" 19 : #include "Parallel/CharmPupable.hpp" 20 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" // IWYU pragma: keep 21 : #include "Utilities/Math.hpp" 22 : #include "Utilities/TMPL.hpp" 23 : 24 : /// \cond 25 : class DataVector; 26 : 27 : /// \endcond 28 : 29 : // IWYU pragma: no_forward_declare Tensor 30 : 31 : namespace EquationsOfState { 32 : 33 : /*! 34 : * \ingroup EquationsOfStateGroup 35 : * \brief An equation of state given by parametrized enthalpy 36 : * 37 : * This equation of state is determined as a function of \f$x = 38 : * \ln(\rho/\rho_0)\f$ where \f$\rho\f$ is the rest mass density and 39 : * \f$\rho_0\f$ is the provided reference density. 40 : * The pseudo-enthalpy \f$h \equiv (p + rho + u)/rho\f$ 41 : * is expanded as 42 : * 43 : * \f{equation} 44 : * h(x) = \sum_i a_i x^i + \sum_j b_j \sin(jkx) + c_j \cos(jkx) 45 : * \f} 46 : * 47 : * This form allows for convenient calculation of thermodynamic 48 : * quantities for a cold equation of state. For example 49 : * 50 : * \f{equation} 51 : * h(x) = \frac{d e} {d \rho} |_{x = \log(\rho/\rho_0)} 52 : * \f} 53 : * 54 : * where \f$e\f$ is the total energy density. At the same time \f$dx = 55 : * d\rho/\rho \f$ so \f$\rho_0 e^x dx = d \rho \f$ Therefore, 56 : * 57 : * \f{equation} 58 : * e(x) - e(x_0) = \int_{x_0}^x h(x') e^{x'} dx ' 59 : * \f} 60 : * 61 : * This can be computed analytically because 62 : * 63 : * \f{equation} 64 : * \int a_i \frac{x^i}{i!} e^{x} dx = \sum_{j \leq i} a_i (-1)^{i-j} 65 : * \frac{(x)^{j}}{j!} 66 : * + C \f} 67 : * 68 : * and 69 : * 70 : * \f{equation} 71 : * \int b_j \sin(j k x) e^x dx = b_j e^x \frac{\sin(jkx) - j k \cos(jkx)}{j^2 72 : * k^2 + 1} \f} 73 : * 74 : * \f{equation} 75 : * \int c_j \cos(j k x) e^x dx = b_j e^x \frac{\cos(jkx) + j k \sin(jkx)}{j^2 76 : * k^2 + 1} \f} 77 : * 78 : * From this most other thermodynamic quantities can be computed 79 : * analytically 80 : * 81 : * The internal energy density 82 : * \f{equation} 83 : * \epsilon(x)\rho(x) = e(x) - \rho(x) 84 : * \f} 85 : * 86 : * The pressure 87 : * \f{equation} 88 : * p(x) = \rho(x) h(x) - e(x) 89 : * \f} 90 : * 91 : * The derivative of the pressure with respect to the rest mass density 92 : * \f{equation} 93 : * \chi(x) = \frac{dp}{d\rho} |_{x = x(\rho)} = \frac{dh}{dx} 94 : * \f} 95 : * 96 : * Below the minimum density, a spectral parameterization 97 : * is used. 98 : * 99 : * 100 : * 101 : */ 102 : template 103 1 : class Enthalpy : public EquationOfState { 104 : private: 105 0 : struct Coefficients { 106 0 : std::vector polynomial_coefficients; 107 0 : std::vector sin_coefficients; 108 0 : std::vector cos_coefficients; 109 0 : double trig_scale; 110 0 : double reference_density; 111 0 : bool has_exponential_prefactor; 112 0 : double exponential_external_constant; 113 0 : Coefficients() = default; 114 0 : ~Coefficients() = default; 115 0 : Coefficients(const Coefficients& coefficients) = default; 116 0 : Coefficients(std::vector in_polynomial_coefficients, 117 : std::vector in_sin_coefficients, 118 : std::vector in_cos_coefficients, double in_trig_scale, 119 : double in_reference_density, 120 : double in_exponential_constant = 121 : std::numeric_limits::quiet_NaN()); 122 : 123 0 : Enthalpy::Coefficients compute_exponential_integral( 124 : const std::pair& initial_condition); 125 0 : Enthalpy::Coefficients compute_derivative(); 126 0 : void pup(PUP::er& p); 127 : }; 128 : 129 : public: 130 0 : static constexpr size_t thermodynamic_dim = 1; 131 0 : static constexpr bool is_relativistic = true; 132 : 133 0 : struct ReferenceDensity { 134 0 : using type = double; 135 0 : static constexpr Options::String help = {"Reference density rho_0"}; 136 0 : static double lower_bound() { return 0.0; } 137 : }; 138 : 139 0 : struct MinimumDensity { 140 0 : using type = double; 141 0 : static constexpr Options::String help = { 142 : "Minimum valid density rho_min," 143 : " for this parametrization"}; 144 0 : static double lower_bound() { return 0.0; } 145 : }; 146 0 : struct MaximumDensity { 147 0 : using type = double; 148 0 : static constexpr Options::String help = {"Maximum density for this EoS"}; 149 0 : static double lower_bound() { return 0.0; } 150 : }; 151 : 152 0 : struct PolynomialCoefficients { 153 0 : using type = std::vector; 154 0 : static constexpr Options::String help = {"Polynomial coefficients a_i"}; 155 : }; 156 : 157 0 : struct TrigScaling { 158 0 : using type = double; 159 0 : static constexpr Options::String help = { 160 : "Fundamental wavenumber of trig " 161 : "functions, k"}; 162 0 : static double lower_bound() { return 0.0; } 163 : }; 164 : 165 0 : struct SinCoefficients { 166 0 : using type = std::vector; 167 0 : static constexpr Options::String help = {"Sine coefficients b_j"}; 168 : }; 169 0 : struct CosCoefficients { 170 0 : using type = std::vector; 171 0 : static constexpr Options::String help = {"Cosine coefficients c_j"}; 172 : }; 173 0 : struct StitchedLowDensityEoS { 174 0 : using type = LowDensityEoS; 175 0 : static std::string name() { 176 : return pretty_type::short_name(); 177 : } 178 0 : static constexpr Options::String help = { 179 : "Low density EoS stitched at the MinimumDensity"}; 180 : }; 181 : 182 0 : struct TransitionDeltaEpsilon { 183 0 : using type = double; 184 0 : static constexpr Options::String help = { 185 : "the change in internal energy across the low-" 186 : "to-high-density transition, generically 0.0"}; 187 0 : static double lower_bound() { return 0.0; } 188 : }; 189 : 190 0 : static constexpr Options::String help = { 191 : "An EoS with a parametrized value h(log(rho/rho_0)) with h the specific " 192 : "enthalpy and rho the baryon rest mass density. The enthalpy is " 193 : "expanded as a sum of polynomial terms and trigonometric corrections. " 194 : "let x = log(rho/rho_0) in" 195 : "h(x) = \\sum_i a_ix^i + \\sum_j b_jsin(k * j * x) + c_jcos(k * j * x) " 196 : "Note that rho(x)(1+epsilon(x)) = int_0^x e^x' h((x') dx' can be " 197 : "computed " 198 : "analytically, and therefore so can " 199 : "P(x) = rho(x) * (h(x) - (1 + epsilon(x))) "}; 200 : 201 0 : using options = 202 : tmpl::list; 205 : 206 0 : Enthalpy() = default; 207 0 : Enthalpy(const Enthalpy&) = default; 208 0 : Enthalpy& operator=(const Enthalpy&) = default; 209 0 : Enthalpy(Enthalpy&&) = default; 210 0 : Enthalpy& operator=(Enthalpy&&) = default; 211 0 : ~Enthalpy() override = default; 212 : 213 0 : Enthalpy(double reference_density, double max_density, double min_density, 214 : double trig_scale, 215 : const std::vector& polynomial_coefficients, 216 : const std::vector& sin_coefficients, 217 : const std::vector& cos_coefficients, 218 : const LowDensityEoS& low_density_eos, 219 : const double transition_delta_epsilon); 220 : 221 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Enthalpy, 1) 222 : 223 0 : WRAPPED_PUPable_decl_base_template( // NOLINT 224 : SINGLE_ARG(EquationOfState), Enthalpy); 225 : 226 : /// The lower bound of the rest mass density that is valid for this EOS 227 1 : double rest_mass_density_lower_bound() const override { return 0.0; } 228 : 229 : /// The upper bound of the rest mass density that is valid for this EOS 230 1 : double rest_mass_density_upper_bound() const override { 231 : return std::numeric_limits::max(); 232 : } 233 : 234 : /// The lower bound of the specific internal energy that is valid for this EOS 235 : /// at the given rest mass density \f$\rho\f$ 236 1 : double specific_internal_energy_lower_bound( 237 : const double /* rest_mass_density */) const override { 238 : return 0.0; 239 : } 240 : 241 : /// The upper bound of the specific internal energy that is valid for this EOS 242 : /// at the given rest mass density \f$\rho\f$ 243 1 : double specific_internal_energy_upper_bound( 244 : const double /* rest_mass_density */) const override { 245 : return std::numeric_limits::max(); 246 : } 247 : 248 : /// The lower bound of the specific enthalpy that is valid for this EOS 249 1 : double specific_enthalpy_lower_bound() const override { return 1.0; } 250 : 251 : private: 252 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1) 253 : 254 : SPECTRE_ALWAYS_INLINE 255 0 : bool in_low_density_domain(const double density) const { 256 : return density < minimum_density_; 257 : } 258 : 259 0 : double x_from_density(const double density) const; 260 0 : double density_from_x(const double x) const; 261 0 : double energy_density_from_log_density(const double x) const; 262 0 : static double evaluate_coefficients( 263 : const Enthalpy::Coefficients& coefficients, const double x); 264 : 265 0 : double chi_from_density(const double density) const; 266 0 : double specific_internal_energy_from_density(const double density) const; 267 0 : double specific_enthalpy_from_density(const double density) const; 268 0 : double pressure_from_density(const double density) const; 269 0 : double pressure_from_log_density(const double x) const; 270 0 : double rest_mass_density_from_enthalpy(const double specific_enthalpy) const; 271 : 272 0 : double reference_density_ = std::numeric_limits::signaling_NaN(); 273 0 : double minimum_density_ = std::numeric_limits::signaling_NaN(); 274 0 : double maximum_density_ = std::numeric_limits::signaling_NaN(); 275 0 : double minimum_enthalpy_ = std::numeric_limits::signaling_NaN(); 276 : 277 0 : LowDensityEoS low_density_eos_; 278 0 : Coefficients coefficients_; 279 0 : Coefficients exponential_integral_coefficients_; 280 0 : Coefficients derivative_coefficients_; 281 : }; 282 : 283 : } // namespace EquationsOfState 

 Generated by: LCOV version 1.14