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 "Evolution/Systems/RelativisticEuler/Valencia/TagsDeclarations.hpp"
      11             : #include "Options/Context.hpp"
      12             : #include "Options/String.hpp"
      13             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// \cond
      18             : class DataVector;
      19             : 
      20             : namespace PUP {
      21             : class er;
      22             : }  // namespace PUP
      23             : /// \endcond
      24             : 
      25             : namespace RelativisticEuler {
      26             : namespace Valencia {
      27             : 
      28             : /*!
      29             :  * \brief Fix conservative variables using the method proposed in F. Foucart's
      30             :  * PhD thesis (Cornell)
      31             :  *
      32             :  * Fix the conservative variables at each grid point according to the following
      33             :  * recipe:
      34             :  *
      35             :  * - If \f$D < D_\text{cutoff}\f$, set \f$D = D_\text{min}\f$
      36             :  * - If \f$\tau < 0\f$, set \f$\tau = 0\f$
      37             :  * - If \f$\tilde S^2 > \tilde S_\text{max}^2 = \tilde\tau(\tilde\tau +
      38             :  *   2\tilde D)\f$, rescale \f$\tilde S_i\f$ by the factor
      39             :  *
      40             :  * \f{align*}
      41             :  * f = \sqrt{\frac{(1 - \epsilon)\tilde S_\text{max}^2}{\tilde S^2 +
      42             :  * 10^{-16}\tilde D^2}} < 1 \f}
      43             :  *
      44             :  * (if \f$ f\geq 1\f$, do not fix \f$\tilde S_i\f$),
      45             :  * where \f$D_\text{min}\f$, \f$D_\text{cutoff}\f$ and \f$\epsilon\f$
      46             :  * are small and positive input parameters, typically chosen in consistency with
      47             :  * some variable fixing applied to the primitives variables
      48             :  * (e.g. VariableFixing::FixToAtmosphere). Note that, in order to fix
      49             :  * \f$\tilde S_i\f$, \f$\epsilon\f$ is chosen so that \f$1 - \epsilon\f$
      50             :  * (the relevant quantity used in the equations)
      51             :  * is sufficiently close (but smaller) than 1.
      52             :  */
      53             : template <size_t Dim>
      54           1 : class FixConservatives {
      55             :  public:
      56             :   /// The minimum value of the rest mass density times the Lorentz factor,
      57             :   /// \f$D\f$
      58           1 :   struct MinimumValueOfD {
      59           0 :     using type = double;
      60           0 :     static constexpr Options::String help = {
      61             :         "Minimum value of rest-mass density times Lorentz factor"};
      62           0 :     static type lower_bound() { return 0.0; }
      63             :   };
      64             : 
      65             :   /// The cutoff below which \f$D\f$ is set to `MinimumValueOfD`
      66           1 :   struct CutoffD {
      67           0 :     using type = double;
      68           0 :     static constexpr Options::String help = {
      69             :         "Cutoff below which D is set to MinimumValueOfD"};
      70           0 :     static type lower_bound() { return 0.0; }
      71             :   };
      72             : 
      73             :   /// The safety factor to fix \f$\tilde S_i\f$
      74           1 :   struct SafetyFactorForS {
      75           0 :     using type = double;
      76           0 :     static constexpr Options::String help = {
      77             :         "Safety factor for momentum density bound."};
      78           0 :     static type lower_bound() { return std::numeric_limits<double>::epsilon(); }
      79           0 :     static type upper_bound() { return 1.0; }
      80             :   };
      81             : 
      82           0 :   using options = tmpl::list<MinimumValueOfD, CutoffD, SafetyFactorForS>;
      83           0 :   static constexpr Options::String help = {
      84             :       "Variable fixing used in Foucart's thesis."};
      85             : 
      86           0 :   FixConservatives() = default;
      87           0 :   FixConservatives(const FixConservatives& /*rhs*/) = default;
      88           0 :   FixConservatives& operator=(const FixConservatives& /*rhs*/) = default;
      89           0 :   FixConservatives(FixConservatives&& /*rhs*/) = default;
      90           0 :   FixConservatives& operator=(FixConservatives&& /*rhs*/) = default;
      91           0 :   ~FixConservatives() = default;
      92             : 
      93           0 :   FixConservatives(double minimum_rest_mass_density_times_lorentz_factor,
      94             :                    double rest_mass_density_times_lorentz_factor_cutoff,
      95             :                    double safety_factor_for_momentum_density,
      96             :                    const Options::Context& context = {});
      97             : 
      98             :   // NOLINTNEXTLINE(google-runtime-references)
      99           0 :   void pup(PUP::er& p);
     100             : 
     101           0 :   using return_tags =
     102             :       tmpl::list<RelativisticEuler::Valencia::Tags::TildeD,
     103             :                  RelativisticEuler::Valencia::Tags::TildeTau,
     104             :                  RelativisticEuler::Valencia::Tags::TildeS<Dim>>;
     105             : 
     106           0 :   using argument_tags =
     107             :       tmpl::list<gr::Tags::InverseSpatialMetric<DataVector, Dim>,
     108             :                  gr::Tags::SqrtDetSpatialMetric<DataVector>>;
     109             : 
     110           0 :   void operator()(
     111             :       gsl::not_null<Scalar<DataVector>*> tilde_d,
     112             :       gsl::not_null<Scalar<DataVector>*> tilde_tau,
     113             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> tilde_s,
     114             :       const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
     115             :       const Scalar<DataVector>& sqrt_det_spatial_metric) const;
     116             : 
     117             :  private:
     118             :   template <size_t SpatialDim>
     119           0 :   friend bool operator==(  // NOLINT(readability-redundant-declaration)
     120             :       const FixConservatives<SpatialDim>& lhs,
     121             :       const FixConservatives<SpatialDim>& rhs);
     122             : 
     123           0 :   double minimum_rest_mass_density_times_lorentz_factor_ =
     124             :       std::numeric_limits<double>::signaling_NaN();
     125           0 :   double rest_mass_density_times_lorentz_factor_cutoff_ =
     126             :       std::numeric_limits<double>::signaling_NaN();
     127           0 :   double one_minus_safety_factor_for_momentum_density_ =
     128             :       std::numeric_limits<double>::signaling_NaN();
     129             : };
     130             : 
     131             : template <size_t Dim>
     132           0 : bool operator!=(const FixConservatives<Dim>& lhs,
     133             :                 const FixConservatives<Dim>& rhs);
     134             : }  // namespace Valencia
     135             : }  // namespace RelativisticEuler