Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 :
9 : #include "DataStructures/DataBox/Tag.hpp"
10 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Domain/FaceNormal.hpp"
13 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
14 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
16 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : /// \cond
21 : class DataVector;
22 : namespace Tags {
23 : template <typename Tag>
24 : struct Normalized;
25 : } // namespace Tags
26 : /// \endcond
27 :
28 : namespace grmhd {
29 : namespace ValenciaDivClean {
30 : /// @{
31 : /*!
32 : * \brief Compute the characteristic speeds for the Valencia formulation of
33 : * GRMHD with divergence cleaning.
34 : *
35 : * Obtaining the exact form of the characteristic speeds involves the solution
36 : * of a nontrivial quartic equation for the fast and slow modes. Here we make
37 : * use of a common approximation in the literature (e.g. \cite Gammie2003)
38 : * where the resulting characteristic speeds are analogous to those of the
39 : * Valencia formulation of the 3-D relativistic Euler system
40 : * (see RelativisticEuler::Valencia::characteristic_speeds),
41 : *
42 : * \f{align*}
43 : * \lambda_2 &= \alpha \Lambda^- - \beta_n,\\
44 : * \lambda_{3, 4, 5, 6, 7} &= \alpha v_n - \beta_n,\\
45 : * \lambda_{8} &= \alpha \Lambda^+ - \beta_n,
46 : * \f}
47 : *
48 : * with the substitution
49 : *
50 : * \f{align*}
51 : * c_s^2 \longrightarrow c_s^2 + v_A^2(1 - c_s^2)
52 : * \f}
53 : *
54 : * in the definition of \f$\Lambda^\pm\f$. Here \f$v_A\f$ is the Alfvén
55 : * speed. In addition, two more speeds corresponding to the divergence cleaning
56 : * mode and the longitudinal magnetic field are added,
57 : *
58 : * \f{align*}
59 : * \lambda_1 = -\alpha - \beta_n,\\
60 : * \lambda_9 = \alpha - \beta_n.
61 : * \f}
62 : *
63 : * \note The ordering assumed here is such that, in the Newtonian limit,
64 : * the exact expressions for \f$\lambda_{2, 8}\f$, \f$\lambda_{3, 7}\f$,
65 : * and \f$\lambda_{4, 6}\f$ should reduce to the
66 : * corresponding fast modes, Alfvén modes, and slow modes, respectively.
67 : * See \cite Dedner2002 for a detailed description of the hyperbolic
68 : * characterization of Newtonian MHD. In terms of the primitive variables:
69 : *
70 : * \f{align*}
71 : * v^2 &= \gamma_{mn} v^m v^n \\
72 : * c_s^2 &= \frac{1}{h} \left[ \left( \frac{\partial p}{\partial \rho}
73 : * \right)_\epsilon +
74 : * \frac{p}{\rho^2} \left(\frac{\partial p}{\partial \epsilon}
75 : * \right)_\rho \right] \\
76 : * v_A^2 &= \frac{b^2}{b^2 + \rho h} \\
77 : * b^2 &= \frac{1}{W^2} \gamma_{mn} B^m B^n + \left( \gamma_{mn} B^m v^n
78 : * \right)^2
79 : * \f}
80 : *
81 : * where \f$\gamma_{mn}\f$ is the spatial metric, \f$\rho\f$ is the rest
82 : * mass density, \f$W = 1/\sqrt{1-v_i v^i}\f$ is the Lorentz factor, \f$h = 1 +
83 : * \epsilon + \frac{p}{\rho}\f$ is the specific enthalpy, \f$v^i\f$ is the
84 : * spatial velocity, \f$\epsilon\f$ is the specific internal energy, \f$p\f$ is
85 : * the pressure, and \f$B^i\f$ is the spatial magnetic field measured by an
86 : * Eulerian observer.
87 : */
88 :
89 : template <size_t ThermodynamicDim>
90 1 : std::array<DataVector, 9> characteristic_speeds_approximate_mhd(
91 : const Scalar<DataVector>& rest_mass_density,
92 : const Scalar<DataVector>& electron_fraction,
93 : const Scalar<DataVector>& specific_internal_energy,
94 : const Scalar<DataVector>& specific_enthalpy,
95 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
96 : const Scalar<DataVector>& lorentz_factor,
97 : const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
98 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, 3>& shift,
99 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
100 : const tnsr::i<DataVector, 3>& unit_normal,
101 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
102 : equation_of_state);
103 :
104 : template <size_t ThermodynamicDim>
105 1 : void characteristic_speeds_approximate_mhd(
106 : gsl::not_null<std::array<DataVector, 9>*> char_speeds,
107 : const Scalar<DataVector>& rest_mass_density,
108 : const Scalar<DataVector>& electron_fraction,
109 : const Scalar<DataVector>& specific_internal_energy,
110 : const Scalar<DataVector>& specific_enthalpy,
111 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
112 : const Scalar<DataVector>& lorentz_factor,
113 : const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
114 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, 3>& shift,
115 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
116 : const tnsr::i<DataVector, 3>& unit_normal,
117 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
118 : equation_of_state);
119 : /// @}
120 :
121 : /*!
122 : * \brief Labels for the characteristic speeds of the relativistic hydrodynamics
123 : * system.
124 : *
125 : * \see `grmhd::ValenciaDivClean::characteristic_speeds_hydro`
126 : */
127 1 : namespace HydroSpeed {
128 : /// The smallest speed, $\lambda_-$
129 1 : constexpr size_t LambdaMinus = 0;
130 : /// The normal velocity, $v_n = n_i v^i$
131 1 : constexpr size_t NormalDotVelocity = 1;
132 : /// The largest speed, $\lambda_+$
133 1 : constexpr size_t LambdaPlus = 2;
134 : } // namespace HydroSpeed
135 :
136 : /// @{
137 : /*!
138 : * \brief Compute the characteristic speeds for the relativistic hydrodynamics
139 : * system in the Eulerian frame.
140 : *
141 : * These are the eigenvalues of the flux Jacobian projected along a spatial
142 : * direction with unit normal \f$ s_i\f$, measured by an Eulerian observer.
143 : * They consist of four degenerate modes and two non-degenerate modes:
144 : *
145 : * \f{align*}
146 : * \lambda_\mathrm{deg} &= v_n ,\\
147 : * \lambda_{\pm} &=
148 : * \frac{ (1 - c_s^2)\,v_n \pm c_s\,d / W^2 }
149 : * { 1 - v^2 c_s^2 }.
150 : * \f}
151 : *
152 : * where
153 : *
154 : * \f{align*}
155 : * d &\equiv
156 : * W\,\sqrt{ 1 - v^2 c_s^2 - v_n^2(1 - c_s^2) } .
157 : * \f}
158 : *
159 : * The variables are defined as:
160 : *
161 : * \f{align*}
162 : * v^2 &= \gamma_{mn} v^m v^n ,\\
163 : * v_n &= v^a s_a ,\\
164 : * W &= \frac{1}{\sqrt{1 - v^a v_a}} ,\\
165 : * h &= 1 + \epsilon + \frac{p}{\rho} ,\\
166 : * c_s^2 &= \frac{1}{h}\left(\chi + \kappa \frac{p}{\rho^2}\right) ,
167 : * \f}
168 : *
169 : * The returned tensor has index 0..2, in increasing order (see `HydroSpeed`):
170 : * - `characteristic_speeds.get(0)`: acoustic speed \f$\lambda_-\f$
171 : * - `characteristic_speeds.get(1)`: degenerate speed \f$v_n\f$
172 : * - `characteristic_speeds.get(2)`: acoustic speed \f$\lambda_+\f$
173 : *
174 : * \see `grmhd::ValenciaDivClean::numerical_characteristics` for how these
175 : * speeds are paired with characteristic modes and characteristic projectors.
176 : */
177 :
178 : template <size_t ThermodynamicDim>
179 1 : void characteristic_speeds_hydro(
180 : gsl::not_null<tnsr::i<DataVector, 3>*> characteristic_speeds,
181 :
182 : /* primitive variables */
183 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
184 : const Scalar<DataVector>& rest_mass_density,
185 : const Scalar<DataVector>& specific_internal_energy,
186 : const Scalar<DataVector>& electron_fraction,
187 :
188 : /* other helpful quantities */
189 : const Scalar<DataVector>& lorentz_factor,
190 : const Scalar<DataVector>& specific_enthalpy,
191 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
192 : const tnsr::i<DataVector, 3>& unit_normal,
193 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
194 : equation_of_state);
195 :
196 : template <size_t ThermodynamicDim>
197 1 : tnsr::i<DataVector, 3> characteristic_speeds_hydro(
198 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
199 : const Scalar<DataVector>& rest_mass_density,
200 : const Scalar<DataVector>& specific_internal_energy,
201 : const Scalar<DataVector>& electron_fraction,
202 : const Scalar<DataVector>& lorentz_factor,
203 : const Scalar<DataVector>& specific_enthalpy,
204 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
205 : const tnsr::i<DataVector, 3>& unit_normal,
206 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
207 : equation_of_state);
208 : /// @}
209 :
210 : /**
211 : * \brief Compute the characteristic matrix (flux Jacobian) for relativistic
212 : * hydrodynamics + composition ($Y_e$), in a given direction.
213 : *
214 : * The characteristic matrix is the Jacobian of the flux along the unit normal
215 : * $n_i$ with respect to the conserved variables $U_n$,
216 : * \begin{equation}
217 : * A_{m}{}^{n} = \frac{\partial (n_i F^i_m)}{\partial U_n},
218 : * \end{equation}
219 : * evaluated with unit lapse and zero shift. With lapse $\alpha$ and shift
220 : * $\beta^i$, the physical characteristic speeds are $\alpha \lambda - \beta_n$,
221 : * where $\lambda$ are the eigenvalues of this matrix and
222 : * $\beta_n = n_i \beta^i$.
223 : *
224 : * The index $m$ labels the rows of the matrix (the flux components), and the
225 : * index $n$ labels its columns (the conserved variables). Both follow the order
226 : * $(D, S_i, \tau, D Y_e)$, which differs from the order of the system's
227 : * `variables_tag`.
228 : *
229 : * \note The normal must be unit with respect to the spatial metric.
230 : *
231 : * \note With a 2D equation of state, the pressure must be positive, since
232 : * $\kappa = \partial p / \partial \epsilon |_\rho$ is recovered from
233 : * $\kappa p / \rho^2$, the only form the 2D interface provides.
234 : *
235 : * \note The specific enthalpy must be consistent with the equation of state,
236 : * $h = 1 + \epsilon + p / \rho$. The matrix uses both $h$ and
237 : * $h - 1 = \epsilon + p / \rho$, which is computed from the equation of state
238 : * to avoid cancellation when $h \approx 1$.
239 : *
240 : * \see `grmhd::ValenciaDivClean::numerical_characteristics` for the
241 : * conventions relating this matrix to characteristic speeds, modes, and
242 : * projectors.
243 : *
244 : * \see `find_real_eigensystem` for the index conventions of the matrix.
245 : */
246 : template <size_t ThermodynamicDim>
247 1 : void flux_jacobian_hydro(
248 : gsl::not_null<tnsr::iJ<DataVector, 6>*> characteristic_matrix,
249 :
250 : /* primitive variables */
251 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
252 : const Scalar<DataVector>& rest_mass_density,
253 : const Scalar<DataVector>& specific_internal_energy,
254 : const Scalar<DataVector>& electron_fraction,
255 :
256 : /* other helpful quantities */
257 : const Scalar<DataVector>& lorentz_factor,
258 : const Scalar<DataVector>& specific_enthalpy,
259 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
260 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
261 : const tnsr::i<DataVector, 3>& unit_normal,
262 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
263 : equation_of_state);
264 :
265 : /**
266 : * \brief Compute numerical characteristic speeds (eigenvalues), characteristic
267 : * projectors (left eigenvectors), and characteristic modes (right eigenvectors)
268 : * in a given direction, from the primitive variables.
269 : *
270 : * \note Currently, this function only supports building the eigensystem for
271 : * relativistic hydrodynamics (no magnetic field) with composition dependence
272 : * ($Y_e$).
273 : *
274 : * The characteristic speeds $\lambda_{(i)}$ are the eigenvalues of the
275 : * characteristic matrix in a given direction, $A_{m}{}^{n}$. The characteristic
276 : * projectors $L^{im}$ and modes $R_{im}$ are its left and right eigenvectors,
277 : * with the index conventions of `find_real_eigensystem`. We provide the
278 : * characteristics in the following types:
279 : * - Eigenvalues $\lambda_{(i)}$: `tnsr::i<DataVector, 6> characteristic_speeds`
280 : * - Right eigenvectors $R_{in}$: `tnsr::ij<DataVector, 6> characteristic_modes`
281 : * - Left eigenvectors $L^{im}$: `tnsr::IJ<DataVector, 6>
282 : * characteristic_projectors`
283 : *
284 : * We refer to the left eigenvectors as "projectors" because they are used to
285 : * project a state into the characteristic basis. Similarly, we refer to the
286 : * right eigenvectors as "modes" because they are used to reconstruct the state.
287 : * For example, consider a state vector $U_{m}$ that can be represented in the
288 : * characteristic basis as
289 : * \begin{equation}
290 : * U_{m} = \sum_i w^{i} R_{im}.
291 : * \end{equation}
292 : * Then, the characteristic fields $w^{i}$ can be obtained by projecting the
293 : * state vector $U_{m}$ into the characteristic basis using the left
294 : * eigenvectors:
295 : * \begin{equation}
296 : * w^{i} = \sum_m L^{im} U_{m}.
297 : * \end{equation}
298 : *
299 : * The eigenvalues are in increasing order. The eigenvectors are normalized so
300 : * that $\sum_m L^{im} R_{jm} = \delta^{i}_{j}$. The right eigenvectors are not
301 : * unit vectors.
302 : *
303 : * \see `find_real_eigensystem` for the index conventions and for how this
304 : * eigensystem is computed at each grid point.
305 : */
306 : template <size_t ThermodynamicDim>
307 1 : void numerical_characteristics(
308 : gsl::not_null<tnsr::i<DataVector, 6>*> characteristic_speeds,
309 : gsl::not_null<tnsr::ij<DataVector, 6>*> characteristic_modes,
310 : gsl::not_null<tnsr::IJ<DataVector, 6>*> characteristic_projectors,
311 :
312 : /* primitive variables */
313 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
314 : const Scalar<DataVector>& rest_mass_density,
315 : const Scalar<DataVector>& specific_internal_energy,
316 : const Scalar<DataVector>& electron_fraction,
317 :
318 : /* other helpful quantities */
319 : const Scalar<DataVector>& lorentz_factor,
320 : const Scalar<DataVector>& specific_enthalpy,
321 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
322 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
323 : const tnsr::i<DataVector, 3>& unit_normal,
324 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
325 : equation_of_state);
326 :
327 1 : namespace Tags {
328 : /// \brief Compute the characteristic speeds for the Valencia formulation of
329 : /// GRMHD with divergence cleaning.
330 : ///
331 : /// \details see grmhd::ValenciaDivClean::characteristic_speeds
332 1 : struct CharacteristicSpeedsCompute : Tags::CharacteristicSpeeds,
333 : db::ComputeTag {
334 0 : using base = Tags::CharacteristicSpeeds;
335 0 : using argument_tags =
336 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
337 : hydro::Tags::ElectronFraction<DataVector>,
338 : hydro::Tags::SpecificInternalEnergy<DataVector>,
339 : hydro::Tags::SpecificEnthalpy<DataVector>,
340 : hydro::Tags::SpatialVelocity<DataVector, 3>,
341 : hydro::Tags::LorentzFactor<DataVector>,
342 : hydro::Tags::MagneticField<DataVector, 3>,
343 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
344 : gr::Tags::SpatialMetric<DataVector, 3>,
345 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<3>>,
346 : hydro::Tags::GrmhdEquationOfState>;
347 :
348 0 : using volume_tags = tmpl::list<hydro::Tags::GrmhdEquationOfState>;
349 :
350 0 : using return_type = std::array<DataVector, 9>;
351 :
352 : template <size_t ThermodynamicDim>
353 0 : void function(gsl::not_null<return_type*> result,
354 : const Scalar<DataVector>& rest_mass_density,
355 : const Scalar<DataVector>& /* electron_fraction */,
356 : const Scalar<DataVector>& specific_internal_energy,
357 : const Scalar<DataVector>& specific_enthalpy,
358 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
359 : const Scalar<DataVector>& lorentz_factor,
360 : const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
361 : const Scalar<DataVector>& lapse,
362 : const tnsr::I<DataVector, 3>& shift,
363 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
364 : const tnsr::i<DataVector, 3>& unit_normal,
365 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
366 : equation_of_state);
367 : };
368 :
369 0 : struct LargestCharacteristicSpeed : db::SimpleTag {
370 0 : using type = double;
371 : };
372 :
373 0 : struct ComputeLargestCharacteristicSpeed : db::ComputeTag,
374 : LargestCharacteristicSpeed {
375 0 : using argument_tags =
376 : tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
377 : gr::Tags::SpatialMetric<DataVector, 3>>;
378 0 : using return_type = double;
379 0 : using base = LargestCharacteristicSpeed;
380 0 : static void function(gsl::not_null<double*> speed,
381 : const Scalar<DataVector>& lapse,
382 : const tnsr::I<DataVector, 3>& shift,
383 : const tnsr::ii<DataVector, 3>& spatial_metric) {
384 : const auto shift_magnitude = magnitude(shift, spatial_metric);
385 : *speed = max(get(shift_magnitude) + get(lapse));
386 : }
387 : };
388 : } // namespace Tags
389 : } // namespace ValenciaDivClean
390 : } // namespace grmhd
|
