Minmod.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstdlib>
8 #include <limits>
9 #include <type_traits>
10 #include <unordered_map>
11 #include <utility>
12 
13 #include "DataStructures/Tags.hpp"
15 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodImpl.hpp"
16 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
17 #include "Options/Options.hpp"
18 #include "Utilities/Gsl.hpp"
19 #include "Utilities/MakeArray.hpp"
20 #include "Utilities/TMPL.hpp"
21 #include "Utilities/TaggedTuple.hpp"
22 
23 /// \cond
24 class DataVector;
25 template <size_t VolumeDim>
26 class Direction;
27 template <size_t Dim, typename T>
28 class DirectionMap;
29 template <size_t VolumeDim>
30 class Element;
31 template <size_t VolumeDim>
32 class ElementId;
33 template <size_t VolumeDim>
34 class Mesh;
35 template <size_t VolumeDim>
36 class OrientationMap;
37 
38 namespace boost {
39 template <class T>
40 struct hash;
41 } // namespace boost
42 
43 namespace PUP {
44 class er;
45 } // namespace PUP
46 
47 namespace Limiters {
48 template <size_t VolumeDim, typename TagsToLimit>
49 class Minmod;
50 
51 namespace Minmod_detail {
52 template <size_t VolumeDim>
53 class BufferWrapper;
54 } // namespace Minmod_detail
55 } // namespace Limiters
56 
57 namespace domain::Tags {
58 template <size_t Dim, typename Frame>
59 struct Coordinates;
60 template <size_t VolumeDim>
61 struct Element;
62 template <size_t VolumeDim>
63 struct Mesh;
64 template <size_t VolumeDim>
65 struct SizeOfElement;
66 } // namespace domain::Tags
67 /// \endcond
68 
69 namespace Limiters {
70 /// \ingroup LimitersGroup
71 /// \brief A general minmod slope limiter
72 ///
73 /// Provides an implementation for the three minmod-based generalized slope
74 /// limiters from \cite Cockburn1999 Sec. 2.4: \f$\Lambda\Pi^1\f$,
75 /// \f$\Lambda\Pi^N\f$, and MUSCL. Below we summarize these three limiters, but
76 /// the reader should refer to the reference for full details. The limiter has a
77 /// general implementation that can work on an arbitrary set of tensors; the
78 /// limiting algorithm is applied to each component of each tensor
79 /// independently.
80 ///
81 /// The MUSCL and \f$\Lambda\Pi^1\f$ limiters are both intended for use on
82 /// piecewise-linear solutions, i.e., on linear-order elements with two points
83 /// per dimension. These limiters operate by reducing the spatial slope of the
84 /// tensor components if the data look like they may contain oscillations.
85 /// Between these two, MUSCL is more dissipative --- it more aggressively
86 /// reduces the slopes of the data, so it may better handle strong shocks, but
87 /// it correspondingly produces the most broadening of features in the solution.
88 ///
89 /// Note that we do not _require_ the MUSCL and \f$\Lambda\Pi^1\f$ limiters to
90 /// be used on linear-order elements. However, when they are used on a
91 /// higher-resolution grid, the limiters act to linearize the solution (by
92 /// discarding higher-order mode content) whether or not the slopes must be
93 /// reduced.
94 ///
95 /// The \f$\Lambda\Pi^N\f$ limiter is intended for use with higher-order
96 /// elements (with more than two points per dimension), where the solution is a
97 /// piecewise polynomial of higher-than-linear order. This limiter generalizes
98 /// \f$\Lambda\Pi^1\f$: the post-limiter solution is the linearized solution of
99 /// \f$\Lambda\Pi^1\f$ in the case that the slopes must be reduced, but is the
100 /// original (higher-order) data in the case that the slopes are acceptable.
101 ///
102 /// For all three types of minmod limiter, the algorithm can be relaxed from
103 /// TVD (total variation diminishing) in the means to TVB (total variation
104 /// bound) in the means. This may avoid limiting away smooth extrema in the
105 /// solution that would otherwise look like spurious oscillations. When this
106 /// correction is enabled, the limiter will not reduce the slope (but may still
107 /// linearize) on elements where the slope is less than \f$m h^2\f$, where
108 /// \f$m\f$ is the TVB constant and \f$h\f$ is the size of the DG element.
109 /// Note the "in the means" qualifier: the limiter controls the oscillation
110 /// between the mean solution values across neighboring cells, but may not
111 /// control oscillations within the cells.
112 ///
113 /// The choice of the TVB constant \f$m\f$ is difficult. Larger values result in
114 /// fewer limiter activations, especially near smooth extrema in the solution
115 /// --- this can help to avoid incorrectly limiting away these smooth extrema,
116 /// but can also result in insufficient limiting of truly spurious oscillations.
117 /// The reference uses a value of 50 when presenting the limiter with simple
118 /// shock tests, but in general the value \f$m\f$ that optimizes between
119 /// robustness and minimal loss of accuracy is problem dependent.
120 ///
121 /// The limiter acts in the `Frame::Logical` coordinates, because in these
122 /// coordinates it is straightforward to formulate the algorithm. This means the
123 /// limiter can operate on generic deformed grids. However, if the grid is too
124 /// strongly deformed, some things can start to break down:
125 /// 1. When an element is deformed so that the Jacobian (from `Frame::Logical`
126 /// to `Frame::Inertial`) varies across the element, then the limiter fails
127 /// to be conservative. In other words, the integral of a tensor `u` over the
128 /// element will change after the limiter activates on `u`. This error is
129 /// typically small.
130 /// 2. When there is a sudden change in the size of the elements (perhaps at an
131 /// h-refinement boundary, or at the boundary between two blocks with very
132 /// different mappings), a smooth solution in `Frame::Inertial` can appear
133 /// to have a kink in `Frame::Logical`. The Minmod implementation includes
134 /// some (untested) tweaks that try to reduce spurious limiter activations
135 /// near these fake kinks.
136 ///
137 /// When an element has multiple neighbors in any direction, an effective mean
138 /// and neighbor size in this direction are computed by averaging over the
139 /// multiple neighbors. This simple generalization of the minmod limiter enables
140 /// it to operate on h-refined grids.
141 ///
142 /// \tparam VolumeDim The number of spatial dimensions.
143 /// \tparam Tags A typelist of tags specifying the tensors to limit.
144 template <size_t VolumeDim, typename... Tags>
145 class Minmod<VolumeDim, tmpl::list<Tags...>> {
146  public:
147  /// \brief The MinmodType
148  ///
149  /// One of `Limiters::MinmodType`. See `Limiters::Minmod`
150  /// documentation for details. Note in particular that on grids with more than
151  /// two points per dimension, the recommended type is `LambdaPiN`.
152  struct Type {
153  using type = MinmodType;
154  static constexpr Options::String help = {"Type of minmod"};
155  };
156  /// \brief The TVB constant
157  ///
158  /// See `Limiters::Minmod` documentation for details. The optimal value of
159  /// this parameter is unfortunately problem-dependent.
160  struct TvbConstant {
161  using type = double;
162  static type lower_bound() noexcept { return 0.0; }
163  static constexpr Options::String help = {"TVB constant 'm'"};
164  };
165  /// \brief Turn the limiter off
166  ///
167  /// This option exists to temporarily disable the limiter for debugging
168  /// purposes. For problems where limiting is not needed, the preferred
169  /// approach is to not compile the limiter into the executable.
170  struct DisableForDebugging {
171  using type = bool;
172  static type suggested_value() noexcept { return false; }
173  static constexpr Options::String help = {"Disable the limiter"};
174  };
175  using options = tmpl::list<Type, TvbConstant, DisableForDebugging>;
176  static constexpr Options::String help = {
177  "A minmod-based slope limiter.\n"
178  "The different types of minmod are more or less aggressive in trying\n"
179  "to reduce slopes. The TVB correction allows the limiter to ignore\n"
180  "'small' slopes, and helps to avoid limiting of smooth extrema in the\n"
181  "solution.\n"};
182 
183  /// \brief Constuct a Minmod slope limiter
184  ///
185  /// \param minmod_type The type of Minmod slope limiter.
186  /// \param tvb_constant The value of the TVB constant.
187  /// \param disable_for_debugging Switch to turn the limiter off.
188  explicit Minmod(MinmodType minmod_type, double tvb_constant,
189  bool disable_for_debugging = false) noexcept;
190 
191  Minmod() noexcept = default;
192  Minmod(const Minmod& /*rhs*/) = default;
193  Minmod& operator=(const Minmod& /*rhs*/) = default;
194  Minmod(Minmod&& /*rhs*/) noexcept = default;
195  Minmod& operator=(Minmod&& /*rhs*/) noexcept = default;
196  ~Minmod() = default;
197 
198  // clang-tidy: google-runtime-references
199  void pup(PUP::er& p) noexcept; // NOLINT
200 
201  // To facilitate testing
202  /// \cond
203  MinmodType minmod_type() const noexcept { return minmod_type_; }
204  /// \endcond
205 
206  /// \brief Data to send to neighbor elements.
207  struct PackagedData {
209  std::array<double, VolumeDim> element_size =
210  make_array<VolumeDim>(std::numeric_limits<double>::signaling_NaN());
211 
212  // clang-tidy: google-runtime-references
213  void pup(PUP::er& p) noexcept { // NOLINT
214  p | means;
215  p | element_size;
216  }
217  };
218 
219  using package_argument_tags =
220  tmpl::list<Tags..., domain::Tags::Mesh<VolumeDim>,
222 
223  /// \brief Package data for sending to neighbor elements.
224  ///
225  /// The following quantities are stored in `PackagedData` and communicated
226  /// between neighboring elements:
227  /// - the cell-averaged mean of each tensor component, and
228  /// - the size of the cell along each logical coordinate direction.
229  ///
230  /// \param packaged_data The data package to fill with this element's values.
231  /// \param tensors The tensors to be averaged and packaged.
232  /// \param mesh The mesh on which the tensor values are measured.
233  /// \param element_size The size of the element in inertial coordinates, along
234  /// each dimension of logical coordinates.
235  /// \param orientation_map The orientation of the neighbor
236  void package_data(
237  gsl::not_null<PackagedData*> packaged_data,
238  const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
239  const std::array<double, VolumeDim>& element_size,
240  const OrientationMap<VolumeDim>& orientation_map) const noexcept;
241 
242  using limit_tags = tmpl::list<Tags...>;
243  using limit_argument_tags =
244  tmpl::list<domain::Tags::Mesh<VolumeDim>,
248 
249  /// \brief Limits the solution on the element.
250  ///
251  /// For each component of each tensor, the limiter will (in general) linearize
252  /// the data, then possibly reduce its slope, dimension-by-dimension, until it
253  /// no longer looks oscillatory.
254  ///
255  /// \param tensors The tensors to be limited.
256  /// \param element The element on which the tensors to limit live.
257  /// \param mesh The mesh on which the tensor values are measured.
258  /// \param logical_coords The logical coordinates of the mesh gridpoints.
259  /// \param element_size The size of the element, in the inertial coordinates.
260  /// \param neighbor_data The data from each neighbor.
261  ///
262  /// \return whether the limiter modified the solution or not.
263  ///
264  /// \note The return value is false if the limiter knows it has not modified
265  /// the solution. True return values can indicate:
266  /// - The solution was limited to reduce the slope, whether by a large factor
267  /// or by a factor only roundoff away from unity.
268  /// - The solution was linearized but not limited.
269  /// - The solution is identical to the input, if the input was a linear
270  /// function on a higher-order mesh, so that the limiter cannot know that
271  /// the linearization step did not actually modify the data. This is
272  /// somewhat contrived and is unlikely to occur outside of code tests or
273  /// test cases with very clean initial data.
274  bool operator()(
276  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
277  const tnsr::I<DataVector, VolumeDim, Frame::Logical>& logical_coords,
278  const std::array<double, VolumeDim>& element_size,
279  const std::unordered_map<
282  neighbor_data) const noexcept;
283 
284  private:
285  template <size_t LocalDim, typename LocalTagList>
286  // NOLINTNEXTLINE(readability-redundant-declaration) false positive
287  friend bool operator==(const Minmod<LocalDim, LocalTagList>& lhs,
288  const Minmod<LocalDim, LocalTagList>& rhs) noexcept;
289 
290  MinmodType minmod_type_;
291  double tvb_constant_;
292  bool disable_for_debugging_;
293 };
294 
295 template <size_t VolumeDim, typename... Tags>
296 Minmod<VolumeDim, tmpl::list<Tags...>>::Minmod(
297  const MinmodType minmod_type, const double tvb_constant,
298  const bool disable_for_debugging) noexcept
299  : minmod_type_(minmod_type),
300  tvb_constant_(tvb_constant),
301  disable_for_debugging_(disable_for_debugging) {
302  ASSERT(tvb_constant >= 0.0, "The TVB constant must be non-negative.");
303 }
304 
305 template <size_t VolumeDim, typename... Tags>
306 void Minmod<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) noexcept {
307  p | minmod_type_;
308  p | tvb_constant_;
309  p | disable_for_debugging_;
310 }
311 
312 template <size_t VolumeDim, typename... Tags>
313 void Minmod<VolumeDim, tmpl::list<Tags...>>::package_data(
314  const gsl::not_null<PackagedData*> packaged_data,
315  const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
316  const std::array<double, VolumeDim>& element_size,
317  const OrientationMap<VolumeDim>& orientation_map) const noexcept {
318  if (UNLIKELY(disable_for_debugging_)) {
319  // Do not initialize packaged_data
320  return;
321  }
322 
323  const auto wrap_compute_means = [&mesh, &packaged_data](
324  auto tag, const auto tensor) noexcept {
325  for (size_t i = 0; i < tensor.size(); ++i) {
326  // Compute the mean using the local orientation of the tensor and mesh:
327  // this avoids the work of reorienting the tensor while giving the same
328  // result.
329  get<::Tags::Mean<decltype(tag)>>(packaged_data->means)[i] =
330  mean_value(tensor[i], mesh);
331  }
332  return '0';
333  };
334  expand_pack(wrap_compute_means(Tags{}, tensors)...);
335  packaged_data->element_size =
336  orientation_map.permute_from_neighbor(element_size);
337 }
338 
339 template <size_t VolumeDim, typename... Tags>
340 bool Minmod<VolumeDim, tmpl::list<Tags...>>::operator()(
342  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
343  const tnsr::I<DataVector, VolumeDim, Frame::Logical>& logical_coords,
344  const std::array<double, VolumeDim>& element_size,
345  const std::unordered_map<
347  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
348  neighbor_data) const noexcept {
349  if (UNLIKELY(disable_for_debugging_)) {
350  // Do not modify input tensors
351  return false;
352  }
353 
354  DataVector u_lin_buffer(mesh.number_of_grid_points());
355  Minmod_detail::BufferWrapper<VolumeDim> buffer(mesh);
356 
357  bool limiter_activated = false;
358  const auto wrap_minmod_impl = [this, &limiter_activated, &element, &mesh,
359  &logical_coords, &element_size, &neighbor_data,
360  &u_lin_buffer, &buffer](
361  auto tag, const auto tensor) noexcept {
362  limiter_activated =
363  Minmod_detail::minmod_impl<VolumeDim, decltype(tag)>(
364  &u_lin_buffer, &buffer, tensor, minmod_type_, tvb_constant_, mesh,
365  element, logical_coords, element_size, neighbor_data) or
366  limiter_activated;
367  return '0';
368  };
369  expand_pack(wrap_minmod_impl(Tags{}, tensors)...);
370  return limiter_activated;
371 }
372 
373 template <size_t LocalDim, typename LocalTagList>
374 bool operator==(const Minmod<LocalDim, LocalTagList>& lhs,
375  const Minmod<LocalDim, LocalTagList>& rhs) noexcept {
376  return lhs.minmod_type_ == rhs.minmod_type_ and
377  lhs.tvb_constant_ == rhs.tvb_constant_ and
378  lhs.disable_for_debugging_ == rhs.disable_for_debugging_;
379 }
380 
381 template <size_t VolumeDim, typename TagList>
382 bool operator!=(const Minmod<VolumeDim, TagList>& lhs,
383  const Minmod<VolumeDim, TagList>& rhs) noexcept {
384  return not(lhs == rhs);
385 }
386 
387 } // namespace Limiters
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:549
Limiters
Things relating to limiting.
Definition: HwenoImpl.hpp:42
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
domain::Tags::Coordinates< VolumeDim, Frame::Logical >
std::pair
Options.hpp
MakeArray.hpp
domain::Tags::Element
Definition: Tags.hpp:97
Tags::Mean
Given the Tag holding a Tensor<DataVector, ...>, swap the DataVector with a double.
Definition: Tags.hpp:23
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
OrientationMap
A mapping of the logical coordinate axes of a host to the logical coordinate axes of a neighbor of th...
Definition: OrientationMap.hpp:34
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
domain::Tags::SizeOfElement
Definition: SizeOfElement.hpp:77
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:50
std::add_pointer_t
array
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Limiters::MinmodType
MinmodType
Possible types of the minmod slope limiter and/or troubled-cell indicator.
Definition: MinmodType.hpp:20
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
mean_value
double mean_value(const DataVector &f, const Mesh< Dim > &mesh) noexcept
Compute the mean value of a function over a manifold.
Definition: MeanValue.hpp:47
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:48
DirectionMap
Definition: DirectionMap.hpp:15
cstdlib
limits
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Tensor.hpp
unordered_map
std::numeric_limits
type_traits
TMPL.hpp
domain::Tags
Tags for the domain.
Definition: FaceNormal.hpp:107
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13