Hllc.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <memory>
7 #include <optional>
8 
11 #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp"
12 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
13 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
14 #include "Options/Options.hpp"
16 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
17 #include "PointwiseFunctions/Hydro/Tags.hpp"
18 #include "Utilities/Gsl.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 
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 class Hllc final : public BoundaryCorrection<Dim> {
171  private:
172  struct InterfaceUnitNormal : db::SimpleTag {
173  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
174  };
175  struct NormalDotVelocity : db::SimpleTag {
176  using type = Scalar<DataVector>;
177  };
178 
179  public:
181  using type = Scalar<DataVector>;
182  };
184  using type = Scalar<DataVector>;
185  };
186 
187  using options = tmpl::list<>;
188  static constexpr Options::String help = {
189  "Computes the HLLC boundary correction term for the "
190  "Newtonian Euler/hydrodynamics system."};
191 
192  Hllc() = default;
193  Hllc(const Hllc&) = default;
194  Hllc& operator=(const Hllc&) = default;
195  Hllc(Hllc&&) = default;
196  Hllc& operator=(Hllc&&) = default;
197  ~Hllc() override = default;
198 
199  /// \cond
200  explicit Hllc(CkMigrateMessage* msg) noexcept;
201  using PUP::able::register_constructor;
203  /// \endcond
204  void pup(PUP::er& p) override; // NOLINT
205 
206  std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const noexcept override;
207 
208  using dg_package_field_tags = tmpl::list<
209  Tags::MassDensityCons, Tags::MomentumDensity<Dim>, Tags::EnergyDensity,
210  Tags::Pressure<DataVector>, ::Tags::NormalDotFlux<Tags::MassDensityCons>,
211  ::Tags::NormalDotFlux<Tags::MomentumDensity<Dim>>,
212  ::Tags::NormalDotFlux<Tags::EnergyDensity>, InterfaceUnitNormal,
213  NormalDotVelocity, LargestOutgoingCharSpeed, LargestIngoingCharSpeed>;
214  using dg_package_data_temporary_tags = tmpl::list<>;
215  using dg_package_data_primitive_tags =
216  tmpl::list<NewtonianEuler::Tags::Velocity<DataVector, Dim>,
217  NewtonianEuler::Tags::SpecificInternalEnergy<DataVector>>;
218  using dg_package_data_volume_tags =
219  tmpl::list<hydro::Tags::EquationOfStateBase>;
220 
221  template <size_t ThermodynamicDim>
222  double dg_package_data(
223  gsl::not_null<Scalar<DataVector>*> packaged_mass_density,
224  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
225  packaged_momentum_density,
226  gsl::not_null<Scalar<DataVector>*> packaged_energy_density,
227  gsl::not_null<Scalar<DataVector>*> packaged_pressure,
228  gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_mass_density,
229  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
230  packaged_normal_dot_flux_momentum_density,
231  gsl::not_null<Scalar<DataVector>*>
232  packaged_normal_dot_flux_energy_density,
233  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
234  packaged_interface_unit_normal,
235  gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_velocity,
236  gsl::not_null<Scalar<DataVector>*> packaged_largest_outgoing_char_speed,
237  gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_char_speed,
238 
239  const Scalar<DataVector>& mass_density,
240  const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density,
241  const Scalar<DataVector>& energy_density,
242 
243  const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_mass_density,
244  const tnsr::IJ<DataVector, Dim, Frame::Inertial>& flux_momentum_density,
245  const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_energy_density,
246 
247  const tnsr::I<DataVector, Dim, Frame::Inertial>& velocity,
248  const Scalar<DataVector>& specific_internal_energy,
249 
250  const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
251  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
252  /*mesh_velocity*/,
253  const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
254  const EquationsOfState::EquationOfState<false, ThermodynamicDim>&
255  equation_of_state) const noexcept;
256 
257  void dg_boundary_terms(
258  gsl::not_null<Scalar<DataVector>*> boundary_correction_mass_density,
259  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
260  boundary_correction_momentum_density,
261  gsl::not_null<Scalar<DataVector>*> boundary_correction_energy_density,
262  const Scalar<DataVector>& mass_density_int,
263  const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_int,
264  const Scalar<DataVector>& energy_density_int,
265  const Scalar<DataVector>& pressure_int,
266  const Scalar<DataVector>& normal_dot_flux_mass_density_int,
267  const tnsr::I<DataVector, Dim, Frame::Inertial>&
268  normal_dot_flux_momentum_density_int,
269  const Scalar<DataVector>& normal_dot_flux_energy_density_int,
270  const tnsr::i<DataVector, Dim, Frame::Inertial>&
271  interface_unit_normal_int,
272  const Scalar<DataVector>& normal_dot_velocity_int,
273  const Scalar<DataVector>& largest_outgoing_char_speed_int,
274  const Scalar<DataVector>& largest_ingoing_char_speed_int,
275  const Scalar<DataVector>& mass_density_ext,
276  const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_ext,
277  const Scalar<DataVector>& energy_density_ext,
278  const Scalar<DataVector>& pressure_ext,
279  const Scalar<DataVector>& normal_dot_flux_mass_density_ext,
280  const tnsr::I<DataVector, Dim, Frame::Inertial>&
281  normal_dot_flux_momentum_density_ext,
282  const Scalar<DataVector>& normal_dot_flux_energy_density_ext,
283  const tnsr::i<DataVector, Dim, Frame::Inertial>&
284  interface_unit_normal_ext,
285  const Scalar<DataVector>& normal_dot_velocity_ext,
286  const Scalar<DataVector>& largest_outgoing_char_speed_ext,
287  const Scalar<DataVector>& largest_ingoing_char_speed_ext,
288  dg::Formulation dg_formulation) const noexcept;
289 };
290 } // namespace NewtonianEuler::BoundaryCorrections
EquationsOfState
Contains all equations of state, including base class.
Definition: DarkEnergyFluid.hpp:26
CharmPupable.hpp
Options.hpp
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
NewtonianEuler
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8
NewtonianEuler::BoundaryCorrections::Hllc::LargestOutgoingCharSpeed
Definition: Hllc.hpp:180
NewtonianEuler::BoundaryCorrections
Boundary corrections/numerical fluxes.
Definition: BoundaryCorrection.hpp:13
dg
Functionality related to discontinuous Galerkin schemes.
Definition: BoundaryData.hpp:17
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
NewtonianEuler::BoundaryCorrections::Hllc
The HLLC (Harten-Lax-van Leer-Contact) Riemann solver for the NewtonianEuler system.
Definition: Hllc.hpp:170
hydro
Items related to hydrodynamic systems.
Definition: SmoothFlow.hpp:19
memory
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
NewtonianEuler::BoundaryCorrections::BoundaryCorrection
The base class used to create boundary corrections from input files and store them in the global cach...
Definition: BoundaryCorrection.hpp:28
Gsl.hpp
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
optional
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
std::unique_ptr
Prefixes.hpp
NewtonianEuler::BoundaryCorrections::Hllc::LargestIngoingCharSpeed
Definition: Hllc.hpp:183
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:11
TMPL.hpp