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