Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <limits>
       8             : 
       9             : #include "DataStructures/Tensor/TypeAliases.hpp"
      10             : #include "Options/String.hpp"
      11             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      12             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      13             : #include "Utilities/GenerateInstantiations.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// \cond
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : class DataVector;
      22             : namespace domain {
      23             : namespace Tags {
      24             : template <size_t Dim, typename Frame>
      25             : struct Coordinates;
      26             : }  // namespace Tags
      27             : }  // namespace domain
      28             : namespace hydro {
      29             : namespace Tags {
      30             : template <typename DataType>
      31             : struct Pressure;
      32             : template <typename DataType>
      33             : struct RestMassDensity;
      34             : template <typename DataType>
      35             : struct SpecificInternalEnergy;
      36             : template <typename DataType>
      37             : struct Temperature;
      38             : template <typename DataType>
      39             : struct ElectronFraction;
      40             : }  // namespace Tags
      41             : }  // namespace hydro
      42             : /// \endcond
      43             : 
      44             : /// Contains all variable fixers.
      45             : namespace VariableFixing {
      46             : 
      47             : /// \ingroup VariableFixingGroup
      48             : /// \brief Applies a pressure and density floor dependent on the distance
      49             : /// to the origin.
      50             : ///
      51             : /// Applies the floors:
      52             : /// \f$\rho(r) \geq \rho_{\mathrm{fl}}(r) = C_\rho r^{k_\rho}\f$
      53             : /// and \f$P(r) \geq P_{\mathrm{fl}}(r) = C_p r^{k_p}\f$
      54             : /// when \f$ r > r_{min}\f$, where \f$C_\rho\f$ is given by the option
      55             : /// `ScaleDensityFloor`, \f$k_\rho\f$ is given by the option
      56             : /// `PowerDensityFloor`,  \f$C_p\f$ is given by the option
      57             : /// `ScalePressureFloor`, \f$k_p\f$ is given by the option
      58             : /// `PowerPressureFloor`, and \f$r_{min}\f$ is given by the option
      59             : /// `MinimumRadius`.
      60             : ///
      61             : /// \note In \cite Porth2016rfi, the following floors are applied:
      62             : /// \f$\rho(r) \geq \rho_{\mathrm{fl}}(r) = 10^{-5}r^{-3/2}\f$
      63             : /// and \f$P(r) \geq P_{\mathrm{fl}}(r) = \frac{1}{3} \times 10^{-7}r^{-5/2}\f$
      64             : ///
      65             : /// \note For 3d equations of state we treat the product \f$\rho T\f$ as the
      66             : /// "pressure" that we floor and then compute the pressure from the
      67             : /// temperature. The floor is applied after flooring the density.
      68             : template <size_t Dim>
      69           1 : class RadiallyFallingFloor {
      70             :  public:
      71             :   /// \brief The minimum radius at which to begin applying the floors on the
      72             :   /// density and pressure.
      73           1 :   struct MinimumRadius {
      74           0 :     static constexpr Options::String help =
      75             :         "The radius at which to begin applying the lower bound.";
      76           0 :     using type = double;
      77           0 :     static double lower_bound() { return 0.0; }
      78             :   };
      79             : 
      80             :   /// \brief The scale of the floor of the rest mass density.
      81           1 :   struct ScaleDensityFloor {
      82           0 :     static constexpr Options::String help =
      83             :         "The rest mass density floor at r = 1.";
      84           0 :     using type = double;
      85           0 :     static double lower_bound() { return 0.0; }
      86             :   };
      87             : 
      88             :   /// \brief The power of the radius of the floor of the rest mass density.
      89           1 :   struct PowerDensityFloor {
      90           0 :     static constexpr Options::String help =
      91             :         "Radial power for the floor of the rest mass density.";
      92           0 :     using type = double;
      93             :   };
      94             : 
      95             :   /// \brief The scale of the floor of the pressure.
      96           1 :   struct ScalePressureFloor {
      97           0 :     static constexpr Options::String help = "The pressure floor at r = 1.";
      98           0 :     using type = double;
      99           0 :     static double lower_bound() { return 0.0; }
     100             :   };
     101             : 
     102             :   /// \brief The power of the radius of the floor of the pressure.
     103           1 :   struct PowerPressureFloor {
     104           0 :     static constexpr Options::String help =
     105             :         "The radial power for the floor of the pressure.";
     106           0 :     using type = double;
     107             :   };
     108             : 
     109           0 :   using options =
     110             :       tmpl::list<MinimumRadius, ScaleDensityFloor, PowerDensityFloor,
     111             :                  ScalePressureFloor, PowerPressureFloor>;
     112           0 :   static constexpr Options::String help = {
     113             :       "Applies a pressure and density floor dependent on the distance to the "
     114             :       "origin."};
     115             : 
     116           0 :   RadiallyFallingFloor(double minimum_radius_at_which_to_apply_floor,
     117             :                        double rest_mass_density_scale,
     118             :                        double rest_mass_density_power, double pressure_scale,
     119             :                        double pressure_power);
     120             : 
     121           0 :   RadiallyFallingFloor() = default;
     122           0 :   RadiallyFallingFloor(const RadiallyFallingFloor& /*rhs*/) = default;
     123           0 :   RadiallyFallingFloor& operator=(const RadiallyFallingFloor& /*rhs*/) =
     124             :       default;
     125           0 :   RadiallyFallingFloor(RadiallyFallingFloor&& /*rhs*/) = default;
     126           0 :   RadiallyFallingFloor& operator=(RadiallyFallingFloor&& /*rhs*/) = default;
     127           0 :   ~RadiallyFallingFloor() = default;
     128             : 
     129             :   // NOLINTNEXTLINE(google-runtime-references)
     130           0 :   void pup(PUP::er& p);
     131             : 
     132           0 :   using return_tags =
     133             :       tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
     134             :                  hydro::Tags::Pressure<DataVector>,
     135             :                  hydro::Tags::SpecificInternalEnergy<DataVector>,
     136             :                  hydro::Tags::Temperature<DataVector>,
     137             :                  hydro::Tags::ElectronFraction<DataVector>>;
     138           0 :   using argument_tags =
     139             :       tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
     140             :                  hydro::Tags::GrmhdEquationOfState>;
     141             : 
     142             :   template <size_t ThermodynamicDim>
     143           0 :   void operator()(
     144             :       gsl::not_null<Scalar<DataVector>*> density,
     145             :       gsl::not_null<Scalar<DataVector>*> pressure,
     146             :       gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     147             :       gsl::not_null<Scalar<DataVector>*> temperature,
     148             :       gsl::not_null<Scalar<DataVector>*> electron_fraction,
     149             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
     150             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     151             :           equation_of_state) const;
     152             : 
     153             :  private:
     154             :   template <size_t LocalDim>
     155             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     156           0 :   friend bool operator==(const RadiallyFallingFloor<LocalDim>& lhs,
     157             :                          const RadiallyFallingFloor<LocalDim>& rhs);
     158             : 
     159           0 :   double minimum_radius_at_which_to_apply_floor_{
     160             :       std::numeric_limits<double>::signaling_NaN()};
     161           0 :   double rest_mass_density_scale_{std::numeric_limits<double>::signaling_NaN()};
     162           0 :   double rest_mass_density_power_{std::numeric_limits<double>::signaling_NaN()};
     163           0 :   double pressure_scale_{std::numeric_limits<double>::signaling_NaN()};
     164           0 :   double pressure_power_{std::numeric_limits<double>::signaling_NaN()};
     165             : };
     166             : 
     167             : template <size_t Dim>
     168           0 : bool operator!=(const RadiallyFallingFloor<Dim>& lhs,
     169             :                 const RadiallyFallingFloor<Dim>& rhs);
     170             : 
     171             : }  // namespace VariableFixing