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 <iosfwd>
8 : #include <limits>
9 : #include <optional>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Options/Auto.hpp"
13 : #include "Options/Context.hpp"
14 : #include "Options/String.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
16 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
17 : #include "PointwiseFunctions/Hydro/Tags.hpp"
18 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 :
25 : namespace Options {
26 : class Option;
27 : template <typename T>
28 : struct create_from_yaml;
29 : } // namespace Options
30 :
31 : namespace PUP {
32 : class er;
33 : } // namespace PUP
34 : /// \endcond
35 :
36 : namespace VariableFixing {
37 : /// \brief How to apply atmosphere to the reconstructed state.
38 0 : enum class FixReconstructedStateToAtmosphere : uint8_t {
39 : Always,
40 : AtDgFdInterfaceOnly,
41 : OnFdOnly,
42 : Never
43 : };
44 :
45 0 : std::ostream& operator<<(std::ostream& os,
46 : const FixReconstructedStateToAtmosphere& t);
47 : } // namespace VariableFixing
48 :
49 : template <>
50 0 : struct Options::create_from_yaml<
51 : VariableFixing::FixReconstructedStateToAtmosphere> {
52 : template <typename Metavariables>
53 0 : static VariableFixing::FixReconstructedStateToAtmosphere create(
54 : const Options::Option& options) {
55 : return create<void>(options);
56 : }
57 : };
58 :
59 : template <>
60 0 : VariableFixing::FixReconstructedStateToAtmosphere
61 : Options::create_from_yaml<VariableFixing::FixReconstructedStateToAtmosphere>::
62 : create<void>(const Options::Option& options);
63 :
64 : namespace VariableFixing {
65 : /*!
66 : * \ingroup VariableFixingGroup
67 : * \brief Fix the primitive variables to an atmosphere in low density regions
68 : *
69 : * If the rest mass density is below \f$\rho_{\textrm{cutoff}}\f$
70 : * (DensityCutoff), it is set to \f$\rho_{\textrm{atm}}\f$
71 : * (DensityOfAtmosphere), and the pressure, and specific internal energy (for
72 : * one-dimensional equations of state) are adjusted to satisfy the equation of
73 : * state. For a two-dimensional equation of state, the specific internal energy
74 : * is set to zero.
75 : */
76 : template <size_t Dim>
77 1 : class FixToAtmosphere {
78 0 : struct Disabled {};
79 0 : struct FromEos {};
80 :
81 : public:
82 : /// \brief Rest mass density of the atmosphere
83 1 : struct DensityOfAtmosphere {
84 0 : using type = double;
85 0 : static type lower_bound() { return 0.0; }
86 0 : static constexpr Options::String help = {"Density of atmosphere"};
87 : };
88 : /// \brief Rest mass density at which to impose the atmosphere. Should be
89 : /// greater than or equal to the density of the atmosphere.
90 1 : struct DensityCutoff {
91 0 : using type = double;
92 0 : static type lower_bound() { return 0.0; }
93 0 : static constexpr Options::String help = {
94 : "Density to impose atmosphere at. Must be >= rho_atm"};
95 : };
96 :
97 : /*!
98 : * \brief Limit the velocity in and near the atmosphere.
99 : *
100 : * Let $v_{\max}$ be the maximum magnitude of the
101 : * velocity near the atmosphere, which we typically set to $10^{-4}$.
102 : * We let $v_{\mathrm{atm}}$ be the maximum magnitude of the velocity
103 : * in the atmosphere, which we typically set to $0$. We then define
104 : * the maximum magnitude of the spatial velocity to be
105 : *
106 : * \f{align*}{
107 : * \tilde{v}
108 : * &=\begin{cases}
109 : * v_{\mathrm{atm}}, & \mathrm{if}\; \rho < \rho_{v_-} \ \
110 : * v_{\mathrm{atm}} + \left(v_{\max} - v_{\mathrm{atm}}\right)
111 : * \frac{\rho - \rho_{v_-}}{\rho_{v_+} - \rho_{v_-}},
112 : * & \mathrm{if}\;\rho_{v_-} \le \rho < \rho_{v_+}
113 : * \end{cases}
114 : * \f}
115 : *
116 : * We then rescale the velocity by
117 : *
118 : * \f{align*}{
119 : * v^i\to v^i\frac{\tilde{v}}{\sqrt{v^i\gamma_{ij}v^j}}.
120 : * \f}
121 : */
122 1 : struct VelocityLimitingOptions {
123 0 : struct AtmosphereMaxVelocity {
124 0 : using type = double;
125 0 : static constexpr Options::String help = {
126 : "The maximum velocity magnitude IN the atmosphere. Typically set to "
127 : "0."};
128 : };
129 :
130 0 : struct NearAtmosphereMaxVelocity {
131 0 : using type = double;
132 0 : static constexpr Options::String help = {
133 : "The maximum velocity magnitude NEAR the atmosphere. Typically set "
134 : "to 1e-4."};
135 : };
136 :
137 0 : struct AtmosphereDensityCutoff {
138 0 : using type = double;
139 0 : static constexpr Options::String help = {
140 : "The rest mass density cutoff below which the velocity magnitude is "
141 : "limited to AtmosphereMaxVelocity. Typically set to "
142 : "(10 or 20)*DensityOfAtmosphere."};
143 : };
144 :
145 0 : struct TransitionDensityBound {
146 0 : using type = double;
147 0 : static constexpr Options::String help = {
148 : "The rest mass density above which no velocity limiting is done. "
149 : "Between "
150 : "this value and AtmosphereDensityCutoff a linear transition in the "
151 : "maximum magnitude of the velocity between AtmosphereMaxVelocity and "
152 : "NearAtmosphereMaxVelocity is done. Typically set to "
153 : "10*AtmosphereDensityCutoff."};
154 : };
155 0 : using options = tmpl::list<AtmosphereMaxVelocity, NearAtmosphereMaxVelocity,
156 : AtmosphereDensityCutoff, TransitionDensityBound>;
157 0 : static constexpr Options::String help = {
158 : "Limit the velocity in and near the atmosphere."};
159 :
160 : // NOLINTNEXTLINE(google-runtime-references)
161 0 : void pup(PUP::er& p);
162 :
163 0 : bool operator==(const VelocityLimitingOptions& rhs) const;
164 0 : bool operator!=(const VelocityLimitingOptions& rhs) const;
165 :
166 0 : double atmosphere_max_velocity{
167 : std::numeric_limits<double>::signaling_NaN()};
168 0 : double near_atmosphere_max_velocity{
169 : std::numeric_limits<double>::signaling_NaN()};
170 0 : double atmosphere_density_cutoff{
171 : std::numeric_limits<double>::signaling_NaN()};
172 0 : double transition_density_bound{
173 : std::numeric_limits<double>::signaling_NaN()};
174 : };
175 0 : struct VelocityLimiting {
176 0 : using type = Options::Auto<VelocityLimitingOptions, Disabled>;
177 0 : static constexpr Options::String help = VelocityLimitingOptions::help;
178 : };
179 :
180 : /*!
181 : * \brief Options for limiting the temperature in the atmosphere by
182 : * effectively limiting the polytropic constant, with a generalization for
183 : * finite temperature equations of state.
184 : *
185 : * The basic density limit is fine for a cold EOS, but when a finite
186 : * temperature EOS is used the temperature needs to be controlled in the
187 : * atmosphere. While the basic bound of
188 : * $T\in[T_{\mathrm{min}},T_{\mathrm{max}}]$ is a necessary condition, it is
189 : * often not sufficient for stability in the atmosphere. We define lower and
190 : * upper bounds $\rho_{\tau_-}$ and $\rho_{\tau_+}$ for $\rho$ where we apply
191 : * different limiting behavior when $\rho<\rho_{\tau_-}$,
192 : * $\rho\in[\rho_{\tau_-},\rho_{\tau_+}]$, and $\rho>\rho_{\tau_+}$.
193 : *
194 : * When $\rho<\rho_{\tau_-}$ we set $T=T_{\mathrm{min}}$ if
195 : * $|T-T_{\mathrm{min}}|>\epsilon_{\kappa_-} |T|$, otherwise we don't change
196 : * $T$. SpEC uses $\epsilon_{\kappa_-}=10^{-3}$.
197 : *
198 : * When $\rho\in[\rho_{\tau_-},\rho_{\tau_+}]$ we apply a more complicated
199 : * algorithm. We define
200 : *
201 : * \f{align*}{
202 : * \kappa_{\mathrm{max}} =
203 : * 1 +
204 : * \epsilon_{\kappa,\mathrm{max}}
205 : * \min\left(\frac{\rho-\rho_{\tau_-}}{\rho_{\tau_+}-\rho_{\tau_-}},1\right)
206 : * \f}
207 : *
208 : * where $\epsilon_{\kappa,\mathrm{max}}\ge0$ (set to 0.01 in SpEC). With
209 : * $\kappa_{\mathrm{max}}$ computed we now limit the temperature. We define
210 : * $p_{\mathrm{min}}=p(\rho,T_{\mathrm{min}},Y_e)$, $p=p(\rho,T,Y_e)$, and
211 : * $p_{\kappa_{\mathrm{max}}}=p_{\mathrm{min}}\kappa_{\mathrm{max}}$. We then
212 : * need to find $T$ such that $p(\rho, T, Y_e) = p_{\kappa_{\mathrm{max}}}$.
213 : * We do this by finding the root of
214 : *
215 : * \f{align*}{
216 : * f(T')=p(\rho, T', Y_e) - p_{\kappa_{\mathrm{max}}}.
217 : * \f}
218 : *
219 : * We bracket the temperature $T'\in[T_{\min}, T]$, where $T$ is the
220 : * unmodified temperature. We solve the equation using the TOMS748 algorithm
221 : * to avoid the need for derivatives of the pressure with respect to the
222 : * temperature. If $T-T_{\min}<10^{-13}$ we set $T'=\frac{1}{2}(T+T_{\min})$
223 : * to avoid any root find.
224 : *
225 : * When $\rho>\rho_{\tau_+}$ we optionally can apply the same procedure as in
226 : * the previous case, but with
227 : * $\kappa_\mathrm{max}=1+\epsilon_{\kappa,\mathrm{max}}$. This is not
228 : * typically done.
229 : */
230 1 : struct KappaLimitingOptions {
231 0 : struct DensityLowerBound {
232 0 : using type = double;
233 0 : static type lower_bound() { return 0.0; }
234 0 : static constexpr Options::String help = {
235 : "Below this value we set T=T_min if |T-T_min|<EplisonKappaMinus "
236 : "|T|. Typically set to about a factor of 20 larger than the "
237 : "atmosphere density."};
238 : };
239 0 : struct EplisonKappaMinus {
240 0 : using type = double;
241 0 : static constexpr Options::String help = {
242 : "Used to limit the temperature in conjunction with "
243 : "DensityLowerBound. See that text for the formula. Typically set "
244 : "to 1.0e-3."};
245 : };
246 0 : struct DensityUpperBound {
247 0 : using type = double;
248 0 : static type lower_bound() { return 0.0; }
249 0 : static constexpr Options::String help = {
250 : "The density below which we limit the temperature by limiting the "
251 : "pressure (at fixed rho, Y_e) to 'p(T)<p_min kappa_max'. Typically "
252 : "set a factor of 10 larger than DensityLowerBound."};
253 : };
254 0 : struct EpsilonKappaMax {
255 0 : using type = double;
256 0 : static constexpr Options::String help = {
257 : "The epsilon factor in the KappaMax computation multiplying the "
258 : "transition factor. Typically set to 0.01."};
259 : };
260 0 : struct MinTemperature {
261 0 : using type = Options::Auto<double, FromEos>;
262 0 : static constexpr Options::String help = {"The minimum temperature."};
263 : };
264 0 : struct LimitAboveDensityUpperBound {
265 0 : using type = bool;
266 0 : static constexpr Options::String help = {
267 : "If true then we limit the temperature using the KappaMax procedure "
268 : "at all densities above DensityUpperBound, but with "
269 : "`KappaMax=1+EplisonKappaMax`. Typically set to False."};
270 : };
271 0 : using options = tmpl::list<DensityLowerBound, EplisonKappaMinus,
272 : DensityUpperBound, EpsilonKappaMax,
273 : MinTemperature, LimitAboveDensityUpperBound>;
274 0 : static constexpr Options::String help = {
275 : "If set then we apply a limiting precodure on the temperature near the "
276 : "atmosphere based on essentially limiting the polytropic constant in a "
277 : "Gamma-law equation of state."};
278 :
279 : // NOLINTNEXTLINE(google-runtime-references)
280 0 : void pup(PUP::er& p);
281 :
282 0 : bool operator==(const KappaLimitingOptions& rhs) const;
283 0 : bool operator!=(const KappaLimitingOptions& rhs) const;
284 :
285 0 : double density_lower_bound{std::numeric_limits<double>::signaling_NaN()};
286 0 : double eplison_kappa_minus{std::numeric_limits<double>::signaling_NaN()};
287 0 : double density_upper_bound{std::numeric_limits<double>::signaling_NaN()};
288 0 : double epsilon_kappa_max{std::numeric_limits<double>::signaling_NaN()};
289 0 : std::optional<double> min_temperature{std::nullopt};
290 0 : bool limit_above_density_upper_bound{false};
291 : };
292 : /// \brief If set then we apply a limiting precodure on the temperature near
293 : /// the atmosphere based on essentially limiting the polytropic constant in a
294 : /// Gamma-law equation of state.
295 1 : struct KappaLimiting {
296 0 : using type = Options::Auto<KappaLimitingOptions, Disabled>;
297 0 : static constexpr Options::String help = KappaLimitingOptions::help;
298 : };
299 :
300 0 : using options = tmpl::list<DensityOfAtmosphere, DensityCutoff,
301 : VelocityLimiting, KappaLimiting>;
302 0 : static constexpr Options::String help = {
303 : "If the rest mass density is below DensityCutoff, it is set\n"
304 : "to DensityOfAtmosphere, and the pressure, and specific internal energy\n"
305 : "(for one-dimensional equations of state) are\n"
306 : "adjusted to satisfy the equation of state. For a two-dimensional\n"
307 : "equation of state, the specific internal energy is set to zero.\n"
308 : "In addition, the spatial velocity is set to zero, and the Lorentz\n"
309 : "factor is set to one.\n"};
310 :
311 0 : FixToAtmosphere(double density_of_atmosphere, double density_cutoff,
312 : std::optional<VelocityLimitingOptions> velocity_limiting,
313 : std::optional<KappaLimitingOptions> kappa_limiting,
314 : const Options::Context& context = {});
315 :
316 0 : FixToAtmosphere() = default;
317 0 : FixToAtmosphere(const FixToAtmosphere& /*rhs*/) = default;
318 0 : FixToAtmosphere& operator=(const FixToAtmosphere& /*rhs*/) = default;
319 0 : FixToAtmosphere(FixToAtmosphere&& /*rhs*/) = default;
320 0 : FixToAtmosphere& operator=(FixToAtmosphere&& /*rhs*/) = default;
321 0 : ~FixToAtmosphere() = default;
322 :
323 : // NOLINTNEXTLINE(google-runtime-references)
324 0 : void pup(PUP::er& p);
325 :
326 0 : using return_tags =
327 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
328 : hydro::Tags::SpecificInternalEnergy<DataVector>,
329 : hydro::Tags::SpatialVelocity<DataVector, Dim>,
330 : hydro::Tags::LorentzFactor<DataVector>,
331 : hydro::Tags::Pressure<DataVector>,
332 : hydro::Tags::Temperature<DataVector>>;
333 0 : using argument_tags = tmpl::list<hydro::Tags::ElectronFraction<DataVector>,
334 : gr::Tags::SpatialMetric<DataVector, Dim>,
335 : hydro::Tags::GrmhdEquationOfState>;
336 :
337 : // for use in `db::mutate_apply`
338 : template <size_t ThermodynamicDim>
339 0 : void operator()(
340 : gsl::not_null<Scalar<DataVector>*> rest_mass_density,
341 : gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
342 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
343 : spatial_velocity,
344 : gsl::not_null<Scalar<DataVector>*> lorentz_factor,
345 : gsl::not_null<Scalar<DataVector>*> pressure,
346 : gsl::not_null<Scalar<DataVector>*> temperature,
347 : const Scalar<DataVector>& electron_fraction,
348 : const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
349 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
350 : equation_of_state) const;
351 :
352 : /// @{
353 : /// \brief Other algorithmic decisions may depend on the atmosphere treatment
354 : /// so provide access to the values.
355 1 : double density_of_atmosphere() const { return density_of_atmosphere_; }
356 1 : double density_cutoff() const { return density_cutoff_; }
357 1 : const std::optional<VelocityLimitingOptions>& velocity_limiting() const {
358 : return velocity_limiting_;
359 : }
360 1 : const std::optional<KappaLimitingOptions>& kappa_limiting() const {
361 : return kappa_limiting_;
362 : }
363 : /// @}
364 :
365 : private:
366 : template <size_t ThermodynamicDim>
367 0 : void set_density_to_atmosphere(
368 : gsl::not_null<Scalar<DataVector>*> rest_mass_density,
369 : gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
370 : gsl::not_null<Scalar<DataVector>*> temperature,
371 : gsl::not_null<Scalar<DataVector>*> pressure,
372 : const Scalar<DataVector>& electron_fraction,
373 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
374 : equation_of_state,
375 : size_t grid_index) const;
376 :
377 0 : void apply_velocity_limit(
378 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
379 : spatial_velocity,
380 : gsl::not_null<Scalar<DataVector>*> lorentz_factor,
381 : const Scalar<DataVector>& rest_mass_density,
382 : const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
383 : size_t grid_index) const;
384 :
385 : template <size_t ThermodynamicDim>
386 0 : bool apply_kappa_limit(
387 : gsl::not_null<Scalar<DataVector>*> temperature,
388 : const Scalar<DataVector>& rest_mass_density,
389 : const Scalar<DataVector>& electron_fraction,
390 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
391 : equation_of_state,
392 : size_t grid_index) const;
393 :
394 : template <size_t SpatialDim>
395 : // NOLINTNEXTLINE(readability-redundant-declaration)
396 0 : friend bool operator==(const FixToAtmosphere<SpatialDim>& lhs,
397 : const FixToAtmosphere<SpatialDim>& rhs);
398 :
399 0 : double density_of_atmosphere_{std::numeric_limits<double>::signaling_NaN()};
400 0 : double density_cutoff_{std::numeric_limits<double>::signaling_NaN()};
401 0 : std::optional<VelocityLimitingOptions> velocity_limiting_{std::nullopt};
402 0 : std::optional<KappaLimitingOptions> kappa_limiting_{std::nullopt};
403 : };
404 :
405 : template <size_t Dim>
406 0 : bool operator!=(const FixToAtmosphere<Dim>& lhs,
407 : const FixToAtmosphere<Dim>& rhs);
408 :
409 : } // namespace VariableFixing
|
