7 #include <boost/functional/hash.hpp>
15 #include "DataStructures/ApplyMatrices.hpp"
16 #include "DataStructures/DataVector.hpp"
19 #include "DataStructures/Tags.hpp"
22 #include "Domain/Structure/DirectionMap.hpp"
23 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
24 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
25 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
27 #include "Utilities/Algorithm.hpp"
28 #include "Utilities/EqualWithinRoundoff.hpp"
33 template <
size_t VolumeDim>
35 template <
size_t VolumeDim>
41 namespace Limiters::Weno_detail {
81 template <
size_t VolumeDim>
82 class ConstrainedFitCache {
93 const Matrix& retrieve_inverse_a_matrix(
101 quadrature_weights_dot_interpolation_matrices;
114 template <
size_t VolumeDim>
115 const ConstrainedFitCache<VolumeDim>& constrained_fit_cache(
120 template <
size_t VolumeDim>
126 quadrature_weights_dot_interpolation_matrices,
145 template <
typename Tag,
size_t VolumeDim,
typename Package>
147 secondary_neighbors_to_exclude_from_fit(
148 const double local_mean,
const size_t tensor_index,
154 primary_neighbor) noexcept {
156 if (
UNLIKELY(neighbor_data.size() == 1)) {
161 const auto mean_difference =
162 [&tensor_index, &local_mean](
164 Package>& neighbor_and_data) noexcept {
166 neighbor_and_data.second.means)[tensor_index] -
172 neighbor_max_difference{};
173 for (
const auto& neighbor_and_data : neighbor_data) {
174 const auto& neighbor = neighbor_and_data.first;
175 if (neighbor == primary_neighbor) {
178 const double difference = mean_difference(neighbor_and_data);
179 if (difference > max_difference) {
180 max_difference = difference;
181 neighbor_max_difference = neighbor;
186 neighbors_to_exclude{{neighbor_max_difference}};
191 for (
const auto& neighbor_and_data : neighbor_data) {
192 const auto& neighbor = neighbor_and_data.first;
193 if (neighbor == primary_neighbor or neighbor == neighbor_max_difference) {
196 const double difference = mean_difference(neighbor_and_data);
198 neighbors_to_exclude.push_back(neighbor);
203 "Logical inconsistency: trying to exclude the primary neighbor.");
205 return neighbors_to_exclude;
210 template <
typename Tag,
size_t VolumeDim,
typename Package>
216 quadrature_weights_dot_interpolation_matrices,
224 neighbors_to_exclude) noexcept {
225 const size_t number_of_grid_points = mesh.number_of_grid_points();
228 for (
const auto& neighbor_and_data : neighbor_data) {
236 const auto& direction = neighbor_and_data.first.first;
237 ASSERT(interpolation_matrices.contains(direction),
238 "interpolation_matrices does not contain key: " << direction);
240 quadrature_weights_dot_interpolation_matrices.contains(direction),
241 "quadrature_weights_dot_interpolation_matrices does not contain key: "
244 const auto& neighbor_mesh = mesh;
247 const auto& quadrature_weights_dot_interpolation_matrix =
248 quadrature_weights_dot_interpolation_matrices.at(direction);
250 const auto& neighbor_tensor_component =
251 get<Tag>(neighbor_and_data.second.volume_data)[tensor_index];
254 if (neighbor_and_data.first == primary_neighbor) {
255 for (
size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
256 for (
size_t s = 0; s < number_of_grid_points; ++s) {
257 b[s] += neighbor_tensor_component[r] *
264 const double quadrature_weights_dot_u = [&]() noexcept {
266 for (
size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
268 neighbor_tensor_component[r] * neighbor_quadrature_weights[r];
272 b += quadrature_weights_dot_u *
273 quadrature_weights_dot_interpolation_matrix;
282 template <
typename Tag,
size_t VolumeDim,
typename Package>
283 void solve_constrained_fit(
294 neighbors_to_exclude) noexcept {
296 "Logical inconsistency: trying to exclude the primary neighbor.");
297 ASSERT(not neighbors_to_exclude.empty() or neighbor_data.size() == 1,
298 "The HWENO constrained fit algorithm expects at least one neighbor \n"
299 "to exclude from the fit (unless if the element has a single \n"
300 "neighbor, which would automatically be the primary neighbor).");
303 const ConstrainedFitCache<VolumeDim>& cache =
304 constrained_fit_cache(element, mesh);
310 [&neighbors_to_exclude]() noexcept {
312 for (
size_t i = 0; i < result.size(); ++i) {
313 result[i] = neighbors_to_exclude[i].first;
318 const DataVector& w = cache.quadrature_weights;
320 cache.interpolation_matrices;
322 cache.quadrature_weights_dot_interpolation_matrices;
326 LIKELY(directions_to_exclude.size() < 2)
327 ? cache.retrieve_inverse_a_matrix(primary_direction,
328 directions_to_exclude)
329 : inverse_a_matrix(mesh, element, w, interp_matrices,
330 w_dot_interp_matrices, primary_direction,
331 directions_to_exclude);
333 const DataVector b = b_vector<Tag>(tensor_index, mesh, w, interp_matrices,
334 w_dot_interp_matrices, neighbor_data,
335 primary_neighbor, neighbors_to_exclude);
351 double numerator = 0.;
352 double denominator = 0.;
354 numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
355 denominator += local_w[s] * inverse_a_times_w[s];
357 return -numerator / denominator;
361 *constrained_fit_result =
362 inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
525 template <
typename Tag,
size_t VolumeDim,
typename PackagedData>
530 modified_neighbor_solution_buffer,
537 neighbor_data) noexcept {
539 modified_neighbor_solution_buffer->size() == neighbor_data.size(),
540 "modified_neighbor_solution_buffer->size() = "
541 << modified_neighbor_solution_buffer->size()
542 <<
"\nneighbor_data.size() = " << neighbor_data.size()
543 <<
"\nmodified_neighbor_solution_buffer was incorrectly initialized "
544 "before calling hweno_impl.");
547 neighbor_data, [&element, &mesh](
const auto& neighbor_and_data) noexcept {
548 ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
549 element, neighbor_and_data.first.first),
550 "Found some amount of h-refinement; this is not supported");
551 ASSERT(neighbor_and_data.second.mesh == mesh,
552 "Found some amount of p-refinement; this is not supported");
555 for (
size_t tensor_index = 0; tensor_index < tensor->size(); ++tensor_index) {
556 const auto& tensor_component = (*tensor)[tensor_index];
557 for (
const auto& neighbor_and_data : neighbor_data) {
558 const auto& primary_neighbor = neighbor_and_data.first;
559 const auto neighbors_to_exclude =
560 secondary_neighbors_to_exclude_from_fit<Tag>(
561 mean_value(tensor_component, mesh), tensor_index, neighbor_data,
565 modified_neighbor_solution_buffer->at(primary_neighbor);
566 solve_constrained_fit<Tag>(
make_not_null(&buffer), tensor_component,
567 tensor_index, mesh, element, neighbor_data,
568 primary_neighbor, neighbors_to_exclude);
572 Weno_detail::reconstruct_from_weighted_sum(
573 make_not_null(&((*tensor)[tensor_index])), neighbor_linear_weight,
574 Weno_detail::DerivativeWeight::PowTwoEllOverEllFactorial, mesh,
575 *modified_neighbor_solution_buffer);