FixConservatives.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <limits>
8 
10 #include "Evolution/Systems/RelativisticEuler/Valencia/TagsDeclarations.hpp"
11 #include "Options/Options.hpp"
12 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
13 #include "Utilities/Gsl.hpp"
14 #include "Utilities/TMPL.hpp"
15 
16 /// \cond
17 class DataVector;
18 
19 namespace PUP {
20 class er;
21 } // namespace PUP
22 /// \endcond
23 
24 namespace RelativisticEuler {
25 namespace Valencia {
26 
27 /*!
28  * \brief Fix conservative variables using the method proposed in F. Foucart's
29  * PhD thesis (Cornell)
30  *
31  * Fix the conservative variables at each grid point according to the following
32  * recipe:
33  *
34  * - If \f$D < D_\text{cutoff}\f$, set \f$D = D_\text{min}\f$
35  * - If \f$\tau < 0\f$, set \f$\tau = 0\f$
36  * - If \f$\tilde S^2 > \tilde S_\text{max}^2 = \tilde\tau(\tilde\tau +
37  * 2\tilde D)\f$, rescale \f$\tilde S_i\f$ by the factor
38  *
39  * \f{align*}
40  * f = \sqrt{\frac{(1 - \epsilon)\tilde S_\text{max}^2}{\tilde S^2 +
41  * 10^{-16}\tilde D^2}} < 1 \f}
42  *
43  * (if \f$ f\geq 1\f$, do not fix \f$\tilde S_i\f$),
44  * where \f$D_\text{min}\f$, \f$D_\text{cutoff}\f$ and \f$\epsilon\f$
45  * are small and positive input parameters, typically chosen in consistency with
46  * some variable fixing applied to the primitives variables
47  * (e.g. VariableFixing::FixToAtmosphere). Note that, in order to fix
48  * \f$\tilde S_i\f$, \f$\epsilon\f$ is chosen so that \f$1 - \epsilon\f$
49  * (the relevant quantity used in the equations)
50  * is sufficiently close (but smaller) than 1.
51  */
52 template <size_t Dim>
54  public:
55  /// The minimum value of the rest mass density times the Lorentz factor,
56  /// \f$D\f$
57  struct MinimumValueOfD {
58  using type = double;
59  static constexpr OptionString help = {
60  "Minimum value of rest-mass density times Lorentz factor"};
61  static type lower_bound() noexcept { return 0.0; }
62  };
63 
64  /// The cutoff below which \f$D\f$ is set to `MinimumValueOfD`
65  struct CutoffD {
66  using type = double;
67  static constexpr OptionString help = {
68  "Cutoff below which D is set to MinimumValueOfD"};
69  static type lower_bound() noexcept { return 0.0; }
70  };
71 
72  /// The safety factor to fix \f$\tilde S_i\f$
74  using type = double;
75  static constexpr OptionString help = {
76  "Safety factor for momentum density bound."};
77  static type lower_bound() noexcept {
79  }
80  static type upper_bound() noexcept { return 1.0; }
81  };
82 
83  using options = tmpl::list<MinimumValueOfD, CutoffD, SafetyFactorForS>;
84  static constexpr OptionString help = {
85  "Variable fixing used in Foucart's thesis."};
86 
87  FixConservatives() = default;
88  FixConservatives(const FixConservatives& /*rhs*/) = default;
89  FixConservatives& operator=(const FixConservatives& /*rhs*/) = default;
90  FixConservatives(FixConservatives&& /*rhs*/) noexcept = default;
91  FixConservatives& operator=(FixConservatives&& /*rhs*/) noexcept = default;
92  ~FixConservatives() = default;
93 
94  FixConservatives(double minimum_rest_mass_density_times_lorentz_factor,
95  double rest_mass_density_times_lorentz_factor_cutoff,
96  double safety_factor_for_momentum_density,
97  const OptionContext& context = {});
98 
99  // NOLINTNEXTLINE(google-runtime-references)
100  void pup(PUP::er& p) noexcept;
101 
102  using return_tags =
106 
107  using argument_tags = tmpl::list<gr::Tags::InverseSpatialMetric<Dim>,
109 
110  void operator()(
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 noexcept;
116 
117  private:
118  template <size_t SpatialDim>
119  friend bool operator==( // NOLINT(readability-redundant-declaration)
120  const FixConservatives<SpatialDim>& lhs,
121  const FixConservatives<SpatialDim>& rhs) noexcept;
122 
123  double minimum_rest_mass_density_times_lorentz_factor_ =
125  double rest_mass_density_times_lorentz_factor_cutoff_ =
127  double one_minus_safety_factor_for_momentum_density_ =
129 };
130 
131 template <size_t Dim>
132 bool operator!=(const FixConservatives<Dim>& lhs,
133  const FixConservatives<Dim>& rhs) noexcept;
134 } // namespace Valencia
135 } // namespace RelativisticEuler
The densitized momentum density .
Definition: Tags.hpp:42
The safety factor to fix .
Definition: FixConservatives.hpp:73
Definition: Strahlkorper.hpp:14
T epsilon(T... args)
T signaling_NaN(T... args)
Items related to evolving the relativistic Euler system.
Definition: EvolveValenciaFwd.hpp:8
Definition: Tags.hpp:49
Fix conservative variables using the method proposed in F. Foucart&#39;s PhD thesis (Cornell) ...
Definition: FixConservatives.hpp:53
Defines classes and functions for making classes creatable from input files.
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:29
The densitized rest-mass density .
Definition: Tags.hpp:29
Defines a list of useful type aliases for tensors.
The cutoff below which is set to MinimumValueOfD
Definition: FixConservatives.hpp:65
Stores a collection of function values.
Definition: DataVector.hpp:42
Information about the nested operations being performed by the parser, for use in printing errors...
Definition: Options.hpp:38
Wraps the template metaprogramming library used (brigand)
The minimum value of the rest mass density times the Lorentz factor, .
Definition: FixConservatives.hpp:57
Defines functions and classes from the GSL.
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
The densitized energy density .
Definition: Tags.hpp:35
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182