SimpleWenoImpl.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 
21 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
22 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
23 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
24 #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
27 #include "Utilities/Gsl.hpp"
28 
29 namespace Limiters::Weno_detail {
30 
31 // Compute the Simple WENO solution for one tensor component
32 //
33 // This interface is intended for use in limiters that check the troubled-cell
34 // indicator independently for each tensor component. These limiters generally
35 // need to limit only a subset of the tensor components.
36 //
37 // When calling `simple_weno_impl`,
38 // - `interpolator_buffer` may be empty
39 // - `modified_neighbor_solution_buffer` should contain one DataVector for each
40 // neighboring element (i.e. for each entry in `neighbor_data`)
41 template <typename Tag, size_t VolumeDim, typename PackagedData>
42 void simple_weno_impl(
47  interpolator_buffer,
51  modified_neighbor_solution_buffer,
53  const double neighbor_linear_weight, const size_t tensor_storage_index,
54  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
55  const std::unordered_map<
58  neighbor_data) noexcept {
59  ASSERT(
60  modified_neighbor_solution_buffer->size() == neighbor_data.size(),
61  "modified_neighbor_solution_buffer->size() = "
62  << modified_neighbor_solution_buffer->size()
63  << "\nneighbor_data.size() = " << neighbor_data.size()
64  << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
65  "before calling simple_weno_impl.");
66 
67  // Compute the modified neighbor solutions.
68  // First extrapolate neighbor data onto local grid points, then shift the
69  // extrapolated data so its mean matches the local mean.
70  DataVector& component_to_limit = (*tensor)[tensor_storage_index];
71  const double local_mean = mean_value(component_to_limit, mesh);
72  for (const auto& neighbor_and_data : neighbor_data) {
73  const auto& neighbor = neighbor_and_data.first;
74  const auto& data = neighbor_and_data.second;
75 
76  if (interpolator_buffer->find(neighbor) == interpolator_buffer->end()) {
77  // No interpolator found => create one
78  const auto& direction = neighbor.first;
79  const auto& source_mesh = data.mesh;
80  const auto target_1d_logical_coords =
81  Weno_detail::local_grid_points_in_neighbor_logical_coords(
82  mesh, source_mesh, element, direction);
83  interpolator_buffer->insert(std::make_pair(
84  neighbor, intrp::RegularGrid<VolumeDim>(source_mesh, mesh,
85  target_1d_logical_coords)));
86  }
87 
88  // Avoid allocations by working directly in the preallocated buffer
89  DataVector& buffer = modified_neighbor_solution_buffer->at(neighbor);
90 
91  interpolator_buffer->at(neighbor).interpolate(
92  make_not_null(&buffer),
93  get<Tag>(data.volume_data)[tensor_storage_index]);
94  const double neighbor_mean = mean_value(buffer, mesh);
95  buffer += (local_mean - neighbor_mean);
96  }
97 
98  // Sum local and modified neighbor polynomials for the WENO reconstruction
99  Weno_detail::reconstruct_from_weighted_sum(
100  make_not_null(&component_to_limit), neighbor_linear_weight,
101  Weno_detail::DerivativeWeight::PowTwoEll, mesh,
102  *modified_neighbor_solution_buffer);
103 }
104 
105 // Compute the Simple WENO solution for one tensor
106 //
107 // This interface is intended for use in limiters that check the troubled-cell
108 // indicator for the whole cell, and apply the limiter to all fields.
109 //
110 // When calling `simple_weno_impl`,
111 // - `interpolator_buffer` may be empty
112 // - `modified_neighbor_solution_buffer` should contain one DataVector for each
113 // neighboring element (i.e. for each entry in `neighbor_data`)
114 template <typename Tag, size_t VolumeDim, typename PackagedData>
115 void simple_weno_impl(
120  interpolator_buffer,
124  modified_neighbor_solution_buffer,
126  const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
127  const Element<VolumeDim>& element,
128  const std::unordered_map<
131  neighbor_data) noexcept {
132  for (size_t tensor_storage_index = 0; tensor_storage_index < tensor->size();
133  ++tensor_storage_index) {
134  simple_weno_impl<Tag>(interpolator_buffer,
135  modified_neighbor_solution_buffer, tensor,
136  neighbor_linear_weight, tensor_storage_index, mesh,
137  element, neighbor_data);
138  }
139 }
140 
141 } // namespace Limiters::Weno_detail
DataBoxTag.hpp
utility
std::pair
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
ElementId.hpp
cstddef
array
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
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:51
Variables.hpp
Element.hpp
Mesh< VolumeDim >
limits
Gsl.hpp
Tensor.hpp
make_not_null
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,...
Definition: Gsl.hpp:880
Direction.hpp
Mesh.hpp
intrp::RegularGrid
Interpolate data from a Mesh onto a regular grid of points.
Definition: RegularGridInterpolant.hpp:46
unordered_map
std::data
T data(T... args)
type_traits
MeanValue.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string