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