HwenoImpl.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 <functional>
10 #include <limits>
11 #include <unordered_map>
12 #include <utility>
13 #include <vector>
14 
16 #include "DataStructures/DataVector.hpp"
17 #include "DataStructures/Index.hpp"
19 #include "DataStructures/Tags.hpp" // IWYU pragma: keep
20 #include "DataStructures/Variables.hpp" // IWYU pragma: keep
21 #include "Domain/Direction.hpp"
22 #include "Domain/DirectionMap.hpp"
23 #include "ErrorHandling/Assert.hpp"
24 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
25 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
26 #include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
27 #include "NumericalAlgorithms/LinearOperators/ApplyMatrices.hpp"
29 #include "Utilities/Algorithm.hpp"
30 #include "Utilities/EqualWithinRoundoff.hpp"
31 #include "Utilities/Gsl.hpp"
32 
33 /// \cond
34 template <size_t VolumeDim>
35 class Element;
36 template <size_t VolumeDim>
37 class ElementId;
38 template <size_t>
39 class Mesh;
40 /// \endcond
41 
42 namespace Limiters {
43 namespace Weno_detail {
44 
45 // Caching class that holds various precomputed terms used in the constrained-
46 // fit algebra on each element.
47 //
48 // The terms to precompute and cache depend on the configuration of neighbors to
49 // the element (how many neighbors, in which directions, of what resolution).
50 // Each instance of the caching class represents one particular neighbor
51 // configuration, and typical simulations will create many instances of
52 // the caching class: one instance for each neighbor-element configuration
53 // present in the computational domain.
54 //
55 // Because the current implementation of the HWENO fitting makes the simplifying
56 // restrictions of no h/p-refinement, the structure of the caching is also
57 // simplified. In particular, with no h/p-refinement, the allowable neighbor
58 // configurations satisfy,
59 // 1. the element has at most one neighbor per dimension, AND
60 // 2. the mesh on every neighbor is the same as on the element.
61 // With these restrictions, the number of independent configurations to be
62 // cached is greatly reduced. Because an element has 2*VolumeDim boundaries,
63 // it has 2^(2*VolumeDim) possible configurations of internal/external
64 // boundaries, and therefore there are 2^(2*VolumeDim) configurations to cache.
65 // Different elements with the same configuration of internal/external
66 // boundaries vs. direction can share the same caching-class instance.
67 //
68 // Each instance of the caching class holds several terms, some of which also
69 // depend on the neighbor configuration. The restriction of no h/p-refinement
70 // therefore simplifies each cache instance, as well as reducing the necessary
71 // number of instances.
72 //
73 // The most complicated term to cache is the A^{-1} matrix. The complexity
74 // arises because the matrix can take many values depending on (runtime) choices
75 // made for each individual HWNEO fit: which element is the primary neighbor,
76 // and which element(s) are the excluded neighbors.
77 // Fits with one excluded neighbor are overwhelmingly the most likely, and the
78 // cache is designed for this particular case. Fits with no excluded neighbors
79 // can arise in elements with only one neighbor (always the primary neighbor);
80 // these cases are also handled. However, the very rare case in which more than
81 // one neighbor is excluded is not handled by the cache; the caller must compute
82 // A^{-1} from scratch if this scenario arises.
83 template <size_t VolumeDim>
84 class ConstrainedFitCache {
85  public:
86  ConstrainedFitCache(const Element<VolumeDim>& element,
87  const Mesh<VolumeDim>& mesh) noexcept;
88 
89  // Valid calls must satisfy these constraints on directions_to_exclude:
90  // - the vector is empty, OR
91  // - the vector contains one element, and this element is a direction
92  // different from the direction to the primary neighbor.
93  // The very rare case where more than one neighbor is excluded from the fit is
94  // not cached. In this case, A^{-1} must be computed from scratch.
95  const Matrix& retrieve_inverse_a_matrix(
96  const Direction<VolumeDim>& primary_direction,
97  const std::vector<Direction<VolumeDim>>& directions_to_exclude) const
98  noexcept;
99 
101  DirectionMap<VolumeDim, Matrix> interpolation_matrices;
103  quadrature_weights_dot_interpolation_matrices;
104  // The many possible values of A^{-1} are stored in a map of maps. The outer
105  // map indexes over the primary neighbor, and the inner map indexes over the
106  // excluded neighbor.
107  // This data structure is not perfect: configurations with only one neighbor
108  // (always the primary neighbor) lead to no excluded neighbors, and there is
109  // no natural place to hold A^{-1} in the maps. For this case, we simply store
110  // the data in the normally-nonsensical slot where
111  // excluded_neighbor == primary_neighbor.
113 };
114 
115 // Return the appropriate cache for the given element and mesh.
116 template <size_t VolumeDim>
117 const ConstrainedFitCache<VolumeDim>& constrained_fit_cache(
118  const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh) noexcept;
119 
120 // Compute the inverse of the matrix A_st for the constrained fit.
121 // See the documentation of `hweno_modified_neighbor_solution`.
122 template <size_t VolumeDim>
123 Matrix inverse_a_matrix(
124  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
126  const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
128  quadrature_weights_dot_interpolation_matrices,
129  const Direction<VolumeDim>& primary_direction,
130  const std::vector<Direction<VolumeDim>>& directions_to_exclude) noexcept;
131 
132 // Not all secondary neighbors (i.e., neighbors that are not the primary) are
133 // included in the HWENO constrained fit. In particular, for the tensor
134 // component specified by `Tag` and `tensor_index`, the secondary neighbor whose
135 // mean is most different from the troubled cell's mean is excluded.
136 //
137 // Zhu2016 indicate that when multiple secondary neighbors share the property of
138 // being "most different" from the troubled cell, then all of these secondary
139 // neighbors should be excluded from the minimization. The probability of this
140 // occuring is exceedingly low, but we handle it anyway. We return the excluded
141 // secondary neighbors in a vector.
142 //
143 // Note that if there is only one neighbor, it is the primary neighbor, and so
144 // there are no secondary neighbors to exclude. We return an empty vector. This
145 // scenario can arise in various test cases, but is unlikely to arise in science
146 // cases.
147 template <typename Tag, size_t VolumeDim, typename Package>
149 secondary_neighbors_to_exclude_from_fit(
150  const double local_mean, const size_t tensor_index,
151  const std::unordered_map<
154  neighbor_data,
156  primary_neighbor) noexcept {
157  // Rare case: with only one neighbor, there is no secondary to exclude
158  if (UNLIKELY(neighbor_data.size() == 1)) {
159  return std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>{};
160  }
161 
162  // Identify element with maximum mean difference
163  const auto mean_difference = [&tensor_index, &local_mean ](
165  Package>& neighbor_and_data) noexcept {
166  return fabs(
167  get<::Tags::Mean<Tag>>(neighbor_and_data.second.means)[tensor_index] -
168  local_mean);
169  };
170 
171  double max_difference = std::numeric_limits<double>::lowest();
172  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>
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) {
177  continue;
178  }
179  const double difference = mean_difference(neighbor_and_data);
180  if (difference > max_difference) {
181  max_difference = difference;
182  neighbor_max_difference = neighbor;
183  }
184  }
185 
186  std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>
187  neighbors_to_exclude{{neighbor_max_difference}};
188 
189  // See if other elements share this maximum mean difference. This loop should
190  // only rarely find other neighbors with the same maximal mean difference to
191  // add to the vector, so it will usually not change the vector.
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) {
195  continue;
196  }
197  const double difference = mean_difference(neighbor_and_data);
198  if (UNLIKELY(equal_within_roundoff(difference, max_difference))) {
199  neighbors_to_exclude.push_back(neighbor);
200  }
201  }
202 
203  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
204  "Logical inconsistency: trying to exclude the primary neighbor.");
205 
206  return neighbors_to_exclude;
207 }
208 
209 // Compute the vector b_s for the constrained fit. For details, see the
210 // documentation of `hweno_modified_neighbor_solution` below.
211 template <typename Tag, size_t VolumeDim, typename Package>
212 DataVector b_vector(
213  const size_t tensor_index, const Mesh<VolumeDim>& mesh,
215  const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
217  quadrature_weights_dot_interpolation_matrices,
218  const std::unordered_map<
219  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
220  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
221  neighbor_data,
222  const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
223  primary_neighbor,
224  const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
225  neighbors_to_exclude) noexcept {
226  const size_t number_of_grid_points = mesh.number_of_grid_points();
227  DataVector b(number_of_grid_points, 0.);
228 
229  for (const auto& neighbor_and_data : neighbor_data) {
230  // Generally, neighbors_to_exclude contains just one neighbor to exclude,
231  // and the loop over neighbor_data will hit this (in 3D) roughly 1/6 times.
232  // So the condition is somewhat, but not overwhelmingly, unlikely:
233  if (UNLIKELY(alg::found(neighbors_to_exclude, neighbor_and_data.first))) {
234  continue;
235  }
236 
237  const auto& direction = neighbor_and_data.first.first;
238  ASSERT(interpolation_matrices.contains(direction),
239  "interpolation_matrices does not contain key: " << direction);
240  ASSERT(
241  quadrature_weights_dot_interpolation_matrices.contains(direction),
242  "quadrature_weights_dot_interpolation_matrices does not contain key: "
243  << direction);
244 
245  const auto& neighbor_mesh = mesh;
246  const auto& neighbor_quadrature_weights = quadrature_weights;
247  const auto& interpolation_matrix = interpolation_matrices.at(direction);
248  const auto& quadrature_weights_dot_interpolation_matrix =
249  quadrature_weights_dot_interpolation_matrices.at(direction);
250 
251  const auto& neighbor_tensor_component =
252  get<Tag>(neighbor_and_data.second.volume_data)[tensor_index];
253 
254  // Add terms from the primary neighbor
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] *
259  neighbor_quadrature_weights[r] * interpolation_matrix(r, s);
260  }
261  }
262  }
263  // Add terms from the secondary neighbors
264  else {
265  const double quadrature_weights_dot_u = [&]() noexcept {
266  double result = 0.;
267  for (size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
268  result +=
269  neighbor_tensor_component[r] * neighbor_quadrature_weights[r];
270  }
271  return result;
272  }
273  ();
274  b += quadrature_weights_dot_u *
275  quadrature_weights_dot_interpolation_matrix;
276  }
277  }
278  return b;
279 }
280 
281 // Solve the constrained fit problem that gives the HWENO modified solution,
282 // for one particular tensor component. For details, see documentation of
283 // `hweno_modified_neighbor_solution` below.
284 template <typename Tag, size_t VolumeDim, typename Package>
285 void solve_constrained_fit(
286  const gsl::not_null<DataVector*> constrained_fit_result,
287  const DataVector& u, const size_t tensor_index, const Mesh<VolumeDim>& mesh,
288  const Element<VolumeDim>& element,
289  const std::unordered_map<
290  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
291  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
292  neighbor_data,
293  const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
294  primary_neighbor,
295  const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
296  neighbors_to_exclude) noexcept {
297  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
298  "Logical inconsistency: trying to exclude the primary neighbor.");
299  ASSERT(not neighbors_to_exclude.empty() or neighbor_data.size() == 1,
300  "The HWENO constrained fit algorithm expects at least one neighbor \n"
301  "to exclude from the fit (unless if the element has a single \n"
302  "neighbor, which would automatically be the primary neighbor).");
303 
304  // Get the cache of linear algebra quantities for this element
305  const ConstrainedFitCache<VolumeDim>& cache =
306  constrained_fit_cache(element, mesh);
307 
308  // Because we don't support h-refinement, the direction is the only piece
309  // of the neighbor information that we actually need.
310  const Direction<VolumeDim> primary_direction = primary_neighbor.first;
312  directions_to_exclude = [&neighbors_to_exclude]() noexcept {
313  std::vector<Direction<VolumeDim>> result(neighbors_to_exclude.size());
314  for (size_t i = 0; i < result.size(); ++i) {
315  result[i] = neighbors_to_exclude[i].first;
316  }
317  return result;
318  }
319  ();
320 
321  const DataVector& w = cache.quadrature_weights;
322  const DirectionMap<VolumeDim, Matrix>& interp_matrices =
323  cache.interpolation_matrices;
324  const DirectionMap<VolumeDim, DataVector>& w_dot_interp_matrices =
325  cache.quadrature_weights_dot_interpolation_matrices;
326 
327  // Use cache if possible, or compute matrix if we are in edge case
328  const Matrix& inverse_a =
329  LIKELY(directions_to_exclude.size() < 2)
330  ? cache.retrieve_inverse_a_matrix(primary_direction,
331  directions_to_exclude)
332  : inverse_a_matrix(mesh, element, w, interp_matrices,
333  w_dot_interp_matrices, primary_direction,
334  directions_to_exclude);
335 
336  const DataVector b = b_vector<Tag>(tensor_index, mesh, w, interp_matrices,
337  w_dot_interp_matrices, neighbor_data,
338  primary_neighbor, neighbors_to_exclude);
339 
340  const size_t number_of_points = b.size();
341  const DataVector inverse_a_times_b = apply_matrices(
343  Index<1>(number_of_points));
344  const DataVector inverse_a_times_w = apply_matrices(
346  Index<1>(number_of_points));
347 
348  // Compute Lagrange multiplier:
349  // Note: we take w as an argument (instead of as a lambda capture), because
350  // some versions of Clang incorrectly warn about capturing w.
351  const double lagrange_multiplier =
352  [&number_of_points, &inverse_a_times_b, &inverse_a_times_w, &
353  u ](const DataVector& local_w) noexcept {
354  double numerator = 0.;
355  double denominator = 0.;
356  for (size_t s = 0; s < number_of_points; ++s) {
357  numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
358  denominator += local_w[s] * inverse_a_times_w[s];
359  }
360  return -numerator / denominator;
361  }
362  (w);
363 
364  // Compute solution:
365  *constrained_fit_result =
366  inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
367 }
368 
369 /*!
370  * \ingroup LimitersGroup
371  * \brief Compute the HWENO solution for one tensor
372  *
373  * The HWENO limiter reconstructs a new solution from the linear combination of
374  * the local DG solution and a "modified" solution from each neighbor element.
375  * This function computes the modified solution for a particular tensor and
376  * neighbor, following Section 3 of \cite Zhu2016.
377  *
378  * The modified solution associated with a particular neighbor (the "primary"
379  * neighbor) is obtained by solving a constrained fit over the local element,
380  * the primary neighbor, and the other ("secondary") neighbors of the local
381  * element. This fit seeks to minimize in a least-squared sense:
382  * 1. The distance between the modified solution and the original solution on
383  * the primary neighbor.
384  * 2. The distance between the cell average of the modified solution and the
385  * cell average of the original solution on each secondary neighbor. Note
386  * however that one secondary neighbor (or, rarely, several) is excluded from
387  * this minimization: the secondary neighbor(s) where the original solution
388  * has the most different cell average from the local element. This helps to
389  * prevent an outlier (e.g., near a shock) from biasing the fit.
390  *
391  * The constraint on the minimization is the following: the cell average of the
392  * modified solution on the local element must equal the cell average of the
393  * local element's original solution.
394  *
395  * Below we give the mathematical form of the constraints described above and
396  * show how these are translated into a numerical algorithm.
397  *
398  * Consider an element \f$I_0\f$ with neighbors \f$I_1, I_2, ...\f$. For a
399  * given tensor component \f$u\f$, the values on each of these elements are
400  * \f$u^{(0)}, u^{(1)}, u^{(2)}, ...\f$. Taking for the sake of example the
401  * primary neighbor to be \f$I_1\f$, the modified solution \f$\phi\f$ must
402  * minimize
403  *
404  * \f[
405  * \chi^2 = \int_{I_1} (\phi - u^{(1)})^2 dV
406  * + \sum_{\ell \in L} \left(
407  * \int_{I_{\ell}} ( \phi - u^{(\ell)} ) dV
408  * \right)^2,
409  * \f]
410  *
411  * subject to the constaint
412  *
413  * \f[
414  * C = \int_{I_0} ( \phi - u^{(0)} ) dV = 0.
415  * \f]
416  *
417  * where \f$\ell\f$ ranges over a subset \f$L\f$ of the secondary neighbors.
418  * \f$L\f$ excludes the one (or more) secondary neighbor(s) where the mean of
419  * \f$u\f$ is the most different from the mean of \f$u^{(0)}\f$. Typically, only
420  * one secondary neighbor is excluded, so \f$L\f$ contains two fewer neighbors
421  * than the total number of neighbors to the element. Note that in 1D, this
422  * implies that \f$L\f$ is the empty set; for each modified solution, one
423  * neighbor is the primary neighbor and the other is the excluded neighbor.
424  *
425  * The integrals are evaluated by quadrature. We denote the quadrature weights
426  * by \f$w_s\f$ and the values of some data \f$X\f$ at the quadrature
427  * nodes by \f$X_s\f$. We use subscripts \f$r,s\f$ to denote quadrature nodes
428  * on the neighbor and local elements, respectively. The minimization becomes
429  *
430  * \f[
431  * \chi^2 = \sum_r w^{(1)}_r ( \phi_r - u^{(1)}_r )^2
432  * + \sum_{\ell \in L} \left(
433  * \sum_r w^{(\ell)}_r ( \phi_r - u^{(\ell)}_r )
434  * \right)^2,
435  * \f]
436  *
437  * subject to the constraint
438  *
439  * \f[
440  * C = \sum_s w^{(0)}_s ( \phi_s - u^{(0)}_s ) = 0.
441  * \f]
442  *
443  * Note that \f$\phi\f$ is a function defined on the local element \f$I_0\f$,
444  * and so is fully represented by its values \f$\phi_s\f$ at the quadrature
445  * points on this element. When evaluating \f$\phi\f$ on element \f$I_{\ell}\f$,
446  * we obtain the function values \f$\phi_r\f$ by polynomial extrapolation,
447  * \f$\phi_r = \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s\f$, where
448  * \f$\mathcal{I}^{(\ell)}_{rs}\f$ is the interpolation/extrapolation matrix
449  * that interpolates data defined at grid points \f$x_s\f$ and evaluates it at
450  * grid points \f$x_r\f$ of \f$I_{\ell}\f$. Thus,
451  *
452  * \f[
453  * \chi^2 = \sum_r w^{(1)}_r
454  * \left(
455  * \sum_s \mathcal{I}^{(1)}_{rs} \phi_s - u^{(1)}_r
456  * \right)^2
457  * + \sum_{\ell \in L} \left(
458  * \sum_r w^{(\ell)}_r
459  * \left(
460  * \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s
461  * - u^{(\ell)}_r
462  * \right)
463  * \right)^2.
464  * \f]
465  *
466  * The solution to this optimization problem is found in the standard way,
467  * using a Lagrange multiplier \f$\lambda\f$ to impose the constraint:
468  *
469  * \f[
470  * 0 = \frac{d}{d \phi_s} \left( \chi^2 + \lambda C \right).
471  * \f]
472  *
473  * Working out the differentiation with respect to \f$\phi_s\f$ leads to the
474  * linear problem that must be inverted to obtain the solution,
475  *
476  * \f[
477  * 0 = A_{st} \phi_t - b_s - \lambda w^{(0)}_s,
478  * \f]
479  *
480  * where
481  *
482  * \f{align*}{
483  * A_{st} &= \sum_r \left( w^{(1)}_r
484  * \mathcal{I}^{(1)}_{rs} \mathcal{I}^{(1)}_{rt}
485  * \right)
486  * + \sum_{\ell \in L} \left(
487  * \sum_r \left(
488  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rt}
489  * \right)
490  * \cdot
491  * \sum_r \left(
492  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
493  * \right)
494  * \right)
495  * \\
496  * b_s &= \sum_r \left( w^{(1)}_r u^{(1)}_r \mathcal{I}^{(1)}_{rs} \right)
497  * + \sum_{\ell \in L} \left(
498  * \sum_r \left( w^{(\ell)}_r u^{(\ell)}_r \right)
499  * \cdot
500  * \sum_r \left(
501  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
502  * \right)
503  * \right).
504  * \f}
505  *
506  * Finally, the solution to the constrained fit is
507  *
508  * \f{align*}{
509  * \lambda &= - \frac{ \sum_s w^{(0)}_s
510  * \left( (A^{-1})_{st} b_t - u^{(0)}_s \right)
511  * }{ \sum_s w^{(0)}_s (A^{-1})_{st} w^{(0)}_t }
512  * \\
513  * \phi_s &= (A^{-1})_{st} ( b_t + \lambda w^{(0)}_t ).
514  * \f}
515  *
516  * Note that the matrix \f$A\f$ does not depend on the values of the tensor
517  * \f$u\f$, so its inverse \f$A^{-1}\f$ can be precomputed and stored.
518  *
519  * \warning
520  * Note also that the implementation currently does not support h- or
521  * p-refinement; this is checked by some assertions. The implementation is
522  * untested for grids where elements are curved, and it should not be expected
523  * to work in these cases.
524  */
525 template <typename Tag, size_t VolumeDim, typename PackagedData>
526 void hweno_impl(
528  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, DataVector,
529  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>*>
530  modified_neighbor_solution_buffer,
531  const gsl::not_null<db::item_type<Tag>*> tensor,
532  const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
533  const Element<VolumeDim>& element,
534  const std::unordered_map<
535  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, PackagedData,
536  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
537  neighbor_data) noexcept {
538  alg::for_each(neighbor_data, [&element, &
539  mesh ](const auto& neighbor_and_data) noexcept {
540  ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
541  element, neighbor_and_data.first.first),
542  "Found some amount of h-refinement; this is not supported");
543  ASSERT(neighbor_and_data.second.mesh == mesh,
544  "Found some amount of p-refinement; this is not supported");
545  });
546 
547  for (size_t tensor_index = 0; tensor_index < tensor->size(); ++tensor_index) {
548  const auto& tensor_component = (*tensor)[tensor_index];
549  for (const auto& neighbor_and_data : neighbor_data) {
550  const auto& primary_neighbor = neighbor_and_data.first;
551  const auto neighbors_to_exclude =
552  secondary_neighbors_to_exclude_from_fit<Tag>(
553  mean_value(tensor_component, mesh), tensor_index, neighbor_data,
554  primary_neighbor);
555 
556  DataVector& buffer =
557  modified_neighbor_solution_buffer->at(primary_neighbor);
558  solve_constrained_fit<Tag>(make_not_null(&buffer), tensor_component,
559  tensor_index, mesh, element, neighbor_data,
560  primary_neighbor, neighbors_to_exclude);
561  }
562 
563  // Sum local and modified neighbor polynomials for the WENO reconstruction
564  Weno_detail::reconstruct_from_weighted_sum(
565  make_not_null(&((*tensor)[tensor_index])), neighbor_linear_weight,
566  Weno_detail::DerivativeWeight::PowTwoEllOverEllFactorial, mesh,
567  *modified_neighbor_solution_buffer);
568  }
569 }
570 
571 } // namespace Weno_detail
572 } // namespace Limiters
Defines class template Direction.
Defines class Matrix.
#define UNLIKELY(x)
Definition: Gsl.hpp:72
An ElementId uniquely labels an Element. It is constructed from the BlockId of the Block to which the...
Definition: ElementId.hpp:36
An optimized map with Direction keys.
Definition: DirectionMap.hpp:15
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:50
void apply_matrices(const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:67
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Defines class template Index.
T lowest(T... args)
A particular Side along a particular coordinate Axis.
Definition: Direction.hpp:23
Matrix interpolation_matrix(const size_t num_points, const T &target_points) noexcept
Matrix used to interpolate to the target_points.
Definition: Spectral.cpp:319
bool found(const Container &c, const T &value)
Convenience wrapper around std::find, returns true if value is in c.
Definition: Algorithm.hpp:210
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
T fabs(T... args)
A spectral element with knowledge of its neighbors.
Definition: Element.hpp:29
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Defines class Variables.
Defines macro ASSERT.
Stores a collection of function values.
Definition: DataVector.hpp:42
An integer multi-index.
Definition: Index.hpp:31
decltype(auto) for_each(const Container &c, UnaryFunction &&f)
Convenience wrapper around std::for_each, returns the result of std::for_each(begin(c), end(c), f).
Definition: Algorithm.hpp:241
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:246
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
constexpr bool equal_within_roundoff(const double a, const double b, const double eps=std::numeric_limits< double >::epsilon() *100.0, const double scale=1.0) noexcept
Checks if two values a and b are equal within roundoff, by comparing abs(a - b) < (max(abs(a)...
Definition: EqualWithinRoundoff.hpp:15
Defines functions mean_value and mean_value_on_boundary.
const DataVector & quadrature_weights(const Mesh< 1 > &mesh) noexcept
Quadrature weights for a one-dimensional mesh.
#define LIKELY(x)
Definition: Gsl.hpp:66
Things relating to limiting.
Definition: HwenoImpl.cpp:109
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Given the Tag holding a Tensor<DataVector, ...>, swap the DataVector with a double.
Definition: Tags.hpp:23
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182