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 <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             : #include <vector>
      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 <typename LowDensityEoS>
     103           1 : class Enthalpy : public EquationOfState<true, 1> {
     104             :  private:
     105           0 :   struct Coefficients {
     106           0 :     std::vector<double> polynomial_coefficients;
     107           0 :     std::vector<double> sin_coefficients;
     108           0 :     std::vector<double> 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<double> in_polynomial_coefficients,
     117             :                  std::vector<double> in_sin_coefficients,
     118             :                  std::vector<double> in_cos_coefficients, double in_trig_scale,
     119             :                  double in_reference_density,
     120             :                  double in_exponential_constant =
     121             :                      std::numeric_limits<double>::quiet_NaN());
     122             : 
     123           0 :     Enthalpy<LowDensityEoS>::Coefficients compute_exponential_integral(
     124             :         const std::pair<double, double>& initial_condition);
     125           0 :     Enthalpy<LowDensityEoS>::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<double>;
     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<double>;
     167           0 :     static constexpr Options::String help = {"Sine coefficients b_j"};
     168             :   };
     169           0 :   struct CosCoefficients {
     170           0 :     using type = std::vector<double>;
     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<LowDensityEoS>();
     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<ReferenceDensity, MaximumDensity, MinimumDensity, TrigScaling,
     203             :                  PolynomialCoefficients, SinCoefficients, CosCoefficients,
     204             :                  StitchedLowDensityEoS, TransitionDeltaEpsilon>;
     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<double>& polynomial_coefficients,
     216             :            const std::vector<double>& sin_coefficients,
     217             :            const std::vector<double>& 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<true, 1>), 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<double>::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<double>::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<double>::signaling_NaN();
     273           0 :   double minimum_density_ = std::numeric_limits<double>::signaling_NaN();
     274           0 :   double maximum_density_ = std::numeric_limits<double>::signaling_NaN();
     275           0 :   double minimum_enthalpy_ = std::numeric_limits<double>::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