WenoHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <unordered_map>
8 #include <utility>
9 
11 #include "DataStructures/Variables.hpp" // IWYU pragma: keep
12 #include "Domain/Direction.hpp" // IWYU pragma: keep
13 #include "Domain/ElementId.hpp" // IWYU pragma: keep
16 #include "Utilities/Gsl.hpp"
17 
18 // IWYU pragma: no_forward_declare Variables
19 
20 /// \cond
21 class DataVector;
22 template <size_t>
23 class Mesh;
24 
25 namespace boost {
26 template <class T>
27 struct hash;
28 } // namespace boost
29 /// \endcond
30 
31 namespace SlopeLimiters {
32 namespace Weno_detail {
33 
34 // Compute the WENO oscillation indicator (also called the smoothness indicator)
35 //
36 // The oscillation indicator measures the amount of variation in the input data,
37 // with larger indicator values corresponding to a larger amount of variation
38 // (either from large monotonic slopes or from oscillations).
39 //
40 // Implements an indicator similar to that of Eq. 23 of Dumbser2007, but with
41 // the necessary adaptations for use on square/cube grids. We favor this
42 // indicator because it is formulated in the reference coordinates, which we
43 // use for the WENO reconstruction, and because it lends itself to an efficient
44 // implementation.
45 //
46 // Where (this reference to be added to References.bib later, when it is cited
47 // from _rendered_ documentation):
48 // Dumbser2007:
49 // Dumbser, M and Kaeser, M
50 // Arbitrary high order non-oscillatory finite volume schemes on unstructured
51 // meshes for linear hyperbolic systems
52 // https://doi.org/10.1016/j.jcp.2006.06.043
53 template <size_t VolumeDim>
54 double oscillation_indicator(const DataVector& data,
55  const Mesh<VolumeDim>& mesh) noexcept;
56 
57 // Compute the unnormalized nonlinear WENO weights. This is a fairly standard
58 // choice of weights; see e.g., Eq. 3.9 of Zhong2013 or Eq. 3.6 of Zhu2016.
59 //
60 // Where (these references to be added to References.bib later, when they are
61 // cited from _rendered_ documentation):
62 // Zhong2013
63 // Zhong, X and Shu, C-W
64 // A simple weighted essentially non-oscillatory limiter for Runge-Kutta
65 // discontinuous Galerkin methods
66 // https://doi.org/10.1016/j.jcp.2012.08.028
67 //
68 // Zhu20016
69 // Zhu, J and Zhong, X and Shu, C and Qiu, J
70 // Runge-Kutta Discontinuous Galerkin Method with a Simple and Compact Hermite
71 // WENO Limiter
72 // https://doi.org/10.4208/cicp.070215.200715a
73 inline double unnormalized_nonlinear_weight(
74  const double linear_weight, const double oscillation_indicator) noexcept {
75  return linear_weight / square(1.e-6 + oscillation_indicator);
76 }
77 
78 // Compute the WENO weighted reconstruction of several polynomials, see e.g.,
79 // Eq. 4.3 of Zhu2016. This is fairly standard, though different references can
80 // differ in their choice of oscillation/smoothness indicator.
81 template <typename Tag, size_t VolumeDim, typename TagsList>
82 void reconstruct_from_weighted_sum(
83  const gsl::not_null<db::item_type<Tag>*> local_tensor,
84  const Mesh<VolumeDim>& mesh, const double neighbor_linear_weight,
85  const std::unordered_map<
87  Variables<TagsList>,
89  neighbor_vars) noexcept {
90  for (size_t tensor_storage_index = 0;
91  tensor_storage_index < local_tensor->size(); ++tensor_storage_index) {
92  auto& local_polynomial = (*local_tensor)[tensor_storage_index];
93 
94  // Store linear weights in `local_weights` and `neighbor_weights`
95  // These weights will have to be generalized for multiple neighbors per
96  // face for use with h-refinement and AMR.
97  double local_weight = 1.;
100  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>
101  neighbor_weights;
102  for (const auto& kv : neighbor_vars) {
103  local_weight -= neighbor_linear_weight;
104  neighbor_weights[kv.first] = neighbor_linear_weight;
105  }
106 
107  // Update `local_weights` and `neighbor_weights` to hold the unnormalized
108  // nonlinear weights.
109  local_weight = unnormalized_nonlinear_weight(
110  local_weight, oscillation_indicator(local_polynomial, mesh));
111  for (const auto& kv : neighbor_vars) {
112  const auto& key = kv.first;
113  const auto& neighbor_tensor_component =
114  get<Tag>(kv.second)[tensor_storage_index];
115  neighbor_weights[key] = unnormalized_nonlinear_weight(
116  neighbor_weights[key],
117  oscillation_indicator(neighbor_tensor_component, mesh));
118  }
119 
120  // Update `local_weights` and `neighbor_weights` to hold the normalized
121  // weights; these are the final weights of the WENO reconstruction.
122  double normalization = local_weight;
123  for (const auto& kv : neighbor_weights) {
124  normalization += kv.second;
125  }
126  local_weight /= normalization;
127  for (auto& kv : neighbor_weights) {
128  kv.second /= normalization;
129  }
130 
131  // Perform reconstruction, by superposition of local and neighbor
132  // polynomials.
133  const double mean = mean_value(local_polynomial, mesh);
134  local_polynomial = mean + local_weight * (local_polynomial - mean);
135  for (const auto& kv : neighbor_vars) {
136  const auto& key = kv.first;
137  const auto& neighbor_polynomial =
138  get<Tag>(kv.second)[tensor_storage_index];
139  const double neighbor_mean = mean_value(neighbor_polynomial, mesh);
140  local_polynomial +=
141  neighbor_weights.at(key) * (neighbor_polynomial - neighbor_mean);
142  }
143  }
144 }
145 
146 } // namespace Weno_detail
147 } // namespace SlopeLimiters
Defines class template Direction.
Definition: BoostMultiArray.hpp:11
An ElementId uniquely labels an Element. It is constructed from the BlockId of the Block to which the...
Definition: ElementId.hpp:36
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Defines class ElementId.
A particular Side along a particular coordinate Axis.
Definition: Direction.hpp:23
Things relating to slope limiting.
Definition: Krivodonova.hpp:52
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
Defines class Variables.
Define simple functions for constant expressions.
Stores a collection of function values.
Definition: DataVector.hpp:42
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 function mean_value and mean_value_on_boundary.
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12