Weno.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <boost/functional/hash.hpp> // IWYU pragma: keep
8 #include <cstddef>
9 #include <limits>
10 #include <string>
11 #include <type_traits>
12 #include <unordered_map>
13 #include <utility>
14 
17 #include "DataStructures/Tags.hpp"
20 #include "Domain/Element.hpp" // IWYU pragma: keep
21 #include "Domain/Mesh.hpp"
22 #include "Domain/OrientationMapHelpers.hpp"
23 #include "Domain/SizeOfElement.hpp"
24 #include "Domain/Tags.hpp" // IWYU pragma: keep
25 #include "ErrorHandling/Assert.hpp"
26 #include "Evolution/DiscontinuousGalerkin/Limiters/HwenoModifiedSolution.hpp"
27 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodTci.hpp"
28 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
29 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
30 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
31 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoType.hpp"
32 #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
34 #include "Options/Options.hpp"
35 #include "Utilities/Algorithm.hpp"
36 #include "Utilities/Gsl.hpp"
37 #include "Utilities/MakeArray.hpp"
38 #include "Utilities/TMPL.hpp"
40 
41 /// \cond
42 template <size_t VolumeDim>
43 class Direction;
44 template <size_t VolumeDim>
45 class ElementId;
46 
47 namespace PUP {
48 class er;
49 } // namespace PUP
50 
51 namespace Limiters {
52 template <size_t VolumeDim, typename TagsToLimit>
53 class Weno;
54 } // namespace Limiters
55 /// \endcond
56 
57 namespace Limiters {
58 /// \ingroup LimitersGroup
59 /// \brief A compact-stencil WENO limiter for DG
60 ///
61 /// Implements the simple WENO limiter of \cite Zhong2013 and the Hermite WENO
62 /// (HWENO) limiter of \cite Zhu2016. These limiters require communication only
63 /// between nearest-neighbor elements, but preserve the full order of the DG
64 /// solution when the solution is smooth. Full volume data is communicated
65 /// between neighbors.
66 ///
67 /// The limiter uses a Minmod-based troubled-cell indicator to identify elements
68 /// that need limiting. The \f$\Lambda\Pi^N\f$ limiter of \cite Cockburn1999 is
69 /// used. Note that the HWENO paper recommends a more sophisticated
70 /// troubled-cell indicator instead, but the specific choice of indicator should
71 /// not be too important for a high-order WENO limiter.
72 ///
73 /// On any identified "troubled" elements, the limited solution is obtained by
74 /// WENO reconstruction --- a linear combination of the local DG solution and a
75 /// "modified" solution from each neighbor element. For the simple WENO limiter,
76 /// the modified solution is obtained by simply extrapolating the neighbor
77 /// solution onto the troubled element. For the HWENO limiter, the modified
78 /// solution is obtained by a least-squares fit to the solution across multiple
79 /// neighboring elements.
80 ///
81 /// To reconstruct the WENO solution from the local solution and the modified
82 /// neighbor solutions, the standard WENO procedure is followed. We use the
83 /// oscillation indicator of \cite Dumbser2007, modified for use on the
84 /// square/cube grids of SpECTRE. We favor this indicator because portions of
85 /// the work can be precomputed, leading to an oscillation measure that is
86 /// efficient to evaluate.
87 ///
88 /// \warning
89 /// Limitations:
90 /// - Does not support h- or p-refinement; this is enforced with ASSERTIONS.
91 /// - Does not support curved elements; this is not enforced.
92 template <size_t VolumeDim, typename... Tags>
93 class Weno<VolumeDim, tmpl::list<Tags...>> {
94  public:
95  /// \brief The WenoType
96  ///
97  /// One of `Limiters::WenoType`. See the `Limiters::Weno`
98  /// documentation for details.
99  struct Type {
100  using type = WenoType;
101  static constexpr OptionString help = {"Type of WENO limiter"};
102  };
103  /// \brief The linear weight given to each neighbor
104  ///
105  /// This linear weight gets combined with the oscillation indicator to
106  /// compute the weight for each WENO estimated solution. Larger values are
107  /// better suited for problems with strong shocks, and smaller values are
108  /// better suited to smooth problems.
109  struct NeighborWeight {
110  using type = double;
111  static type default_value() noexcept { return 0.001; }
112  static type lower_bound() noexcept { return 1e-6; }
113  static type upper_bound() noexcept { return 0.1; }
114  static constexpr OptionString help = {
115  "Linear weight for each neighbor element's solution"};
116  };
117  /// \brief Turn the limiter off
118  ///
119  /// This option exists to temporarily disable the limiter for debugging
120  /// purposes. For problems where limiting is not needed, the preferred
121  /// approach is to not compile the limiter into the executable.
122  struct DisableForDebugging {
123  using type = bool;
124  static type default_value() noexcept { return false; }
125  static constexpr OptionString help = {"Disable the limiter"};
126  };
127  using options = tmpl::list<Type, NeighborWeight, DisableForDebugging>;
128  static constexpr OptionString help = {"A WENO limiter for DG"};
129 
130  Weno(WenoType weno_type, double neighbor_linear_weight,
131  bool disable_for_debugging = false) noexcept;
132 
133  Weno() noexcept = default;
134  Weno(const Weno& /*rhs*/) = default;
135  Weno& operator=(const Weno& /*rhs*/) = default;
136  Weno(Weno&& /*rhs*/) noexcept = default;
137  Weno& operator=(Weno&& /*rhs*/) noexcept = default;
138  ~Weno() = default;
139 
140  // NOLINTNEXTLINE(google-runtime-references)
141  void pup(PUP::er& p) noexcept;
142 
143  /// \brief Data to send to neighbor elements
144  struct PackagedData {
145  Variables<tmpl::list<Tags...>> volume_data;
147  Mesh<VolumeDim> mesh;
148  std::array<double, VolumeDim> element_size =
149  make_array<VolumeDim>(std::numeric_limits<double>::signaling_NaN());
150 
151  // NOLINTNEXTLINE(google-runtime-references)
152  void pup(PUP::er& p) noexcept {
153  p | volume_data;
154  p | means;
155  p | mesh;
156  p | element_size;
157  }
158  };
159 
160  using package_argument_tags = tmpl::list<Tags..., ::Tags::Mesh<VolumeDim>,
162 
163  /// \brief Package data for sending to neighbor elements
164  void package_data(gsl::not_null<PackagedData*> packaged_data,
165  const db::item_type<Tags>&... tensors,
166  const Mesh<VolumeDim>& mesh,
167  const std::array<double, VolumeDim>& element_size,
168  const OrientationMap<VolumeDim>& orientation_map) const
169  noexcept;
170 
171  using limit_tags = tmpl::list<Tags...>;
172  using limit_argument_tags =
173  tmpl::list<::Tags::Element<VolumeDim>, ::Tags::Mesh<VolumeDim>,
174  ::Tags::SizeOfElement<VolumeDim>>;
175 
176  /// \brief Limit the solution on the element
177  bool operator()(
179  const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
180  const std::array<double, VolumeDim>& element_size,
181  const std::unordered_map<
184  neighbor_data) const noexcept;
185 
186  private:
187  template <size_t LocalDim, typename LocalTagList>
188  // NOLINTNEXTLINE(readability-redundant-declaration) false positive
189  friend bool operator==(const Weno<LocalDim, LocalTagList>& lhs,
190  const Weno<LocalDim, LocalTagList>& rhs) noexcept;
191 
192  WenoType weno_type_;
193  double neighbor_linear_weight_;
194  bool disable_for_debugging_;
195 };
196 
197 template <size_t VolumeDim, typename... Tags>
198 Weno<VolumeDim, tmpl::list<Tags...>>::Weno(
199  const WenoType weno_type, const double neighbor_linear_weight,
200  const bool disable_for_debugging) noexcept
201  : weno_type_(weno_type),
202  neighbor_linear_weight_(neighbor_linear_weight),
203  disable_for_debugging_(disable_for_debugging) {}
204 
205 template <size_t VolumeDim, typename... Tags>
206 // NOLINTNEXTLINE(google-runtime-references)
207 void Weno<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) noexcept {
208  p | weno_type_;
209  p | neighbor_linear_weight_;
210  p | disable_for_debugging_;
211 }
212 
213 template <size_t VolumeDim, typename... Tags>
214 void Weno<VolumeDim, tmpl::list<Tags...>>::package_data(
215  const gsl::not_null<PackagedData*> packaged_data,
216  const db::item_type<Tags>&... tensors, const Mesh<VolumeDim>& mesh,
217  const std::array<double, VolumeDim>& element_size,
218  const OrientationMap<VolumeDim>& orientation_map) const noexcept {
219  if (UNLIKELY(disable_for_debugging_)) {
220  // Do not initialize packaged_data
221  return;
222  }
223 
224  const auto wrap_compute_means =
225  [&mesh, &packaged_data ](auto tag, const auto& tensor) noexcept {
226  for (size_t i = 0; i < tensor.size(); ++i) {
227  // Compute the mean using the local orientation of the tensor and mesh.
228  get<::Tags::Mean<decltype(tag)>>(packaged_data->means)[i] =
229  mean_value(tensor[i], mesh);
230  }
231  return '0';
232  };
233  expand_pack(wrap_compute_means(Tags{}, tensors)...);
234 
235  packaged_data->element_size =
236  orientation_map.permute_from_neighbor(element_size);
237 
238  (packaged_data->volume_data).initialize(mesh.number_of_grid_points());
239  const auto wrap_copy_tensor = [&packaged_data](auto tag,
240  const auto& tensor) noexcept {
241  get<decltype(tag)>(packaged_data->volume_data) = tensor;
242  return '0';
243  };
244  expand_pack(wrap_copy_tensor(Tags{}, tensors)...);
245  packaged_data->volume_data = orient_variables(
246  packaged_data->volume_data, mesh.extents(), orientation_map);
247 
248  packaged_data->mesh = orientation_map(mesh);
249 }
250 
251 template <size_t VolumeDim, typename... Tags>
252 bool Weno<VolumeDim, tmpl::list<Tags...>>::operator()(
254  const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
255  const std::array<double, VolumeDim>& element_size,
256  const std::unordered_map<
258  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
259  neighbor_data) const noexcept {
260  if (UNLIKELY(disable_for_debugging_)) {
261  // Do not modify input tensors
262  return false;
263  }
264 
265  // Enforce restrictions on h-refinement, p-refinement
266  if (UNLIKELY(alg::any_of(element.neighbors(),
267  [](const auto& direction_neighbors) noexcept {
268  return direction_neighbors.second.size() != 1;
269  }))) {
270  ERROR("The Weno limiter does not yet support h-refinement");
271  // Removing this limitation will require:
272  // - Generalizing the computation of the modified neighbor solutions.
273  // - Generalizing the WENO weighted sum for multiple neighbors in each
274  // direction.
275  }
276  alg::for_each(neighbor_data, [&mesh](const auto& neighbor_and_data) noexcept {
277  if (UNLIKELY(neighbor_and_data.second.mesh != mesh)) {
278  ERROR("The Weno limiter does not yet support p-refinement");
279  // Removing this limitation will require generalizing the
280  // computation of the modified neighbor solutions.
281  }
282  });
283 
284  // Troubled-cell detection for WENO flags the element for limiting if any
285  // component of any tensor needs limiting.
286  const auto minmod_tci_type = MinmodType::LambdaPiN;
287  const double minmod_tci_tvbm_constant = 0.0;
288  const bool cell_is_troubled =
289  Minmod_detail::troubled_cell_indicator<VolumeDim, PackagedData, Tags...>(
290  (*tensors)..., neighbor_data, minmod_tci_type,
291  minmod_tci_tvbm_constant, element, mesh, element_size);
292 
293  if (not cell_is_troubled) {
294  // No limiting is needed
295  return false;
296  }
297 
298  // Compute the modified solutions from each neighbor, for each tensor
299  // component. For this step, each WenoType requires a different treatment.
301  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>,
302  Variables<tmpl::list<Tags...>>,
303  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>
304  modified_neighbor_solutions;
305 
306  if (weno_type_ == WenoType::Hweno) {
307  // For each neighbor, the HWENO fits are done one tensor at a time.
308  for (const auto& neighbor_and_data : neighbor_data) {
309  modified_neighbor_solutions[neighbor_and_data.first].initialize(
310  mesh.number_of_grid_points());
311  }
312  const auto wrap_hweno_neighbor_solution_one_tensor =
313  [&element, &mesh, &neighbor_data, &
314  modified_neighbor_solutions ](auto tag, const auto& tensor) noexcept {
315  for (const auto& neighbor_and_data : neighbor_data) {
316  const auto& primary_neighbor = neighbor_and_data.first;
317  auto& modified_tensor = get<decltype(tag)>(
318  modified_neighbor_solutions.at(primary_neighbor));
319  hweno_modified_neighbor_solution<decltype(tag)>(
320  make_not_null(&modified_tensor), *tensor, element, mesh,
321  neighbor_data, primary_neighbor);
322  }
323  return '0';
324  };
325  expand_pack(wrap_hweno_neighbor_solution_one_tensor(Tags{}, tensors)...);
326  } else if (weno_type_ == WenoType::SimpleWeno) {
327  // For each neighbor, the simple WENO data is simply extrapolated. This
328  // can be done on the entire Variables of tags to limit.
329  for (const auto& neighbor_and_data : neighbor_data) {
330  const auto& neighbor = neighbor_and_data.first;
331  const auto& direction = neighbor.first;
332  const auto& data = neighbor_and_data.second;
333  const auto& source_mesh = data.mesh;
334  const auto target_1d_logical_coords =
335  Weno_detail::local_grid_points_in_neighbor_logical_coords(
336  mesh, source_mesh, element, direction);
337  const intrp::RegularGrid<VolumeDim> interpolant(source_mesh, mesh,
338  target_1d_logical_coords);
339  modified_neighbor_solutions.insert(
340  std::make_pair(neighbor, interpolant.interpolate(data.volume_data)));
341  }
342  } else {
343  ERROR("WENO limiter not implemented for WenoType: " << weno_type_);
344  }
345 
346  // Reconstruct WENO solution from local solution and modified neighbor
347  // solutions.
348  const auto wrap_reconstruct_one_tensor =
349  [ this, &mesh, &
350  modified_neighbor_solutions ](auto tag, const auto& tensor) noexcept {
351  Weno_detail::reconstruct_from_weighted_sum<decltype(tag)>(
352  tensor, mesh, neighbor_linear_weight_, modified_neighbor_solutions);
353  return '0';
354  };
355  expand_pack(wrap_reconstruct_one_tensor(Tags{}, tensors)...);
356  return true; // cell_is_troubled
357 }
358 
359 template <size_t LocalDim, typename LocalTagList>
360 bool operator==(const Weno<LocalDim, LocalTagList>& lhs,
361  const Weno<LocalDim, LocalTagList>& rhs) noexcept {
362  return lhs.weno_type_ == rhs.weno_type_ and
363  lhs.neighbor_linear_weight_ == rhs.neighbor_linear_weight_ and
364  lhs.disable_for_debugging_ == rhs.disable_for_debugging_;
365 }
366 
367 template <size_t VolumeDim, typename TagList>
368 bool operator!=(const Weno<VolumeDim, TagList>& lhs,
369  const Weno<VolumeDim, TagList>& rhs) noexcept {
370  return not(lhs == rhs);
371 }
372 
373 } // namespace Limiters
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
decltype(auto) any_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::any_of.
Definition: Algorithm.hpp:148
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
Variables< TagsList > orient_variables(const Variables< TagsList > &variables, const Index< VolumeDim > &extents, const OrientationMap< VolumeDim > &orientation_of_neighbor) noexcept
Orient variables to the data-storage order of a neighbor element with the given orientation.
Definition: OrientationMapHelpers.hpp:84
Definition: Digraph.hpp:11
#define UNLIKELY(x)
Definition: Gsl.hpp:72
T signaling_NaN(T... args)
Defines function make_array.
An ElementId uniquely labels an Element. It is constructed from the BlockId of the Block to which the...
Definition: ElementId.hpp:36
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
WenoType
Possible types of the WENO limiter.
Definition: WenoType.hpp:19
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:29
A particular Side along a particular coordinate Axis.
Definition: Direction.hpp:23
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
double mean_value(const DataVector &f, const Mesh< Dim > &mesh) noexcept
Compute the mean value of a grid-function over a manifold. .
Definition: MeanValue.hpp:38
A spectral element with knowledge of its neighbors.
Definition: Element.hpp:29
Defines class Variables.
Interpolate a Variables from a Mesh onto a regular grid of points.
Definition: RegularGridInterpolant.hpp:46
Definition: DataBoxTag.hpp:29
Defines the class template Mesh.
Defines classes for Tensor.
void interpolate(gsl::not_null< Variables< TagsList > *> result, const Variables< TagsList > &vars) const noexcept
Interpolate Variables onto new mesh.
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:75
decltype(auto) for_each(const Container &c, UnaryFunction &&f)
Convenience wrapper around std::for_each, returns the result of std::for_each(begin(c), end(c), f).
Definition: Algorithm.hpp:241
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:436
Defines functions and classes from the GSL.
The inertial-coordinate size of an element along each of its logical directions.
Definition: SizeOfElement.hpp:53
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Defines tags related to domain quantities.
Defines function mean_value and mean_value_on_boundary.
Things relating to limiting.
Definition: HwenoModifiedSolution.cpp:109
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
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 class Element.