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"
28 #include "Utilities/Algorithm.hpp"
29 #include "Utilities/EqualWithinRoundoff.hpp"
34 template <
size_t VolumeDim>
36 template <
size_t VolumeDim>
42 namespace Limiters::Weno_detail {
82 template <
size_t VolumeDim>
83 class ConstrainedFitCache {
94 const Matrix& retrieve_inverse_a_matrix(
102 quadrature_weights_dot_interpolation_matrices;
115 template <
size_t VolumeDim>
116 const ConstrainedFitCache<VolumeDim>& constrained_fit_cache(
121 template <
size_t VolumeDim>
127 quadrature_weights_dot_interpolation_matrices,
146 template <
typename Tag,
size_t VolumeDim,
typename Package>
148 secondary_neighbors_to_exclude_from_fit(
149 const double local_mean,
const size_t tensor_index,
155 primary_neighbor) noexcept {
157 if (
UNLIKELY(neighbor_data.size() == 1)) {
162 const auto mean_difference =
163 [&tensor_index, &local_mean](
165 Package>& neighbor_and_data) noexcept {
167 neighbor_and_data.second.means)[tensor_index] -
173 neighbor_max_difference{};
174 for (
const auto& neighbor_and_data : neighbor_data) {
175 const auto& neighbor = neighbor_and_data.first;
176 if (neighbor == primary_neighbor) {
179 const double difference = mean_difference(neighbor_and_data);
180 if (difference > max_difference) {
181 max_difference = difference;
182 neighbor_max_difference = neighbor;
187 neighbors_to_exclude{{neighbor_max_difference}};
192 for (
const auto& neighbor_and_data : neighbor_data) {
193 const auto& neighbor = neighbor_and_data.first;
194 if (neighbor == primary_neighbor or neighbor == neighbor_max_difference) {
197 const double difference = mean_difference(neighbor_and_data);
199 neighbors_to_exclude.push_back(neighbor);
204 "Logical inconsistency: trying to exclude the primary neighbor.");
206 return neighbors_to_exclude;
211 template <
typename Tag,
size_t VolumeDim,
typename Package>
217 quadrature_weights_dot_interpolation_matrices,
225 neighbors_to_exclude) noexcept {
226 const size_t number_of_grid_points = mesh.number_of_grid_points();
229 for (
const auto& neighbor_and_data : neighbor_data) {
237 const auto& direction = neighbor_and_data.first.first;
239 "interpolation_matrices does not contain key: " << direction);
241 quadrature_weights_dot_interpolation_matrices.
contains(direction),
242 "quadrature_weights_dot_interpolation_matrices does not contain key: "
245 const auto& neighbor_mesh = mesh;
248 const auto& quadrature_weights_dot_interpolation_matrix =
249 quadrature_weights_dot_interpolation_matrices.at(direction);
251 const auto& neighbor_tensor_component =
252 get<Tag>(neighbor_and_data.second.volume_data)[tensor_index];
255 if (neighbor_and_data.first == primary_neighbor) {
256 for (
size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
257 for (
size_t s = 0; s < number_of_grid_points; ++s) {
258 b[s] += neighbor_tensor_component[r] *
265 const double quadrature_weights_dot_u = [&]() noexcept {
267 for (
size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
269 neighbor_tensor_component[r] * neighbor_quadrature_weights[r];
273 b += quadrature_weights_dot_u *
274 quadrature_weights_dot_interpolation_matrix;
283 template <
typename Tag,
size_t VolumeDim,
typename Package>
284 void solve_constrained_fit(
295 neighbors_to_exclude) noexcept {
297 "Logical inconsistency: trying to exclude the primary neighbor.");
298 ASSERT(not neighbors_to_exclude.empty() or neighbor_data.size() == 1,
299 "The HWENO constrained fit algorithm expects at least one neighbor \n"
300 "to exclude from the fit (unless if the element has a single \n"
301 "neighbor, which would automatically be the primary neighbor).");
304 const ConstrainedFitCache<VolumeDim>&
cache =
305 constrained_fit_cache(mesh, element);
311 [&neighbors_to_exclude]() noexcept {
313 for (
size_t i = 0; i < result.size(); ++i) {
314 result[i] = neighbors_to_exclude[i].first;
321 cache.interpolation_matrices;
323 cache.quadrature_weights_dot_interpolation_matrices;
327 LIKELY(directions_to_exclude.size() < 2)
328 ?
cache.retrieve_inverse_a_matrix(primary_direction,
329 directions_to_exclude)
330 : inverse_a_matrix(mesh, element, w, interp_matrices,
331 w_dot_interp_matrices, primary_direction,
332 directions_to_exclude);
334 const DataVector b = b_vector<Tag>(tensor_index, mesh, w, interp_matrices,
335 w_dot_interp_matrices, neighbor_data,
336 primary_neighbor, neighbors_to_exclude);
352 double numerator = 0.;
353 double denominator = 0.;
355 numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
356 denominator += local_w[s] * inverse_a_times_w[s];
358 return -numerator / denominator;
362 *constrained_fit_result =
363 inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
526 template <
typename Tag,
size_t VolumeDim,
typename PackagedData>
531 modified_neighbor_solution_buffer,
538 neighbor_data) noexcept {
543 ASSERT(mesh.basis() == make_array<VolumeDim>(Spectral::Basis::Legendre),
544 "Unsupported basis: " << mesh);
545 ASSERT(mesh.quadrature() ==
546 make_array<VolumeDim>(Spectral::Quadrature::GaussLobatto) or
548 make_array<VolumeDim>(Spectral::Quadrature::Gauss),
549 "Unsupported quadrature: " << mesh);
552 modified_neighbor_solution_buffer->size() == neighbor_data.size(),
553 "modified_neighbor_solution_buffer->size() = "
554 << modified_neighbor_solution_buffer->size()
555 <<
"\nneighbor_data.size() = " << neighbor_data.size()
556 <<
"\nmodified_neighbor_solution_buffer was incorrectly initialized "
557 "before calling hweno_impl.");
560 neighbor_data, [&mesh, &element](
const auto& neighbor_and_data) noexcept {
561 ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
562 element, neighbor_and_data.first.first),
563 "Found some amount of h-refinement; this is not supported");
564 ASSERT(neighbor_and_data.second.mesh == mesh,
565 "Found some amount of p-refinement; this is not supported");
568 for (
size_t tensor_index = 0; tensor_index < tensor->size(); ++tensor_index) {
569 const auto& tensor_component = (*tensor)[tensor_index];
570 for (
const auto& neighbor_and_data : neighbor_data) {
571 const auto& primary_neighbor = neighbor_and_data.first;
572 const auto neighbors_to_exclude =
573 secondary_neighbors_to_exclude_from_fit<Tag>(
574 mean_value(tensor_component, mesh), tensor_index, neighbor_data,
578 modified_neighbor_solution_buffer->at(primary_neighbor);
579 solve_constrained_fit<Tag>(
make_not_null(&buffer), tensor_component,
580 tensor_index, mesh, element, neighbor_data,
581 primary_neighbor, neighbors_to_exclude);
585 Weno_detail::reconstruct_from_weighted_sum(
586 make_not_null(&((*tensor)[tensor_index])), neighbor_linear_weight,
587 Weno_detail::DerivativeWeight::PowTwoEllOverEllFactorial, mesh,
588 *modified_neighbor_solution_buffer);