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  // Check that basis is LGL or LG
59  // Note that the SimpleWeno implementation should generalize well to other
60  // bases beyond LGL or LG, once the oscillation_indicator function has been
61  // generalized.
62  ASSERT(mesh.basis() == make_array<VolumeDim>(Spectral::Basis::Legendre),
63  "Unsupported basis: " << mesh);
64  ASSERT(mesh.quadrature() ==
65  make_array<VolumeDim>(Spectral::Quadrature::GaussLobatto) or
66  mesh.quadrature() ==
67  make_array<VolumeDim>(Spectral::Quadrature::Gauss),
68  "Unsupported quadrature: " << mesh);
69 
70  ASSERT(
71  modified_neighbor_solution_buffer->size() == neighbor_data.size(),
72  "modified_neighbor_solution_buffer->size() = "
73  << modified_neighbor_solution_buffer->size()
74  << "\nneighbor_data.size() = " << neighbor_data.size()
75  << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
76  "before calling simple_weno_impl.");
77 
78  // Compute the modified neighbor solutions.
79  // First extrapolate neighbor data onto local grid points, then shift the
80  // extrapolated data so its mean matches the local mean.
81  DataVector& component_to_limit = (*tensor)[tensor_storage_index];
82  const double local_mean = mean_value(component_to_limit, mesh);
83  for (const auto& neighbor_and_data : neighbor_data) {
84  const auto& neighbor = neighbor_and_data.first;
85  const auto& data = neighbor_and_data.second;
86 
87  if (interpolator_buffer->find(neighbor) == interpolator_buffer->end()) {
88  // No interpolator found => create one
89  const auto& direction = neighbor.first;
90  const auto& source_mesh = data.mesh;
91  const auto target_1d_logical_coords =
92  Weno_detail::local_grid_points_in_neighbor_logical_coords(
93  mesh, source_mesh, element, direction);
94  interpolator_buffer->insert(std::make_pair(
95  neighbor, intrp::RegularGrid<VolumeDim>(source_mesh, mesh,
96  target_1d_logical_coords)));
97  }
98 
99  // Avoid allocations by working directly in the preallocated buffer
100  DataVector& buffer = modified_neighbor_solution_buffer->at(neighbor);
101 
102  interpolator_buffer->at(neighbor).interpolate(
103  make_not_null(&buffer),
104  get<Tag>(data.volume_data)[tensor_storage_index]);
105  const double neighbor_mean = mean_value(buffer, mesh);
106  buffer += (local_mean - neighbor_mean);
107  }
108 
109  // Sum local and modified neighbor polynomials for the WENO reconstruction
110  Weno_detail::reconstruct_from_weighted_sum(
111  make_not_null(&component_to_limit), neighbor_linear_weight,
112  Weno_detail::DerivativeWeight::PowTwoEll, mesh,
113  *modified_neighbor_solution_buffer);
114 }
115 
116 // Compute the Simple WENO solution for one tensor
117 //
118 // This interface is intended for use in limiters that check the troubled-cell
119 // indicator for the whole cell, and apply the limiter to all fields.
120 //
121 // When calling `simple_weno_impl`,
122 // - `interpolator_buffer` may be empty
123 // - `modified_neighbor_solution_buffer` should contain one DataVector for each
124 // neighboring element (i.e. for each entry in `neighbor_data`)
125 template <typename Tag, size_t VolumeDim, typename PackagedData>
126 void simple_weno_impl(
131  interpolator_buffer,
135  modified_neighbor_solution_buffer,
137  const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
138  const Element<VolumeDim>& element,
139  const std::unordered_map<
142  neighbor_data) noexcept {
143  for (size_t tensor_storage_index = 0; tensor_storage_index < tensor->size();
144  ++tensor_storage_index) {
145  simple_weno_impl<Tag>(interpolator_buffer,
146  modified_neighbor_solution_buffer, tensor,
147  neighbor_linear_weight, tensor_storage_index, mesh,
148  element, neighbor_data);
149  }
150 }
151 
152 } // namespace Limiters::Weno_detail
utility
std::pair
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:50
ElementId.hpp
cstddef
array
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
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
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
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
string