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/Matrix.hpp"
11 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp" // IWYU pragma: keep
12 : #include "DataStructures/Tensor/Tensor.hpp" // for get
13 : #include "DataStructures/Tensor/TypeAliases.hpp" // IWYU pragma: keep
14 : #include "Domain/FaceNormal.hpp"
15 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
16 : #include "PointwiseFunctions/Hydro/Tags.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : /// \cond
21 : class DataVector;
22 : /// \endcond
23 :
24 : // IWYU pragma: no_forward_declare Tensor
25 :
26 : namespace NewtonianEuler {
27 :
28 : namespace detail {
29 : // Compute the flux Jacobian for the NewtonianEuler system
30 : //
31 : // The flux Jacobian is \f$n_i A^i = n_i \partial F^i / \partial U\f$.
32 : // The input `b_times_theta` is \f$b\theta = \kappa/\rho (v^2 - h) + c_s^2\f$.
33 : //
34 : // This is used for:
35 : // - testing the analytic characteristic transformation
36 : // - as input for the numerical characteristic transformation
37 : template <size_t Dim>
38 : Matrix flux_jacobian(const tnsr::I<double, Dim>& velocity,
39 : double kappa_over_density, double b_times_theta,
40 : double specific_enthalpy,
41 : const tnsr::i<double, Dim>& unit_normal);
42 : } // namespace detail
43 :
44 : /// @{
45 : /*!
46 : * \brief Compute the characteristic speeds of NewtonianEuler system
47 : *
48 : * The principal symbol of the system is diagonalized so that the elements of
49 : * the diagonal matrix are the characteristic speeds
50 : *
51 : * \f{align*}
52 : * \lambda_1 &= v_n - c_s,\\
53 : * \lambda_{i + 1} &= v_n,\\
54 : * \lambda_{\text{Dim} + 2} &= v_n + c_s,
55 : * \f}
56 : *
57 : * where \f$i = 1,...,\text{Dim}\f$,
58 : * \f$v_n = n_i v^i\f$ is the velocity projected onto the normal,
59 : * and \f$c_s\f$ is the sound speed.
60 : */
61 : template <size_t Dim>
62 1 : void characteristic_speeds(
63 : gsl::not_null<std::array<DataVector, Dim + 2>*> char_speeds,
64 : const tnsr::I<DataVector, Dim>& velocity,
65 : const Scalar<DataVector>& sound_speed,
66 : const tnsr::i<DataVector, Dim>& normal);
67 :
68 : template <size_t Dim>
69 1 : std::array<DataVector, Dim + 2> characteristic_speeds(
70 : const tnsr::I<DataVector, Dim>& velocity,
71 : const Scalar<DataVector>& sound_speed,
72 : const tnsr::i<DataVector, Dim>& normal);
73 : /// @}
74 :
75 : /// @{
76 : /*!
77 : * \brief Compute the transform matrices between the conserved variables and
78 : * the characteristic variables of the NewtonianEuler system.
79 : *
80 : * Let \f$u\f$ be the conserved (i.e., evolved) variables of the Newtonian Euler
81 : * system, and \f$w\f$ the characteristic variables of this system with respect
82 : * to a unit normal one form \f$n_i\f$. The function `left_eigenvectors`
83 : * computes the matrix \f$\Omega_{L}\f$ corresponding to the transform
84 : * \f$w = \Omega_{L} u\f$. The function `right_eigenvectors` computes the matrix
85 : * \f$\Omega_{R}\f$ corresponding to the inverse transform
86 : * \f$u = \Omega_{R} w\f$. Here the components of \f$u\f$ are ordered as
87 : * \f$u = \{\rho, \rho v_x, \rho v_y, \rho v_z, e\}\f$ in 3D, and the components
88 : * of \f$w\f$ are ordered by their corresponding eigenvalues
89 : * (i.e., characteristic speeds)
90 : * \f$\lambda = \{v_n - c_s, v_n, v_n, v_n, v_n + c_s\}\f$. In these
91 : * expressions, \f$\rho\f$ is the fluid mass density, \f$v_{x,y,z}\f$ are the
92 : * components of the fluid velocity, \f$e\f$ is the total energy density,
93 : * \f$v_n\f$ is the component of the velocity along the unit normal \f$n_i\f$,
94 : * and \f$c_s\f$ is the sound speed.
95 : *
96 : * For a short discussion of the characteristic transformation and the matrices
97 : * \f$\Omega_{L}\f$ and \f$\Omega_{R}\f$, see \cite Kulikovskii2000 Chapter 3.
98 : *
99 : * Here we briefly summarize the procedure. With \f$F^x(u)\f$ the Newtonian
100 : * Euler flux in direction \f$x\f$, then the flux Jacobian along \f$x\f$ is the
101 : * matrix \f$A_x = \partial F^x_{\beta}(u) / \partial u_{\alpha}\f$. The indices
102 : * \f$\alpha, \beta\f$ range over the different evolved fields. In higher
103 : * dimensions, the flux Jacobian along the unit normal \f$n_i\f$ is
104 : * \f$A = n_x A_x + n_y A_y + n_z A_z\f$.
105 : * This matrix can be diagonalized as \f$A = \Omega_{R} \Lambda \Omega_{L}\f$.
106 : * Here \f$\Lambda = \mathrm{diag}(v_n - c_s, v_n, v_n, v_n, v_n + c_s)\f$
107 : * is a diagonal matrix containing the characteristic speeds; \f$\Omega_{R}\f$
108 : * is a matrix whose columns are the right eigenvectors of \f$A\f$;
109 : * \f$\Omega_{L}\f$ is the inverse of \f$R\f$.
110 : */
111 : template <size_t Dim>
112 1 : Matrix right_eigenvectors(const tnsr::I<double, Dim>& velocity,
113 : const Scalar<double>& sound_speed_squared,
114 : const Scalar<double>& specific_enthalpy,
115 : const Scalar<double>& kappa_over_density,
116 : const tnsr::i<double, Dim>& unit_normal);
117 :
118 : template <size_t Dim>
119 1 : Matrix left_eigenvectors(const tnsr::I<double, Dim>& velocity,
120 : const Scalar<double>& sound_speed_squared,
121 : const Scalar<double>& specific_enthalpy,
122 : const Scalar<double>& kappa_over_density,
123 : const tnsr::i<double, Dim>& unit_normal);
124 : /// @}
125 :
126 : /*!
127 : * \brief Compute the transform matrices between the conserved variables and
128 : * the characteristic variables of the NewtonianEuler system.
129 : *
130 : * See `right_eigenvectors` and `left_eigenvectors` for more details.
131 : *
132 : * However, note that this function computes the transformation (i.e., the
133 : * eigenvectors of the flux Jacobian) numerically, instead of using the analytic
134 : * expressions. This is useful as a proof-of-concept for more complicated
135 : * systems where the analytic expressions may not be known.
136 : */
137 : template <size_t Dim>
138 1 : std::pair<DataVector, std::pair<Matrix, Matrix>> numerical_eigensystem(
139 : const tnsr::I<double, Dim>& velocity,
140 : const Scalar<double>& sound_speed_squared,
141 : const Scalar<double>& specific_enthalpy,
142 : const Scalar<double>& kappa_over_density,
143 : const tnsr::i<double, Dim>& unit_normal);
144 :
145 1 : namespace Tags {
146 :
147 : template <size_t Dim>
148 0 : struct CharacteristicSpeedsCompute : CharacteristicSpeeds<Dim>, db::ComputeTag {
149 0 : using base = CharacteristicSpeeds<Dim>;
150 0 : using argument_tags =
151 : tmpl::list<hydro::Tags::SpatialVelocity<DataVector, Dim>,
152 : SoundSpeed<DataVector>,
153 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>>;
154 :
155 0 : using return_type = std::array<DataVector, Dim + 2>;
156 :
157 0 : static constexpr void function(const gsl::not_null<return_type*> result,
158 : const tnsr::I<DataVector, Dim>& velocity,
159 : const Scalar<DataVector>& sound_speed,
160 : const tnsr::i<DataVector, Dim>& normal) {
161 : characteristic_speeds<Dim>(result, velocity, sound_speed, normal);
162 : }
163 : };
164 :
165 0 : struct LargestCharacteristicSpeed : db::SimpleTag {
166 0 : using type = double;
167 : };
168 :
169 : template <size_t Dim>
170 0 : struct ComputeLargestCharacteristicSpeed : LargestCharacteristicSpeed,
171 : db::ComputeTag {
172 0 : using argument_tags =
173 : tmpl::list<hydro::Tags::SpatialVelocity<DataVector, Dim>,
174 : Tags::SoundSpeed<DataVector>>;
175 0 : using return_type = double;
176 0 : using base = LargestCharacteristicSpeed;
177 0 : static void function(const gsl::not_null<double*> speed,
178 : const tnsr::I<DataVector, Dim>& velocity,
179 : const Scalar<DataVector>& sound_speed) {
180 : *speed = max(get(magnitude(velocity)) + get(sound_speed));
181 : }
182 : };
183 : } // namespace Tags
184 : } // namespace NewtonianEuler
|
