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 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
9 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservativeOptions.hpp"
10 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"
11 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
12 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
13 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
14 : #include "Utilities/TMPL.hpp"
15 :
16 : /// \cond
17 : namespace gsl {
18 : template <typename T>
19 : class not_null;
20 : } // namespace gsl
21 :
22 : class DataVector;
23 : /// \endcond
24 :
25 : namespace grmhd {
26 : namespace ValenciaDivClean {
27 : /*!
28 : * \brief Compute the primitive variables from the conservative variables
29 : *
30 : * For the Valencia formulation of the GRMHD system with divergence cleaning,
31 : * the conversion of the evolved conserved variables to the primitive variables
32 : * cannot be expressed in closed analytic form and requires a root find.
33 : *
34 : * [Siegel {\em et al}, The Astrophysical Journal 859:71(2018)]
35 : * (http://iopscience.iop.org/article/10.3847/1538-4357/aabcc5/meta)
36 : * compares several inversion methods.
37 : *
38 : * If `ErrorOnFailure` is `false` then the returned `bool` will be `false` if
39 : * recovery failed and `true` if it succeeded.
40 : *
41 : * If `EnforcePhysicality` is `false` then the hydrodynamic inversion will
42 : * return with an error if the input conservatives are unphysical, i.e., if no
43 : * solution exits. The exact behavior is governed by `ErrorOnFailure`.
44 : *
45 : * \note
46 : * Here, we outline the algebra for computing the spatial velocity
47 : * \f$ v^i \f$:
48 : *
49 : * We start from definition of the momentum density \f$ S_i \f$:
50 : * \f[
51 : * \begin{aligned}
52 : * S_i &= (\rho h)^* W^2 v_i - \alpha b^0 b_i \\
53 : * &= \left(\rho h + \frac{B^2}{W^2} + (v \cdot B)^2\right) W^2 v_i
54 : * - \alpha \left(\frac{W}{\alpha}(v \cdot B)\right)
55 : * \left(\frac{B_i}{W} + (v \cdot B) W v_i\right) \\
56 : * &= \left(\rho h + \frac{B^2}{W^2}\right) W^2 v_i - (v \cdot B) B_i
57 : * \end{aligned}
58 : * \f]
59 : *
60 : * From the above, we obtain:
61 : * \f[
62 : * \begin{aligned}
63 : * S^i &= \left(\rho h + \frac{B^2}{W^2}\right) W^2 v^i - (v \cdot B) B^i \\
64 : * S \cdot B &= \left(\rho h + \frac{B^2}{W^2}\right) W^2 (v \cdot B)
65 : * - (v \cdot B) B^2 \\
66 : * &= \rho h W^2 (v \cdot B) \\
67 : * v^i &= \frac{1}{\sqrt{\gamma} (\rho h W^2 + B^2)} \tilde{S}^i
68 : * + \frac{S \cdot B}{\rho h W^2 (\rho h W^2 + B^2)} B^i
69 : * \end{aligned}
70 : * \f]
71 : */
72 : template <typename OrderedListOfPrimitiveRecoverySchemes,
73 : bool ErrorOnFailure = true>
74 1 : struct PrimitiveFromConservative {
75 0 : using return_tags =
76 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
77 : hydro::Tags::ElectronFraction<DataVector>,
78 : hydro::Tags::SpecificInternalEnergy<DataVector>,
79 : hydro::Tags::SpatialVelocity<DataVector, 3>,
80 : hydro::Tags::MagneticField<DataVector, 3>,
81 : hydro::Tags::DivergenceCleaningField<DataVector>,
82 : hydro::Tags::LorentzFactor<DataVector>,
83 : hydro::Tags::Pressure<DataVector>,
84 : hydro::Tags::Temperature<DataVector>>;
85 :
86 0 : using argument_tags = tmpl::list<
87 : grmhd::ValenciaDivClean::Tags::TildeD,
88 : grmhd::ValenciaDivClean::Tags::TildeYe,
89 : grmhd::ValenciaDivClean::Tags::TildeTau,
90 : grmhd::ValenciaDivClean::Tags::TildeS<>,
91 : grmhd::ValenciaDivClean::Tags::TildeB<>,
92 : grmhd::ValenciaDivClean::Tags::TildePhi,
93 : gr::Tags::SpatialMetric<DataVector, 3>,
94 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
95 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
96 : hydro::Tags::GrmhdEquationOfState,
97 : grmhd::ValenciaDivClean::Tags::PrimitiveFromConservativeOptions>;
98 :
99 : template <bool EnforcePhysicality = true>
100 0 : static bool apply(
101 : gsl::not_null<Scalar<DataVector>*> rest_mass_density,
102 : gsl::not_null<Scalar<DataVector>*> electron_fraction,
103 : gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
104 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> spatial_velocity,
105 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> magnetic_field,
106 : gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
107 : gsl::not_null<Scalar<DataVector>*> lorentz_factor,
108 : gsl::not_null<Scalar<DataVector>*> pressure,
109 : gsl::not_null<Scalar<DataVector>*> temperature,
110 : const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
111 : const Scalar<DataVector>& tilde_tau,
112 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
113 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
114 : const Scalar<DataVector>& tilde_phi,
115 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
116 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
117 : const Scalar<DataVector>& sqrt_det_spatial_metric,
118 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state,
119 : const grmhd::ValenciaDivClean::PrimitiveFromConservativeOptions&
120 : primitive_from_conservative_options);
121 :
122 : private:
123 : template <bool EnforcePhysicality, typename EosType>
124 0 : static bool impl(
125 : gsl::not_null<Scalar<DataVector>*> rest_mass_density,
126 : gsl::not_null<Scalar<DataVector>*> electron_fraction,
127 : gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
128 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> spatial_velocity,
129 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> magnetic_field,
130 : gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
131 : gsl::not_null<Scalar<DataVector>*> lorentz_factor,
132 : gsl::not_null<Scalar<DataVector>*> pressure,
133 : gsl::not_null<Scalar<DataVector>*> temperature,
134 : const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
135 : const Scalar<DataVector>& tilde_tau,
136 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
137 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
138 : const Scalar<DataVector>& tilde_phi,
139 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
140 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
141 : const Scalar<DataVector>& sqrt_det_spatial_metric,
142 : const EosType& equation_of_state,
143 : const grmhd::ValenciaDivClean::PrimitiveFromConservativeOptions&
144 : primitive_from_conservative_options);
145 :
146 : // Use Kastaun hydro inversion if B is dynamically unimportant
147 0 : static constexpr bool use_hydro_optimization = true;
148 : };
149 : } // namespace ValenciaDivClean
150 : } // namespace grmhd
|
