Hll.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <cstddef>
8 
12 #include "Options/Options.hpp"
13 #include "Utilities/Gsl.hpp"
15 #include "Utilities/TMPL.hpp"
16 
17 namespace dg {
18 namespace NumericalFluxes {
19 
20 /*!
21  * \ingroup DiscontinuousGalerkinGroup
22  * \ingroup NumericalFluxesGroup
23  * \brief Compute the HLL numerical flux.
24  *
25  * Let \f$U\f$ and \f$F^i(U)\f$ be the state vector of the system and its
26  * corresponding volume flux, respectively. Let \f$n_i\f$ be the unit normal to
27  * the interface. Defining \f$F_n(U) := n_i F^i(U)\f$, and denoting the
28  * corresponding projections of the numerical fluxes as \f${F_n}^*(U)\f$, the
29  * HLL flux \ref hll_ref "[1]" for each variable is
30  *
31  * \f{align*}
32  * {F_n}^*(U) = \frac{c_\text{max} F_n(U_\text{int}) -
33  * c_\text{min} F_n(U_\text{ext})}{c_\text{max} - c_\text{min}}
34  * - \frac{c_\text{min}c_\text{max}}{c_\text{max} - c_\text{min}}
35  * \left(U_\text{int} - U_\text{ext}\right)
36  * \f}
37  *
38  * where "int" and "ext" stand for interior and exterior, respectively.
39  * \f$c_\text{min}\f$ and \f$c_\text{max}\f$ are estimates on the minimum
40  * and maximum signal velocities bounding the interior-moving and
41  * exterior-moving wavespeeds that arise when solving the Riemann problem.
42  * Here we use the simple estimates \ref estimates_ref "[2]"
43  *
44  * \f{align*}
45  * c_\text{min} &= \text{min}\left( \lambda_1(U_\text{int}; n_\text{int}),
46  * \lambda_1(U_\text{ext}; n_\text{ext}), 0\right),\\
47  * c_\text{max} &= \text{max}\left( \lambda_N(U_\text{int}; n_\text{int}),
48  * \lambda_N(U_\text{ext}; n_\text{ext}), 0\right),
49  * \f}
50  *
51  * where \f$\lambda_1(U; n) \leq \lambda_2(U; n) \leq ... \leq
52  * \lambda_N(U; n)\f$ are the (ordered) characteristic speeds of the system.
53  * Note that the definitions above ensure that \f$c_\text{min} \leq 0\f$ and
54  * \f$c_\text{max} \geq 0\f$. Also, for either \f$c_\text{min} = 0\f$ or
55  * \f$c_\text{max} = 0\f$ (i.e. all characteristics move in the same direction)
56  * the HLL flux reduces to pure upwinding.
57  *
58  * \anchor hll_ref [1] A. Harten, P. D. Lax, B. van Leer, On Upstream
59  * Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws,
60  * SIAM Rev. [25 (1983) 35](https://doi.org/10.1137/1025002)
61  *
62  * \anchor estimates_ref [2] S. F. Davis, Simplified Second-Order Godunov-Type
63  * Methods, SIAM J. Sci. Stat. Comput.
64  * [9 (1988) 445](https://doi.org/10.1137/0909030)
65  */
66 template <typename System>
67 struct Hll {
68  private:
69  using char_speeds_tag = typename System::char_speeds_tag;
70  using variables_tag = typename System::variables_tag;
71 
72  public:
73  /// The minimum signal velocity bounding
74  /// the wavespeeds on one side of the interface.
76  using type = Scalar<DataVector>;
77  static std::string name() noexcept { return "MinSignalSpeed"; }
78  };
79 
80  /// The maximum signal velocity bounding
81  /// the wavespeeds on one side of the interface.
83  using type = Scalar<DataVector>;
84  static std::string name() noexcept { return "MaxSignalSpeed"; }
85  };
86 
87  using package_tags = tmpl::append<
90  tmpl::list<MinSignalSpeed, MaxSignalSpeed>>;
91 
92  using argument_tags =
93  tmpl::push_back<tmpl::append<db::split_tag<db::add_tag_prefix<
94  ::Tags::NormalDotFlux, variables_tag>>,
95  db::split_tag<variables_tag>>,
96  char_speeds_tag>;
97 
98  private:
99  template <typename VariablesTagList, typename NormalDoFluxTagList>
100  struct package_data_helper;
101 
102  template <typename... VariablesTags, typename... NormalDotFluxTags>
103  struct package_data_helper<tmpl::list<VariablesTags...>,
104  tmpl::list<NormalDotFluxTags...>> {
105  static void function(
106  const gsl::not_null<Variables<package_tags>*> packaged_data,
107  const db::item_type<NormalDotFluxTags>&... n_dot_f_to_package,
108  const db::item_type<VariablesTags>&... u_to_package,
109  const db::item_type<char_speeds_tag>& characteristic_speeds) noexcept {
110  expand_pack((get<VariablesTags>(*packaged_data) = u_to_package)...);
111  expand_pack(
112  (get<NormalDotFluxTags>(*packaged_data) = n_dot_f_to_package)...);
113 
114  get<MinSignalSpeed>(*packaged_data) = make_with_value<Scalar<DataVector>>(
115  characteristic_speeds[0],
117  get<MaxSignalSpeed>(*packaged_data) = make_with_value<Scalar<DataVector>>(
118  characteristic_speeds[0],
120 
121  // This finds the min and max characteristic speeds at each grid point,
122  // which are used as estimates of the min and max signal speeds.
123  for (size_t s = 0; s < characteristic_speeds[0].size(); ++s) {
124  // This ensures that local_min_signal_speed <= 0.0
125  const double local_min_signal_speed = (*std::min_element(
126  characteristic_speeds.begin(), characteristic_speeds.end(),
127  [&s](const auto& a, const auto& b) { return a[s] < b[s]; }))[s];
128  get(get<MinSignalSpeed>(*packaged_data))[s] =
129  std::min(local_min_signal_speed, 0.0);
130 
131  // Likewise, local_max_signal_speed >= 0.0
132  const double local_max_signal_speed = (*std::max_element(
133  characteristic_speeds.begin(), characteristic_speeds.end(),
134  [&s](const auto& a, const auto& b) { return a[s] < b[s]; }))[s];
135  get(get<MaxSignalSpeed>(*packaged_data))[s] =
136  std::max(local_max_signal_speed, 0.0);
137  }
138  }
139  };
140 
141  template <typename NormalDotNumericalFluxTagList, typename VariablesTagList,
142  typename NormalDotFluxTagList>
143  struct call_operator_helper;
144  template <typename... NormalDotNumericalFluxTags, typename... VariablesTags,
145  typename... NormalDotFluxTags>
146  struct call_operator_helper<tmpl::list<NormalDotNumericalFluxTags...>,
147  tmpl::list<VariablesTags...>,
148  tmpl::list<NormalDotFluxTags...>> {
149  static void function(
150  const gsl::not_null<
151  db::item_type<NormalDotNumericalFluxTags>*>... n_dot_numerical_f,
152  const db::item_type<NormalDotFluxTags>&... n_dot_f_interior,
153  const db::item_type<VariablesTags>&... u_interior,
154  const db::item_type<MinSignalSpeed>& min_signal_speed_interior,
155  const db::item_type<MaxSignalSpeed>& max_signal_speed_interior,
156  const db::item_type<NormalDotFluxTags>&... minus_n_dot_f_exterior,
157  const db::item_type<VariablesTags>&... u_exterior,
158  const db::item_type<MinSignalSpeed>& min_signal_speed_exterior,
160  max_signal_speed_exterior) noexcept {
161  auto min_signal_speed = make_with_value<db::item_type<MinSignalSpeed>>(
162  min_signal_speed_interior,
164  auto max_signal_speed = make_with_value<db::item_type<MaxSignalSpeed>>(
165  max_signal_speed_interior,
167  for (size_t s = 0; s < min_signal_speed.begin()->size(); ++s) {
168  get(min_signal_speed)[s] = std::min(get(min_signal_speed_interior)[s],
169  get(min_signal_speed_exterior)[s]);
170  get(max_signal_speed)[s] = std::max(get(max_signal_speed_interior)[s],
171  get(max_signal_speed_exterior)[s]);
172  }
173  const DataVector one_over_cp_minus_cm =
174  1.0 / (get(max_signal_speed) - get(min_signal_speed));
175  const auto assemble_numerical_flux =
176  [&min_signal_speed, &max_signal_speed, &one_over_cp_minus_cm ](
177  const auto n_dot_num_f, const auto& n_dot_f_in, const auto& u_in,
178  const auto& minus_n_dot_f_ex, const auto& u_ex) noexcept {
179  for (size_t i = 0; i < n_dot_num_f->size(); ++i) {
180  (*n_dot_num_f)[i] = one_over_cp_minus_cm *
181  (get(max_signal_speed) * n_dot_f_in[i] +
182  get(min_signal_speed) * minus_n_dot_f_ex[i] -
183  get(max_signal_speed) * get(min_signal_speed) *
184  (u_in[i] - u_ex[i]));
185  }
186  return nullptr;
187  };
188  expand_pack(assemble_numerical_flux(n_dot_numerical_f, n_dot_f_interior,
189  u_interior, minus_n_dot_f_exterior,
190  u_exterior)...);
191  }
192  };
193 
194  public:
195  using options = tmpl::list<>;
196  static constexpr OptionString help = {"Computes the HLL numerical flux."};
197 
198  // clang-tidy: google-runtime-references
199  void pup(PUP::er& /*p*/) noexcept {} // NOLINT
200 
201  template <class... Args>
202  void package_data(const Args&... args) const noexcept {
203  package_data_helper<
204  db::split_tag<variables_tag>,
205  db::split_tag<db::add_tag_prefix<::Tags::NormalDotFlux,
206  variables_tag>>>::function(args...);
207  }
208 
209  template <class... Args>
210  void operator()(const Args&... args) const noexcept {
211  call_operator_helper<
214  db::split_tag<variables_tag>,
215  db::split_tag<db::add_tag_prefix<::Tags::NormalDotFlux,
216  variables_tag>>>::function(args...);
217  }
218 };
219 
220 } // namespace NumericalFluxes
221 } // namespace dg
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Wrap Tag in Prefix<_, Args...>, also wrapping variables tags if Tag is a Tags::Variables.
Definition: DataBoxTag.hpp:533
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:76
T signaling_NaN(T... args)
Definition: InitializeElement.hpp:63
Compute the HLL numerical flux.
Definition: Hll.hpp:67
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:65
tmpl::conditional_t< tmpl::size< typename Subitems< TagList, Tag >::type >::value==0, tmpl::list< Tag >, typename Subitems< TagList, Tag >::type > split_tag
Split a tag into its subitems. Tag cannot be a base tag.
Definition: DataBoxTag.hpp:657
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
The maximum signal velocity bounding the wavespeeds on one side of the interface. ...
Definition: Hll.hpp:82
The minimum signal velocity bounding the wavespeeds on one side of the interface. ...
Definition: Hll.hpp:75
Defines classes for Tensor.
Stores a collection of function values.
Definition: DataVector.hpp:46
Wraps the template metaprogramming library used (brigand)
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
Defines functions and classes from the GSL.
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:545
Defines make_with_value.