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