Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
9 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"
10 : #include "Options/Context.hpp"
11 : #include "Options/String.hpp"
12 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
13 : #include "PointwiseFunctions/Hydro/MagneticFieldTreatment.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 : /// \endcond
23 :
24 : namespace grmhd {
25 : namespace ValenciaDivClean {
26 :
27 : /*!
28 : * \ingroup VariableFixingGroup
29 : * \brief Fix conservative variables using method developed by Foucart.
30 : *
31 : * Adjusts the conservative variables as follows:
32 : * - If the electron fraction \f$Y_e\f$ is below the value of the option
33 : * `CutoffYe`, change \f$\tilde{Y}_e\f$ to \f$\tilde{D}\f$ times the value of
34 : * the option `MinimumValueOfYe`.
35 : * - Changes \f${\tilde D}\f$, the generalized mass-energy density, such
36 : * that \f$D\f$, the product of the rest mass density \f$\rho\f$ and the
37 : * Lorentz factor \f$W\f$, is set to value of the option `MinimumValueOfD`,
38 : * whenever \f$D\f$ is below the value of the option `CutoffD`.
39 : * - Increases \f${\tilde \tau}\f$, the generalized internal energy density,
40 : * such that
41 : * \f${\tilde B}^2 \leq 2 \sqrt{\gamma} (1 - \epsilon_B) {\tilde \tau}\f$,
42 : * where \f${\tilde B}^i\f$ is the generalized magnetic field,
43 : * \f$\gamma\f$ is the determinant of the spatial metric, and
44 : * \f$\epsilon_B\f$ is the option `SafetyFactorForB`.
45 : * - Decreases \f${\tilde S}_i\f$, the generalized momentum density, such that
46 : * \f${\tilde S}^2 \leq (1 - \epsilon_S) {\tilde S}^2_{max}\f$, where
47 : * \f$\epsilon_S\f$ is the option `SafetyFactorForS`, and
48 : * \f${\tilde S}^2_{max}\f$ is a complicated function of the conservative
49 : * variables which can only be found through root finding. There are
50 : * sufficient conditions for a set of conservative variables to satisfy the
51 : * inequality, which can be used to avoid root finding at most points.
52 : *
53 : * \note The routine currently assumes the minimum specific enthalpy is one.
54 : *
55 : * For more details see Appendix B from the [thesis of Francois
56 : * Foucart](https://ecommons.cornell.edu/handle/1813/30652)
57 : *
58 : * You can plot the function whose root we are finding using:
59 : *
60 : * \code{.py}
61 : *
62 : * import numpy as np
63 : * import matplotlib.pyplot as plt
64 : *
65 : * upper_bound = 1.000119047987896748e+00
66 : * lower_bound = 1.000000000000000000e+00
67 : * s_tilde_squared = 5.513009056734747750e-30
68 : * d_tilde = 1.131468709980503465e-12
69 : * sqrt_det_g: 1.131468709980503418e+00
70 : * tau_tilde = 1.346990732914080573e-16
71 : * b_tilde_squared = 3.048155733848927391e-16
72 : * b_squared_over_d = 2.380959757934347320e-04
73 : * tau_over_d = 1.190479878968363843e-04
74 : * normalized_s_dot_b = -9.999999082462245337e-01
75 : *
76 : *
77 : * def function_of_w(lorentz_factor):
78 : * return ((lorentz_factor + b_squared_over_d - tau_over_d - 1.0) *
79 : * (lorentz_factor**2 + b_squared_over_d * normalized_s_dot_b**2 *
80 : * (b_squared_over_d + 2.0 * lorentz_factor)) -
81 : * 0.5 * b_squared_over_d -
82 : * 0.5 * b_squared_over_d * normalized_s_dot_b**2 *
83 : * (lorentz_factor**2 - 1.0 +
84 : * 2.0 * lorentz_factor * b_squared_over_d + b_squared_over_d**2))
85 : *
86 : *
87 : * lorentz_factor = np.linspace(lower_bound, upper_bound, num=10000)
88 : *
89 : * plt.plot(lorentz_factor, function_of_w(lorentz_factor))
90 : * plt.show()
91 : *
92 : * \endcode
93 : */
94 1 : class FixConservatives {
95 : public:
96 : /// \brief Minimum value of rest-mass density times lorentz factor
97 1 : struct MinimumValueOfD {
98 0 : using type = double;
99 0 : static type lower_bound() { return 0.0; }
100 0 : static constexpr Options::String help = {
101 : "Minimum value of rest-mass density times lorentz factor"};
102 : };
103 : /// \brief Cutoff below which \f$D = \rho W\f$ is set to MinimumValueOfD
104 1 : struct CutoffD {
105 0 : using type = double;
106 0 : static type lower_bound() { return 0.0; }
107 0 : static constexpr Options::String help = {
108 : "Cutoff below which D is set to MinimumValueOfD"};
109 : };
110 : /// \brief Minimum value of electron fraction \f$Y_e\f$
111 1 : struct MinimumValueOfYe {
112 0 : using type = double;
113 0 : static type lower_bound() { return 0.0; }
114 0 : static constexpr Options::String help = {
115 : "Minimum value of electron fraction"};
116 : };
117 : /// \brief Cutoff below which \f$Y_e\f$ is set to MinimumValueOfYe
118 1 : struct CutoffYe {
119 0 : using type = double;
120 0 : static type lower_bound() { return 0.0; }
121 0 : static constexpr Options::String help = {
122 : "Cutoff below which Y_e is set to MinimumValueOfYe"};
123 : };
124 : /// \brief Safety factor \f$\epsilon_B\f$.
125 1 : struct SafetyFactorForB {
126 0 : using type = double;
127 0 : static type lower_bound() { return 0.0; }
128 0 : static constexpr Options::String help = {
129 : "Safety factor for magnetic field bound."};
130 : };
131 : /// \brief Safety factor \f$\epsilon_S\f$.
132 1 : struct SafetyFactorForS {
133 0 : using type = double;
134 0 : static type lower_bound() { return 0.0; }
135 0 : static constexpr Options::String help = {
136 : "Safety factor for momentum density bound above density cutoff."};
137 : };
138 : /// \brief Cutoff in \f$\rho_0 W\f$ below which we use a stricter safety
139 : /// factor for the magnitude of S.
140 1 : struct SafetyFactorForSCutoffD {
141 0 : using type = double;
142 0 : static type lower_bound() { return 0.0; }
143 0 : static constexpr Options::String help = {
144 : "Below this value of rest mass density time Lorentz factor, limit S "
145 : "more agressively."};
146 : };
147 :
148 : /// \brief Below SafetyFactorForSCutoffD, reduce \f$\epsilon_S\f$ by
149 : /// SafetyFactorForSSlope times
150 : /// \f$\log_{10}(\rho_0 W / SafetyFactorForSCutoffD)\f$
151 1 : struct SafetyFactorForSSlope {
152 0 : using type = double;
153 0 : static type lower_bound() { return 0.0; }
154 0 : static constexpr Options::String help = {
155 : "Slope of safety factor for momentum density bound below "
156 : "SafetyFactorForSCutoffD, express as a function of log10(rho*W)."};
157 : };
158 : /// Whether or not the limiting is enabled
159 1 : struct Enable {
160 0 : using type = bool;
161 0 : static constexpr Options::String help = {
162 : "If true then the limiting is applied."};
163 : };
164 : /// How to treat the magnetic field
165 1 : struct MagneticField {
166 0 : using type = hydro::MagneticFieldTreatment;
167 0 : static constexpr Options::String help = {
168 : "How to treat the magnetic field."};
169 : };
170 :
171 0 : using options =
172 : tmpl::list<MinimumValueOfD, CutoffD, MinimumValueOfYe, CutoffYe,
173 : SafetyFactorForB, SafetyFactorForS, SafetyFactorForSCutoffD,
174 : SafetyFactorForSSlope, Enable, MagneticField>;
175 0 : static constexpr Options::String help = {
176 : "Variable fixing used in Foucart's thesis.\n"};
177 :
178 0 : FixConservatives(double minimum_rest_mass_density_times_lorentz_factor,
179 : double rest_mass_density_times_lorentz_factor_cutoff,
180 : double minimum_electron_fraction,
181 : double electron_fraction_cutoff,
182 : double safety_factor_for_magnetic_field,
183 : double safety_factor_for_momentum_density,
184 : double safety_factor_for_momentum_density_cutoff_d,
185 : double safety_factor_for_momentum_density_slope, bool enable,
186 : hydro::MagneticFieldTreatment magnetic_field_treatment,
187 : const Options::Context& context = {});
188 :
189 0 : FixConservatives() = default;
190 0 : FixConservatives(const FixConservatives& /*rhs*/) = default;
191 0 : FixConservatives& operator=(const FixConservatives& /*rhs*/) = default;
192 0 : FixConservatives(FixConservatives&& /*rhs*/) = default;
193 0 : FixConservatives& operator=(FixConservatives&& /*rhs*/) = default;
194 0 : ~FixConservatives() = default;
195 :
196 : // NOLINTNEXTLINE(google-runtime-references)
197 0 : void pup(PUP::er& p);
198 :
199 0 : using return_tags = tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
200 : grmhd::ValenciaDivClean::Tags::TildeYe,
201 : grmhd::ValenciaDivClean::Tags::TildeTau,
202 : grmhd::ValenciaDivClean::Tags::TildeS<>>;
203 0 : using argument_tags =
204 : tmpl::list<grmhd::ValenciaDivClean::Tags::TildeB<>,
205 : gr::Tags::SpatialMetric<DataVector, 3>,
206 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
207 : gr::Tags::SqrtDetSpatialMetric<DataVector>>;
208 :
209 : /// Returns `true` if any variables were fixed.
210 1 : bool operator()(
211 : gsl::not_null<Scalar<DataVector>*> tilde_d,
212 : gsl::not_null<Scalar<DataVector>*> tilde_ye,
213 : gsl::not_null<Scalar<DataVector>*> tilde_tau,
214 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> tilde_s,
215 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
216 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
217 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
218 : const Scalar<DataVector>& sqrt_det_spatial_metric) const;
219 :
220 : private:
221 0 : friend bool operator==(const FixConservatives& lhs,
222 : const FixConservatives& rhs);
223 :
224 0 : double minimum_rest_mass_density_times_lorentz_factor_{
225 : std::numeric_limits<double>::signaling_NaN()};
226 0 : double rest_mass_density_times_lorentz_factor_cutoff_{
227 : std::numeric_limits<double>::signaling_NaN()};
228 0 : double minimum_electron_fraction_{
229 : std::numeric_limits<double>::signaling_NaN()};
230 0 : double electron_fraction_cutoff_{
231 : std::numeric_limits<double>::signaling_NaN()};
232 0 : double one_minus_safety_factor_for_magnetic_field_{
233 : std::numeric_limits<double>::signaling_NaN()};
234 0 : double one_minus_safety_factor_for_momentum_density_{
235 : std::numeric_limits<double>::signaling_NaN()};
236 0 : double safety_factor_for_momentum_density_cutoff_d_{
237 : std::numeric_limits<double>::signaling_NaN()};
238 0 : double safety_factor_for_momentum_density_slope_{
239 : std::numeric_limits<double>::signaling_NaN()};
240 0 : bool enable_{true};
241 0 : hydro::MagneticFieldTreatment magnetic_field_treatment_{
242 : hydro::MagneticFieldTreatment::AssumeNonZero};
243 : };
244 :
245 0 : bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs);
246 :
247 : } // namespace ValenciaDivClean
248 : } // namespace grmhd
|
