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 
9 #include "DataStructures/DataBox/PrefixHelpers.hpp"
10 #include "DataStructures/DataBox/Prefixes.hpp"
11 #include "DataStructures/DataBox/Tag.hpp"
12 #include "DataStructures/Tensor/Tensor.hpp"
13 #include "DataStructures/VariablesTag.hpp"
14 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
15 #include "Options/Options.hpp"
16 #include "Utilities/Gsl.hpp"
17 #include "Utilities/MakeWithValue.hpp"
18 #include "Utilities/ProtocolHelpers.hpp"
19 #include "Utilities/TMPL.hpp"
20 
21 namespace dg {
22 namespace NumericalFluxes {
23 
24 /*!
25  * \ingroup DiscontinuousGalerkinGroup
26  * \ingroup NumericalFluxesGroup
27  * \brief Compute the HLL numerical flux.
28  *
29  * Let \f$U\f$ be the state vector of the system and \f$F^i\f$ the corresponding
30  * volume fluxes. Let \f$n_i\f$ be the unit normal to
31  * the interface. Denoting \f$F := n_i F^i\f$, the HLL flux is \cite Harten1983
32  *
33  * \f{align*}
34  * G_\text{HLL} = \frac{S_\text{max} F_\text{int} -
35  * S_\text{min} F_\text{ext}}{S_\text{max} - S_\text{min}}
36  * - \frac{S_\text{min}S_\text{max}}{S_\text{max} - S_\text{min}}
37  * \left(U_\text{int} - U_\text{ext}\right)
38  * \f}
39  *
40  * where "int" and "ext" stand for interior and exterior, respectively.
41  * \f$S_\text{min}\f$ and \f$S_\text{max}\f$ are estimates on the minimum
42  * and maximum signal velocities bounding the ingoing and outgoing wavespeeds
43  * that arise when solving the Riemann problem. Here we use the simple
44  * estimates \cite Davis1988
45  *
46  * \f{align*}
47  * S_\text{min} &=
48  * \text{min}\left(\{\lambda_\text{int}\},\{\lambda_\text{ext}\}, 0\right)\\
49  * S_\text{max} &=
50  * \text{max}\left(\{\lambda_\text{int}\},\{\lambda_\text{ext}\}, 0\right),
51  * \f}
52  *
53  * where \f$\{\lambda\}\f$ is the set of all the characteristic speeds along a
54  * given normal. Note that for either \f$S_\text{min} = 0\f$ or
55  * \f$S_\text{max} = 0\f$ (i.e. all characteristics move in the same direction)
56  * the HLL flux reduces to pure upwinding.
57  */
58 template <typename System>
59 struct Hll : tt::ConformsTo<dg::protocols::NumericalFlux> {
60  private:
61  using char_speeds_tag = typename System::char_speeds_tag;
62  using variables_tag = typename System::variables_tag;
63 
64  public:
65  /// Estimate for the largest ingoing signal speed
66  struct LargestIngoingSpeed : db::SimpleTag {
67  using type = Scalar<DataVector>;
68  };
69 
70  /// Estimate for the largest outgoing signal speed
71  struct LargestOutgoingSpeed : db::SimpleTag {
72  using type = Scalar<DataVector>;
73  };
74 
75  using variables_tags = typename System::variables_tag::tags_list;
76 
77  using package_field_tags =
78  tmpl::append<db::wrap_tags_in<::Tags::NormalDotFlux, variables_tags>,
79  variables_tags,
80  tmpl::list<LargestIngoingSpeed, LargestOutgoingSpeed>>;
81  using package_extra_tags = tmpl::list<>;
82 
83  using argument_tags = tmpl::push_back<
84  tmpl::append<db::wrap_tags_in<::Tags::NormalDotFlux, variables_tags>,
85  variables_tags>,
86  char_speeds_tag>;
87 
88  private:
89  template <typename VariablesTagList, typename NormalDotFluxTagList>
90  struct package_data_helper;
91 
92  template <typename... VariablesTags, typename... NormalDotFluxTags>
93  struct package_data_helper<tmpl::list<VariablesTags...>,
94  tmpl::list<NormalDotFluxTags...>> {
95  static void function(
96  const gsl::not_null<
97  typename NormalDotFluxTags::type*>... packaged_n_dot_f,
98  const gsl::not_null<typename VariablesTags::type*>... packaged_u,
99  const gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_speed,
100  const gsl::not_null<Scalar<DataVector>*>
101  packaged_largest_outgoing_speed,
102  const typename NormalDotFluxTags::type&... n_dot_f_to_package,
103  const typename VariablesTags::type&... u_to_package,
104  const typename char_speeds_tag::type& characteristic_speeds) noexcept {
105  expand_pack((*packaged_u = u_to_package)...);
106  expand_pack((*packaged_n_dot_f = n_dot_f_to_package)...);
107 
108  *packaged_largest_ingoing_speed = make_with_value<Scalar<DataVector>>(
109  characteristic_speeds[0],
110  std::numeric_limits<double>::signaling_NaN());
111  *packaged_largest_outgoing_speed = make_with_value<Scalar<DataVector>>(
112  characteristic_speeds[0],
113  std::numeric_limits<double>::signaling_NaN());
114 
115  // When packaging interior data, LargestIngoingSpeed and
116  // LargestOutgoingSpeed will hold the min and max char speeds,
117  // respectively. On the other hand, when packaging exterior data, the
118  // characteristic speeds will be computed along *minus* the exterior
119  // normal, so LargestIngoingSpeed will hold *minus* the max speed, while
120  // LargestOutgoingSpeed will store *minus* the min speed.
121  for (size_t s = 0; s < characteristic_speeds[0].size(); ++s) {
122  get(*packaged_largest_ingoing_speed)[s] = (*std::min_element(
123  characteristic_speeds.begin(), characteristic_speeds.end(),
124  [&s](const DataVector& a, const DataVector& b) noexcept {
125  return a[s] < b[s];
126  }))[s];
127  get(*packaged_largest_outgoing_speed)[s] = (*std::max_element(
128  characteristic_speeds.begin(), characteristic_speeds.end(),
129  [&s](const DataVector& a, const DataVector& b) noexcept {
130  return a[s] < b[s];
131  }))[s];
132  }
133  }
134  };
135 
136  template <typename NormalDotNumericalFluxTagList, typename VariablesTagList,
137  typename NormalDotFluxTagList>
138  struct call_operator_helper;
139  template <typename... NormalDotNumericalFluxTags, typename... VariablesTags,
140  typename... NormalDotFluxTags>
141  struct call_operator_helper<tmpl::list<NormalDotNumericalFluxTags...>,
142  tmpl::list<VariablesTags...>,
143  tmpl::list<NormalDotFluxTags...>> {
144  static void function(
145  const gsl::not_null<
146  typename NormalDotNumericalFluxTags::type*>... n_dot_numerical_f,
147  const typename NormalDotFluxTags::type&... n_dot_f_interior,
148  const typename VariablesTags::type&... u_interior,
149  const Scalar<DataVector>& largest_ingoing_speed_interior,
150  const typename LargestOutgoingSpeed::type&
151  largest_outgoing_speed_interior,
152  const typename NormalDotFluxTags::type&... minus_n_dot_f_exterior,
153  const typename VariablesTags::type&... u_exterior,
154  // names are inverted w.r.t. interior data. See package_data()
155  const typename LargestIngoingSpeed::type&
156  minus_largest_outgoing_speed_exterior,
157  const typename LargestOutgoingSpeed::type&
158  minus_largest_ingoing_speed_exterior) noexcept {
159  const auto number_of_grid_points =
160  get(largest_ingoing_speed_interior).size();
161  auto largest_ingoing_speed = Scalar<DataVector>{number_of_grid_points};
162  auto largest_outgoing_speed = Scalar<DataVector>{number_of_grid_points};
163  for (size_t s = 0; s < largest_ingoing_speed.begin()->size(); ++s) {
164  get(largest_ingoing_speed)[s] =
165  std::min({get(largest_ingoing_speed_interior)[s],
166  -get(minus_largest_ingoing_speed_exterior)[s], 0.0});
167  get(largest_outgoing_speed)[s] =
168  std::max({get(largest_outgoing_speed_interior)[s],
169  -get(minus_largest_outgoing_speed_exterior)[s], 0.0});
170  }
171  ASSERT(
172  min(get(largest_outgoing_speed) - get(largest_ingoing_speed)) > 0.0,
173  "Max, min speeds are the same:\n"
174  " largest_outgoing_speed = "
175  << get(largest_outgoing_speed)
176  << "\n"
177  " largest_ingoing_speed = "
178  << get(largest_ingoing_speed));
179  const DataVector one_over_sp_minus_sm =
180  1.0 / (get(largest_outgoing_speed) - get(largest_ingoing_speed));
181  const auto assemble_numerical_flux =
182  [
183  &largest_ingoing_speed, &largest_outgoing_speed,
184  &one_over_sp_minus_sm
185  ](const auto n_dot_num_f, const auto& n_dot_f_in, const auto& u_in,
186  const auto& minus_n_dot_f_ex, const auto& u_ex) noexcept {
187  for (size_t i = 0; i < n_dot_num_f->size(); ++i) {
188  (*n_dot_num_f)[i] =
189  one_over_sp_minus_sm *
190  (get(largest_outgoing_speed) * n_dot_f_in[i] +
191  get(largest_ingoing_speed) * minus_n_dot_f_ex[i] -
192  get(largest_outgoing_speed) * get(largest_ingoing_speed) *
193  (u_in[i] - u_ex[i]));
194  }
195  return nullptr;
196  };
197  expand_pack(assemble_numerical_flux(n_dot_numerical_f, n_dot_f_interior,
198  u_interior, minus_n_dot_f_exterior,
199  u_exterior)...);
200  }
201  };
202 
203  public:
204  using options = tmpl::list<>;
205  static constexpr Options::String help = {"Computes the HLL numerical flux."};
206 
207  // clang-tidy: google-runtime-references
208  void pup(PUP::er& /*p*/) noexcept {} // NOLINT
209 
210  template <class... Args>
211  void package_data(const Args&... args) const noexcept {
212  package_data_helper<variables_tags,
213  db::wrap_tags_in<::Tags::NormalDotFlux,
214  variables_tags>>::function(args...);
215  }
216 
217  template <class... Args>
218  void operator()(const Args&... args) const noexcept {
219  call_operator_helper<
220  db::wrap_tags_in<::Tags::NormalDotNumericalFlux, variables_tags>,
221  variables_tags,
222  db::wrap_tags_in<::Tags::NormalDotFlux,
223  variables_tags>>::function(args...);
224  }
225 };
226 
227 } // namespace NumericalFluxes
228 } // namespace dg
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:585
Options.hpp
std::size
T size(T... args)
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::NormalDotFlux
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:96
algorithm
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyMassMatrix.hpp:13
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
cstddef
MakeWithValue.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
std::min
T min(T... args)
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
dg::NumericalFluxes::Hll
Compute the HLL numerical flux.
Definition: Hll.hpp:59
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Tensor.hpp
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
dg::NumericalFluxes::Hll::LargestOutgoingSpeed
Estimate for the largest outgoing signal speed.
Definition: Hll.hpp:71
dg::NumericalFluxes::Hll::LargestIngoingSpeed
Estimate for the largest ingoing signal speed.
Definition: Hll.hpp:66
Prefixes.hpp
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
© Copyright 2017 - 2021 SXS Collaboration, MIT License