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/Context.hpp"
      11             : #include "Options/String.hpp"
      12             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      13             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      14             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      15             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : #include "Utilities/TMPL.hpp"
      18             : 
      19             : /// \cond
      20             : class DataVector;
      21             : 
      22             : namespace PUP {
      23             : class er;
      24             : }  // namespace PUP
      25             : /// \endcond
      26             : 
      27             : namespace VariableFixing {
      28             : 
      29             : /*!
      30             :  * \ingroup VariableFixingGroup
      31             :  * \brief Fix the primitive variables to an atmosphere in low density regions
      32             :  *
      33             :  * If the rest mass density is below \f$\rho_{\textrm{cutoff}}\f$
      34             :  * (DensityCutoff), it is set to \f$\rho_{\textrm{atm}}\f$
      35             :  * (DensityOfAtmosphere), and the pressure, and specific internal energy (for
      36             :  * one-dimensional equations of state) are adjusted to satisfy the equation of
      37             :  * state.  For a two-dimensional equation of state, the specific internal energy
      38             :  * is set to zero. In addition, the spatial velocity is set to zero, and the
      39             :  * Lorentz factor is set to one.
      40             :  *
      41             :  * If the rest mass density is above \f$\rho_{\textrm{cutoff}}\f$ but below
      42             :  * \f$\rho_{\textrm{transition}}\f$ (TransitionDensityCutoff) then the velocity
      43             :  * is rescaled such that
      44             :  *
      45             :  * \f{align*}{
      46             :  * \sqrt{v^i v_i}\le \frac{(\rho-\rho_{\textrm{cutoff}})}
      47             :  * {(\rho_{\textrm{transition}} - \rho_{\textrm{cutoff}})} v_{\max}
      48             :  * \f}
      49             :  *
      50             :  * where \f$v_{\max}\f$ (MaxVelocityMagnitude) is the maximum allowed magnitude
      51             :  * of the velocity. This prescription follows Appendix 2.d of
      52             :  * \cite Muhlberger2014pja Note that we require
      53             :  * \f$\rho_{\textrm{transition}}\in(\rho_{\textrm{cutoff}},
      54             :  *  10\rho_{\textrm{atm}}]\f$
      55             :  */
      56             : template <size_t Dim>
      57           1 : class FixToAtmosphere {
      58             :  public:
      59             :   /// \brief Rest mass density of the atmosphere
      60           1 :   struct DensityOfAtmosphere {
      61           0 :     using type = double;
      62           0 :     static type lower_bound() { return 0.0; }
      63           0 :     static constexpr Options::String help = {"Density of atmosphere"};
      64             :   };
      65             :   /// \brief Rest mass density at which to impose the atmosphere. Should be
      66             :   /// greater than or equal to the density of the atmosphere.
      67           1 :   struct DensityCutoff {
      68           0 :     using type = double;
      69           0 :     static type lower_bound() { return 0.0; }
      70           0 :     static constexpr Options::String help = {
      71             :         "Density to impose atmosphere at. Must be >= rho_atm"};
      72             :   };
      73             :   /// \brief For densities between DensityOfAtmosphere and
      74             :   /// TransitionDensityCutoff the velocity is transitioned away from atmosphere
      75             :   /// to avoid abrupt cutoffs.
      76             :   ///
      77             :   /// This value must not be larger than `10 * DensityOfAtmosphere`.
      78           1 :   struct TransitionDensityCutoff {
      79           0 :     using type = double;
      80           0 :     static type lower_bound() { return 0.0; }
      81           0 :     static constexpr Options::String help = {
      82             :         "For densities between DensityOfAtmosphere and TransitionDensityCutoff "
      83             :         "the velocity is transitioned away from atmosphere to avoid abrupt "
      84             :         "cutoffs.\n\n"
      85             :         "This value must not be larger than 10 * DensityOfAtmosphere."};
      86             :   };
      87             :   /// \brief The maximum magnitude of the velocity when the density is below
      88             :   /// `TransitionDensityCutoff`
      89           1 :   struct MaxVelocityMagnitude {
      90           0 :     using type = double;
      91           0 :     static type lower_bound() { return 0.0; }
      92           0 :     static type upper_bound() { return 1.0; }
      93           0 :     static constexpr Options::String help = {
      94             :         "The maximum sqrt(v^i v^j gamma_{ij}) allowed when the density is "
      95             :         "below TransitionDensityCutoff."};
      96             :   };
      97             : 
      98           0 :   using options =
      99             :       tmpl::list<DensityOfAtmosphere, DensityCutoff, TransitionDensityCutoff,
     100             :                  MaxVelocityMagnitude>;
     101           0 :   static constexpr Options::String help = {
     102             :       "If the rest mass density is below DensityCutoff, it is set\n"
     103             :       "to DensityOfAtmosphere, and the pressure, and specific internal energy\n"
     104             :       "(for one-dimensional equations of state) are\n"
     105             :       "adjusted to satisfy the equation of state. For a two-dimensional\n"
     106             :       "equation of state, the specific internal energy is set to zero.\n"
     107             :       "In addition, the spatial velocity is set to zero, and the Lorentz\n"
     108             :       "factor is set to one.\n"};
     109             : 
     110           0 :   FixToAtmosphere(double density_of_atmosphere, double density_cutoff,
     111             :                   double transition_density_cutoff,
     112             :                   double max_velocity_magnitude,
     113             :                   const Options::Context& context = {});
     114             : 
     115           0 :   FixToAtmosphere() = default;
     116           0 :   FixToAtmosphere(const FixToAtmosphere& /*rhs*/) = default;
     117           0 :   FixToAtmosphere& operator=(const FixToAtmosphere& /*rhs*/) = default;
     118           0 :   FixToAtmosphere(FixToAtmosphere&& /*rhs*/) = default;
     119           0 :   FixToAtmosphere& operator=(FixToAtmosphere&& /*rhs*/) = default;
     120           0 :   ~FixToAtmosphere() = default;
     121             : 
     122             :   // NOLINTNEXTLINE(google-runtime-references)
     123           0 :   void pup(PUP::er& p);
     124             : 
     125           0 :   using return_tags =
     126             :       tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
     127             :                  hydro::Tags::SpecificInternalEnergy<DataVector>,
     128             :                  hydro::Tags::SpatialVelocity<DataVector, Dim>,
     129             :                  hydro::Tags::LorentzFactor<DataVector>,
     130             :                  hydro::Tags::Pressure<DataVector>,
     131             :                  hydro::Tags::Temperature<DataVector>>;
     132           0 :   using argument_tags = tmpl::list<hydro::Tags::ElectronFraction<DataVector>,
     133             :                                    gr::Tags::SpatialMetric<DataVector, Dim>,
     134             :                                    hydro::Tags::EquationOfStateBase>;
     135             : 
     136             :   // for use in `db::mutate_apply`
     137             :   template <size_t ThermodynamicDim>
     138           0 :   void operator()(
     139             :       gsl::not_null<Scalar<DataVector>*> rest_mass_density,
     140             :       gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     141             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     142             :           spatial_velocity,
     143             :       gsl::not_null<Scalar<DataVector>*> lorentz_factor,
     144             :       gsl::not_null<Scalar<DataVector>*> pressure,
     145             :       gsl::not_null<Scalar<DataVector>*> temperature,
     146             :       const Scalar<DataVector>& electron_fraction,
     147             :       const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
     148             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     149             :           equation_of_state) const;
     150             : 
     151             :  private:
     152             :   template <size_t ThermodynamicDim>
     153           0 :   void set_density_to_atmosphere(
     154             :       gsl::not_null<Scalar<DataVector>*> rest_mass_density,
     155             :       gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     156             :       gsl::not_null<Scalar<DataVector>*> temperature,
     157             :       gsl::not_null<Scalar<DataVector>*> pressure,
     158             :       const Scalar<DataVector>& electron_fraction,
     159             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     160             :           equation_of_state,
     161             :       size_t grid_index) const;
     162             : 
     163           0 :   void set_to_magnetic_free_transition(
     164             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     165             :           spatial_velocity,
     166             :       gsl::not_null<Scalar<DataVector>*> lorentz_factor,
     167             :       const Scalar<DataVector>& rest_mass_density,
     168             :       const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
     169             :       size_t grid_index) const;
     170             : 
     171             :   template <size_t SpatialDim>
     172             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     173           0 :   friend bool operator==(const FixToAtmosphere<SpatialDim>& lhs,
     174             :                          const FixToAtmosphere<SpatialDim>& rhs);
     175             : 
     176           0 :   double density_of_atmosphere_{std::numeric_limits<double>::signaling_NaN()};
     177           0 :   double density_cutoff_{std::numeric_limits<double>::signaling_NaN()};
     178           0 :   double transition_density_cutoff_{
     179             :       std::numeric_limits<double>::signaling_NaN()};
     180           0 :   double max_velocity_magnitude_{std::numeric_limits<double>::signaling_NaN()};
     181             : };
     182             : 
     183             : template <size_t Dim>
     184           0 : bool operator!=(const FixToAtmosphere<Dim>& lhs,
     185             :                 const FixToAtmosphere<Dim>& rhs);
     186             : 
     187             : }  // namespace VariableFixing