Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <memory>
7 : #include <optional>
8 :
9 : #include "DataStructures/DataBox/Prefixes.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Evolution/BoundaryCorrection.hpp"
12 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
13 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
14 : #include "Options/String.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
16 : #include "PointwiseFunctions/Hydro/Tags.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/Serialization/CharmPupable.hpp"
19 : #include "Utilities/TMPL.hpp"
20 :
21 : /// \cond
22 : class DataVector;
23 : namespace gsl {
24 : template <typename T>
25 : class not_null;
26 : } // namespace gsl
27 : namespace PUP {
28 : class er;
29 : } // namespace PUP
30 : /// \endcond
31 :
32 : namespace NewtonianEuler::BoundaryCorrections {
33 : /*!
34 : * \brief The HLLC (Harten-Lax-van Leer-Contact) Riemann solver for the
35 : * NewtonianEuler system
36 : *
37 : * Let \f$U\f$ be the evolved variable, \f$F^i\f$ the flux, \f$v^i\f$ the
38 : * spatial velocity, and \f$n_i\f$ be the outward directed unit normal to the
39 : * interface. Denoting \f$F := n_i F^i\f$ and \f$v:=n_iv^i\f$, the HLLC boundary
40 : * correction is \cite Toro1994
41 : *
42 : * \f{align*}
43 : * G_\text{HLLC} = \left\{\begin{array}{lcl}
44 : * F_\text{int} & \text{if} & 0 \leq \lambda_\text{min} \\
45 : * F_\text{int} + \lambda_\text{min}(U_\text{*int} - U_\text{int}) & \text{if} &
46 : * \lambda_\text{min} \leq 0 \leq \lambda_\text{*} \\
47 : * -F_\text{ext} + \lambda_\text{max}(U_\text{*ext} - U_\text{ext}) & \text{if}
48 : * & \lambda_\text{*} \leq 0 \leq \lambda_\text{max} \\
49 : * -F_\text{ext} & \text{if} & \lambda_\text{max} \leq 0
50 : * \end{array}\right\}
51 : * \f}
52 : *
53 : * where "int" and "ext" stand for interior and exterior.
54 : *
55 : * Intermediate ('star') states are given by
56 : *
57 : * \f{align*}
58 : * U_\text{*int} = \left(\frac{\lambda_\text{min} - v_\text{int}}
59 : * {\lambda_\text{min} - \lambda_*}\right)
60 : * \left[\begin{array}{c}
61 : * \displaystyle \rho_\text{int} \\
62 : * \displaystyle \rho_\text{int}[v_\text{int}^x + (\lambda_* - v_\text{int})
63 : * n_x^\text{int}] \\
64 : * \displaystyle \rho_\text{int}[v_\text{int}^y + (\lambda_* - v_\text{int})
65 : * n_y^\text{int}] \\
66 : * \displaystyle \rho_\text{int}[v_\text{int}^z + (\lambda_* - v_\text{int})
67 : * n_z^\text{int}] \\
68 : * \displaystyle E_\text{int} + p_\text{int} \frac{\lambda_* - v_\text{int}}
69 : * {\lambda_\text{min} - v_\text{int}} + \rho_\text{int}\lambda_*(\lambda_* -
70 : * v_\text{int}) \end{array}\right]
71 : * \f}
72 : *
73 : * and
74 : *
75 : * \f{align*}
76 : * U_\text{*ext} = \left(\frac{\lambda_\text{max} + v_\text{ext}}
77 : * {\lambda_\text{max} - \lambda_*}\right)
78 : * \left[\begin{array}{c}
79 : * \displaystyle \rho_\text{ext} \\
80 : * \displaystyle \rho_\text{ext}[-v_\text{ext}^x - (\lambda_* + v_\text{ext})
81 : * n_x^\text{ext}] \\
82 : * \displaystyle \rho_\text{ext}[-v_\text{ext}^y - (\lambda_* + v_\text{ext})
83 : * n_y^\text{ext}] \\
84 : * \displaystyle \rho_\text{ext}[-v_\text{ext}^z - (\lambda_* + v_\text{ext})
85 : * n_z^\text{ext}] \\
86 : * \displaystyle E_\text{ext} + p_\text{ext} \frac{\lambda_* + v_\text{ext}}
87 : * {\lambda_\text{max} + v_\text{ext}} + \rho_\text{ext}\lambda_*(\lambda_* +
88 : * v_\text{ext}) \end{array}\right].
89 : * \f}
90 : *
91 : * The contact wave speed \f$\lambda_*\f$ is \cite Toro2009
92 : *
93 : * \f{align*}
94 : * \lambda_* = \frac
95 : * { p_\text{ext} - p_\text{int} +
96 : * \rho_\text{int}v_\text{int}(\lambda_\text{min}-v_\text{int})
97 : * + \rho_\text{ext}v_\text{ext}(\lambda_\text{max}+v_\text{ext})}
98 : * { \rho_\text{int}(\lambda_\text{min}-v_\text{int}) -
99 : * \rho_\text{ext}(\lambda_\text{max} + v_\text{ext})}.
100 : * \f}
101 : *
102 : * \f$\lambda_\text{min}\f$ and \f$\lambda_\text{max}\f$ are estimated by
103 : * \cite Davis1988
104 : *
105 : * \f{align*}
106 : * \lambda_\text{min} &=
107 : * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}\right) \\
108 : * \lambda_\text{max} &=
109 : * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}\right)
110 : * \f}
111 : *
112 : * where \f$\lambda^{+}\f$ (\f$\lambda^{-}\f$) is the largest characteristic
113 : * speed in the outgoing (ingoing) direction for each domain.
114 : *
115 : * Note the minus signs in front of \f$\lambda^{\pm}_\text{ext}\f$, which is
116 : * because an outgoing speed w.r.t. the neighboring element is an ingoing speed
117 : * w.r.t. the local element, and vice versa. Similarly, the \f$F_{\text{ext}}\f$
118 : * term in \f$G_\text{HLLC}\f$ and the \f$v_\text{ext}\f$ term in
119 : * \f$U_\text{*ext}\f$ have a positive sign because the outward directed normal
120 : * of the neighboring element has the opposite sign, i.e.
121 : * \f$n_i^{\text{ext}}=-n_i^{\text{int}}\f$.
122 : *
123 : * For the NewtonianEuler system, \f$\lambda^\pm\f$ are given as
124 : *
125 : * \f{align*}
126 : * \lambda^\pm = v \pm c_s
127 : * \f}
128 : *
129 : * where \f$c_s\f$ is the sound speed.
130 : *
131 : * \note
132 : * - In the strong form the `dg_boundary_terms` function returns \f$G -
133 : * F_\text{int}\f$
134 : * - In the implementation, we use
135 : *
136 : * \f{align*}
137 : * G_\text{HLLC} = \left\{\begin{array}{lcl}
138 : * F_\text{int} + \lambda_\text{min}(U_\text{*int} - U_\text{int}) & \text{if}
139 : * & 0 \leq \lambda_\text{*} \\
140 : * -F_\text{ext} + \lambda_\text{max}(U_\text{*ext} - U_\text{ext}) &
141 : * \text{if} & \lambda_\text{*} \leq 0 \\
142 : * \end{array}\right\},
143 : * \f}
144 : *
145 : * with
146 : *
147 : * \f{align*}
148 : * \lambda_\text{min} &=
149 : * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}, 0\right) \\
150 : * \lambda_\text{max} &=
151 : * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}, 0\right).
152 : * \f}
153 : *
154 : * Provided that \f$\lambda_*\f$ falls in the correct range i.e
155 : *
156 : * \f{align*}
157 : * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}\right)
158 : * < \lambda_* <
159 : * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}\right),
160 : * \f}
161 : *
162 : * this prescription recovers the original HLLC boundary correction for all
163 : * four cases. For either \f$\lambda_\text{min} = 0\f$ or
164 : * \f$\lambda_\text{max} = 0\f$ (i.e. all characteristics move in the same
165 : * direction), boundary correction reduces to pure upwinding.
166 : * - Some references use \f$S\f$ instead of \f$\lambda\f$ for the
167 : * signal/characteristic speeds.
168 : */
169 : template <size_t Dim>
170 1 : class Hllc final : public evolution::BoundaryCorrection {
171 : private:
172 0 : struct InterfaceUnitNormal : db::SimpleTag {
173 0 : using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
174 : };
175 0 : struct NormalDotVelocity : db::SimpleTag {
176 0 : using type = Scalar<DataVector>;
177 : };
178 :
179 : public:
180 0 : struct LargestOutgoingCharSpeed : db::SimpleTag {
181 0 : using type = Scalar<DataVector>;
182 : };
183 0 : struct LargestIngoingCharSpeed : db::SimpleTag {
184 0 : using type = Scalar<DataVector>;
185 : };
186 :
187 0 : using options = tmpl::list<>;
188 0 : static constexpr Options::String help = {
189 : "Computes the HLLC boundary correction term for the "
190 : "Newtonian Euler/hydrodynamics system."};
191 :
192 0 : Hllc() = default;
193 0 : Hllc(const Hllc&) = default;
194 0 : Hllc& operator=(const Hllc&) = default;
195 0 : Hllc(Hllc&&) = default;
196 0 : Hllc& operator=(Hllc&&) = default;
197 0 : ~Hllc() override = default;
198 :
199 : /// \cond
200 : explicit Hllc(CkMigrateMessage* msg);
201 : using PUP::able::register_constructor;
202 : WRAPPED_PUPable_decl_template(Hllc); // NOLINT
203 : /// \endcond
204 0 : void pup(PUP::er& p) override; // NOLINT
205 :
206 0 : std::unique_ptr<BoundaryCorrection> get_clone() const override;
207 :
208 0 : using dg_package_field_tags =
209 : tmpl::list<Tags::MassDensityCons, Tags::MomentumDensity<Dim>,
210 : Tags::EnergyDensity, hydro::Tags::Pressure<DataVector>,
211 : ::Tags::NormalDotFlux<Tags::MassDensityCons>,
212 : ::Tags::NormalDotFlux<Tags::MomentumDensity<Dim>>,
213 : ::Tags::NormalDotFlux<Tags::EnergyDensity>,
214 : InterfaceUnitNormal, NormalDotVelocity,
215 : LargestOutgoingCharSpeed, LargestIngoingCharSpeed>;
216 0 : using dg_package_data_temporary_tags = tmpl::list<>;
217 0 : using dg_package_data_primitive_tags =
218 : tmpl::list<hydro::Tags::SpatialVelocity<DataVector, Dim>,
219 : hydro::Tags::SpecificInternalEnergy<DataVector>>;
220 0 : using dg_package_data_volume_tags =
221 : tmpl::list<hydro::Tags::EquationOfState<false, 2>>;
222 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
223 :
224 0 : double dg_package_data(
225 : gsl::not_null<Scalar<DataVector>*> packaged_mass_density,
226 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
227 : packaged_momentum_density,
228 : gsl::not_null<Scalar<DataVector>*> packaged_energy_density,
229 : gsl::not_null<Scalar<DataVector>*> packaged_pressure,
230 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_mass_density,
231 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
232 : packaged_normal_dot_flux_momentum_density,
233 : gsl::not_null<Scalar<DataVector>*>
234 : packaged_normal_dot_flux_energy_density,
235 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
236 : packaged_interface_unit_normal,
237 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_velocity,
238 : gsl::not_null<Scalar<DataVector>*> packaged_largest_outgoing_char_speed,
239 : gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_char_speed,
240 :
241 : const Scalar<DataVector>& mass_density,
242 : const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density,
243 : const Scalar<DataVector>& energy_density,
244 :
245 : const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_mass_density,
246 : const tnsr::IJ<DataVector, Dim, Frame::Inertial>& flux_momentum_density,
247 : const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_energy_density,
248 :
249 : const tnsr::I<DataVector, Dim, Frame::Inertial>& velocity,
250 : const Scalar<DataVector>& specific_internal_energy,
251 :
252 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
253 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
254 : /*mesh_velocity*/,
255 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
256 : const EquationsOfState::EquationOfState<false, 2>& equation_of_state)
257 : const;
258 :
259 0 : void dg_boundary_terms(
260 : gsl::not_null<Scalar<DataVector>*> boundary_correction_mass_density,
261 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
262 : boundary_correction_momentum_density,
263 : gsl::not_null<Scalar<DataVector>*> boundary_correction_energy_density,
264 : const Scalar<DataVector>& mass_density_int,
265 : const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_int,
266 : const Scalar<DataVector>& energy_density_int,
267 : const Scalar<DataVector>& pressure_int,
268 : const Scalar<DataVector>& normal_dot_flux_mass_density_int,
269 : const tnsr::I<DataVector, Dim, Frame::Inertial>&
270 : normal_dot_flux_momentum_density_int,
271 : const Scalar<DataVector>& normal_dot_flux_energy_density_int,
272 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
273 : interface_unit_normal_int,
274 : const Scalar<DataVector>& normal_dot_velocity_int,
275 : const Scalar<DataVector>& largest_outgoing_char_speed_int,
276 : const Scalar<DataVector>& largest_ingoing_char_speed_int,
277 : const Scalar<DataVector>& mass_density_ext,
278 : const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_ext,
279 : const Scalar<DataVector>& energy_density_ext,
280 : const Scalar<DataVector>& pressure_ext,
281 : const Scalar<DataVector>& normal_dot_flux_mass_density_ext,
282 : const tnsr::I<DataVector, Dim, Frame::Inertial>&
283 : normal_dot_flux_momentum_density_ext,
284 : const Scalar<DataVector>& normal_dot_flux_energy_density_ext,
285 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
286 : interface_unit_normal_ext,
287 : const Scalar<DataVector>& normal_dot_velocity_ext,
288 : const Scalar<DataVector>& largest_outgoing_char_speed_ext,
289 : const Scalar<DataVector>& largest_ingoing_char_speed_ext,
290 : dg::Formulation dg_formulation) const;
291 : };
292 : } // namespace NewtonianEuler::BoundaryCorrections
|
