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