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::Weno_detail {
43 
44 // Caching class that holds various precomputed terms used in the constrained-
45 // fit algebra on each element.
46 //
47 // The terms to precompute and cache depend on the configuration of neighbors to
48 // the element (how many neighbors, in which directions, of what resolution).
49 // Each instance of the caching class represents one particular neighbor
50 // configuration, and typical simulations will create many instances of
51 // the caching class: one instance for each neighbor-element configuration
52 // present in the computational domain.
53 //
54 // Because the current implementation of the HWENO fitting makes the simplifying
55 // restrictions of no h/p-refinement, the structure of the caching is also
56 // simplified. In particular, with no h/p-refinement, the allowable neighbor
57 // configurations satisfy,
58 // 1. the element has at most one neighbor per dimension, AND
59 // 2. the mesh on every neighbor is the same as on the element.
60 // With these restrictions, the number of independent configurations to be
61 // cached is greatly reduced. Because an element has 2*VolumeDim boundaries,
62 // it has 2^(2*VolumeDim) possible configurations of internal/external
63 // boundaries, and therefore there are 2^(2*VolumeDim) configurations to cache.
64 // Different elements with the same configuration of internal/external
65 // boundaries vs. direction can share the same caching-class instance.
66 //
67 // Each instance of the caching class holds several terms, some of which also
68 // depend on the neighbor configuration. The restriction of no h/p-refinement
69 // therefore simplifies each cache instance, as well as reducing the necessary
70 // number of instances.
71 //
72 // The most complicated term to cache is the A^{-1} matrix. The complexity
73 // arises because the matrix can take many values depending on (runtime) choices
74 // made for each individual HWNEO fit: which element is the primary neighbor,
75 // and which element(s) are the excluded neighbors.
76 // Fits with one excluded neighbor are overwhelmingly the most likely, and the
77 // cache is designed for this particular case. Fits with no excluded neighbors
78 // can arise in elements with only one neighbor (always the primary neighbor);
79 // these cases are also handled. However, the very rare case in which more than
80 // one neighbor is excluded is not handled by the cache; the caller must compute
81 // A^{-1} from scratch if this scenario arises.
82 template <size_t VolumeDim>
83 class ConstrainedFitCache {
84  public:
85  ConstrainedFitCache(const Element<VolumeDim>& element,
86  const Mesh<VolumeDim>& mesh) noexcept;
87 
88  // Valid calls must satisfy these constraints on directions_to_exclude:
89  // - the vector is empty, OR
90  // - the vector contains one element, and this element is a direction
91  // different from the direction to the primary neighbor.
92  // The very rare case where more than one neighbor is excluded from the fit is
93  // not cached. In this case, A^{-1} must be computed from scratch.
94  const Matrix& retrieve_inverse_a_matrix(
95  const Direction<VolumeDim>& primary_direction,
96  const std::vector<Direction<VolumeDim>>& directions_to_exclude) const
97  noexcept;
98 
100  DirectionMap<VolumeDim, Matrix> interpolation_matrices;
102  quadrature_weights_dot_interpolation_matrices;
103  // The many possible values of A^{-1} are stored in a map of maps. The outer
104  // map indexes over the primary neighbor, and the inner map indexes over the
105  // excluded neighbor.
106  // This data structure is not perfect: configurations with only one neighbor
107  // (always the primary neighbor) lead to no excluded neighbors, and there is
108  // no natural place to hold A^{-1} in the maps. For this case, we simply store
109  // the data in the normally-nonsensical slot where
110  // excluded_neighbor == primary_neighbor.
112 };
113 
114 // Return the appropriate cache for the given element and mesh.
115 template <size_t VolumeDim>
116 const ConstrainedFitCache<VolumeDim>& constrained_fit_cache(
117  const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh) noexcept;
118 
119 // Compute the inverse of the matrix A_st for the constrained fit.
120 // See the documentation of `hweno_modified_neighbor_solution`.
121 template <size_t VolumeDim>
122 Matrix inverse_a_matrix(
123  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
125  const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
127  quadrature_weights_dot_interpolation_matrices,
128  const Direction<VolumeDim>& primary_direction,
129  const std::vector<Direction<VolumeDim>>& directions_to_exclude) noexcept;
130 
131 // Not all secondary neighbors (i.e., neighbors that are not the primary) are
132 // included in the HWENO constrained fit. In particular, for the tensor
133 // component specified by `Tag` and `tensor_index`, the secondary neighbor whose
134 // mean is most different from the troubled cell's mean is excluded.
135 //
136 // Zhu2016 indicate that when multiple secondary neighbors share the property of
137 // being "most different" from the troubled cell, then all of these secondary
138 // neighbors should be excluded from the minimization. The probability of this
139 // occuring is exceedingly low, but we handle it anyway. We return the excluded
140 // secondary neighbors in a vector.
141 //
142 // Note that if there is only one neighbor, it is the primary neighbor, and so
143 // there are no secondary neighbors to exclude. We return an empty vector. This
144 // scenario can arise in various test cases, but is unlikely to arise in science
145 // cases.
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,
150  const std::unordered_map<
153  neighbor_data,
155  primary_neighbor) noexcept {
156  // Rare case: with only one neighbor, there is no secondary to exclude
157  if (UNLIKELY(neighbor_data.size() == 1)) {
159  }
160 
161  // Identify element with maximum mean difference
162  const auto mean_difference =
163  [&tensor_index, &local_mean](
165  Package>& neighbor_and_data) noexcept {
166  return fabs(get<::Tags::Mean<Tag>>(
167  neighbor_and_data.second.means)[tensor_index] -
168  local_mean);
169  };
170 
171  double max_difference = std::numeric_limits<double>::lowest();
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 
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<
221  neighbor_data,
223  primary_neighbor,
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  b += quadrature_weights_dot_u *
274  quadrature_weights_dot_interpolation_matrix;
275  }
276  }
277  return b;
278 }
279 
280 // Solve the constrained fit problem that gives the HWENO modified solution,
281 // for one particular tensor component. For details, see documentation of
282 // `hweno_modified_neighbor_solution` below.
283 template <typename Tag, size_t VolumeDim, typename Package>
284 void solve_constrained_fit(
285  const gsl::not_null<DataVector*> constrained_fit_result,
286  const DataVector& u, const size_t tensor_index, const Mesh<VolumeDim>& mesh,
287  const Element<VolumeDim>& element,
288  const std::unordered_map<
291  neighbor_data,
293  primary_neighbor,
295  neighbors_to_exclude) noexcept {
296  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
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).");
302 
303  // Get the cache of linear algebra quantities for this element
304  const ConstrainedFitCache<VolumeDim>& cache =
305  constrained_fit_cache(element, mesh);
306 
307  // Because we don't support h-refinement, the direction is the only piece
308  // of the neighbor information that we actually need.
309  const Direction<VolumeDim> primary_direction = primary_neighbor.first;
310  const std::vector<Direction<VolumeDim>> directions_to_exclude =
311  [&neighbors_to_exclude]() noexcept {
312  std::vector<Direction<VolumeDim>> result(neighbors_to_exclude.size());
313  for (size_t i = 0; i < result.size(); ++i) {
314  result[i] = neighbors_to_exclude[i].first;
315  }
316  return result;
317  }();
318 
319  const DataVector& w = cache.quadrature_weights;
320  const DirectionMap<VolumeDim, Matrix>& interp_matrices =
321  cache.interpolation_matrices;
322  const DirectionMap<VolumeDim, DataVector>& w_dot_interp_matrices =
323  cache.quadrature_weights_dot_interpolation_matrices;
324 
325  // Use cache if possible, or compute matrix if we are in edge case
326  const Matrix& inverse_a =
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);
333 
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);
337 
338  const size_t number_of_points = b.size();
339  const DataVector inverse_a_times_b = apply_matrices(
341  Index<1>(number_of_points));
342  const DataVector inverse_a_times_w = apply_matrices(
344  Index<1>(number_of_points));
345 
346  // Compute Lagrange multiplier:
347  // Note: we take w as an argument (instead of as a lambda capture), because
348  // some versions of Clang incorrectly warn about capturing w.
349  const double lagrange_multiplier = [&number_of_points, &inverse_a_times_b,
350  &inverse_a_times_w,
351  &u](const DataVector& local_w) noexcept {
352  double numerator = 0.;
353  double denominator = 0.;
354  for (size_t s = 0; s < number_of_points; ++s) {
355  numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
356  denominator += local_w[s] * inverse_a_times_w[s];
357  }
358  return -numerator / denominator;
359  }(w);
360 
361  // Compute solution:
362  *constrained_fit_result =
363  inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
364 }
365 
366 /*!
367  * \ingroup LimitersGroup
368  * \brief Compute the HWENO solution for one tensor
369  *
370  * The HWENO limiter reconstructs a new solution from the linear combination of
371  * the local DG solution and a "modified" solution from each neighbor element.
372  * This function computes the modified solution for a particular tensor and
373  * neighbor, following Section 3 of \cite Zhu2016.
374  *
375  * The modified solution associated with a particular neighbor (the "primary"
376  * neighbor) is obtained by solving a constrained fit over the local element,
377  * the primary neighbor, and the other ("secondary") neighbors of the local
378  * element. This fit seeks to minimize in a least-squared sense:
379  * 1. The distance between the modified solution and the original solution on
380  * the primary neighbor.
381  * 2. The distance between the cell average of the modified solution and the
382  * cell average of the original solution on each secondary neighbor. Note
383  * however that one secondary neighbor (or, rarely, several) is excluded from
384  * this minimization: the secondary neighbor(s) where the original solution
385  * has the most different cell average from the local element. This helps to
386  * prevent an outlier (e.g., near a shock) from biasing the fit.
387  *
388  * The constraint on the minimization is the following: the cell average of the
389  * modified solution on the local element must equal the cell average of the
390  * local element's original solution.
391  *
392  * Below we give the mathematical form of the constraints described above and
393  * show how these are translated into a numerical algorithm.
394  *
395  * Consider an element \f$I_0\f$ with neighbors \f$I_1, I_2, ...\f$. For a
396  * given tensor component \f$u\f$, the values on each of these elements are
397  * \f$u^{(0)}, u^{(1)}, u^{(2)}, ...\f$. Taking for the sake of example the
398  * primary neighbor to be \f$I_1\f$, the modified solution \f$\phi\f$ must
399  * minimize
400  *
401  * \f[
402  * \chi^2 = \int_{I_1} (\phi - u^{(1)})^2 dV
403  * + \sum_{\ell \in L} \left(
404  * \int_{I_{\ell}} ( \phi - u^{(\ell)} ) dV
405  * \right)^2,
406  * \f]
407  *
408  * subject to the constaint
409  *
410  * \f[
411  * C = \int_{I_0} ( \phi - u^{(0)} ) dV = 0.
412  * \f]
413  *
414  * where \f$\ell\f$ ranges over a subset \f$L\f$ of the secondary neighbors.
415  * \f$L\f$ excludes the one (or more) secondary neighbor(s) where the mean of
416  * \f$u\f$ is the most different from the mean of \f$u^{(0)}\f$. Typically, only
417  * one secondary neighbor is excluded, so \f$L\f$ contains two fewer neighbors
418  * than the total number of neighbors to the element. Note that in 1D, this
419  * implies that \f$L\f$ is the empty set; for each modified solution, one
420  * neighbor is the primary neighbor and the other is the excluded neighbor.
421  *
422  * The integrals are evaluated by quadrature. We denote the quadrature weights
423  * by \f$w_s\f$ and the values of some data \f$X\f$ at the quadrature
424  * nodes by \f$X_s\f$. We use subscripts \f$r,s\f$ to denote quadrature nodes
425  * on the neighbor and local elements, respectively. The minimization becomes
426  *
427  * \f[
428  * \chi^2 = \sum_r w^{(1)}_r ( \phi_r - u^{(1)}_r )^2
429  * + \sum_{\ell \in L} \left(
430  * \sum_r w^{(\ell)}_r ( \phi_r - u^{(\ell)}_r )
431  * \right)^2,
432  * \f]
433  *
434  * subject to the constraint
435  *
436  * \f[
437  * C = \sum_s w^{(0)}_s ( \phi_s - u^{(0)}_s ) = 0.
438  * \f]
439  *
440  * Note that \f$\phi\f$ is a function defined on the local element \f$I_0\f$,
441  * and so is fully represented by its values \f$\phi_s\f$ at the quadrature
442  * points on this element. When evaluating \f$\phi\f$ on element \f$I_{\ell}\f$,
443  * we obtain the function values \f$\phi_r\f$ by polynomial extrapolation,
444  * \f$\phi_r = \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s\f$, where
445  * \f$\mathcal{I}^{(\ell)}_{rs}\f$ is the interpolation/extrapolation matrix
446  * that interpolates data defined at grid points \f$x_s\f$ and evaluates it at
447  * grid points \f$x_r\f$ of \f$I_{\ell}\f$. Thus,
448  *
449  * \f[
450  * \chi^2 = \sum_r w^{(1)}_r
451  * \left(
452  * \sum_s \mathcal{I}^{(1)}_{rs} \phi_s - u^{(1)}_r
453  * \right)^2
454  * + \sum_{\ell \in L} \left(
455  * \sum_r w^{(\ell)}_r
456  * \left(
457  * \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s
458  * - u^{(\ell)}_r
459  * \right)
460  * \right)^2.
461  * \f]
462  *
463  * The solution to this optimization problem is found in the standard way,
464  * using a Lagrange multiplier \f$\lambda\f$ to impose the constraint:
465  *
466  * \f[
467  * 0 = \frac{d}{d \phi_s} \left( \chi^2 + \lambda C \right).
468  * \f]
469  *
470  * Working out the differentiation with respect to \f$\phi_s\f$ leads to the
471  * linear problem that must be inverted to obtain the solution,
472  *
473  * \f[
474  * 0 = A_{st} \phi_t - b_s - \lambda w^{(0)}_s,
475  * \f]
476  *
477  * where
478  *
479  * \f{align*}{
480  * A_{st} &= \sum_r \left( w^{(1)}_r
481  * \mathcal{I}^{(1)}_{rs} \mathcal{I}^{(1)}_{rt}
482  * \right)
483  * + \sum_{\ell \in L} \left(
484  * \sum_r \left(
485  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rt}
486  * \right)
487  * \cdot
488  * \sum_r \left(
489  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
490  * \right)
491  * \right)
492  * \\
493  * b_s &= \sum_r \left( w^{(1)}_r u^{(1)}_r \mathcal{I}^{(1)}_{rs} \right)
494  * + \sum_{\ell \in L} \left(
495  * \sum_r \left( w^{(\ell)}_r u^{(\ell)}_r \right)
496  * \cdot
497  * \sum_r \left(
498  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
499  * \right)
500  * \right).
501  * \f}
502  *
503  * Finally, the solution to the constrained fit is
504  *
505  * \f{align*}{
506  * \lambda &= - \frac{ \sum_s w^{(0)}_s
507  * \left( (A^{-1})_{st} b_t - u^{(0)}_s \right)
508  * }{ \sum_s w^{(0)}_s (A^{-1})_{st} w^{(0)}_t }
509  * \\
510  * \phi_s &= (A^{-1})_{st} ( b_t + \lambda w^{(0)}_t ).
511  * \f}
512  *
513  * Note that the matrix \f$A\f$ does not depend on the values of the tensor
514  * \f$u\f$, so its inverse \f$A^{-1}\f$ can be precomputed and stored.
515  *
516  * \warning
517  * Note also that the implementation currently does not support h- or
518  * p-refinement; this is checked by some assertions. The implementation is
519  * untested for grids where elements are curved, and it should not be expected
520  * to work in these cases.
521  *
522  * When calling `hweno_impl`, `modified_neighbor_solution_buffer` should contain
523  * one DataVector for each neighboring element (i.e. for each entry in
524  * `neighbor_data`).
525  */
526 template <typename Tag, size_t VolumeDim, typename PackagedData>
527 void hweno_impl(
531  modified_neighbor_solution_buffer,
532  const gsl::not_null<db::item_type<Tag>*> tensor,
533  const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
534  const Element<VolumeDim>& element,
535  const std::unordered_map<
538  neighbor_data) noexcept {
539  ASSERT(
540  modified_neighbor_solution_buffer->size() == neighbor_data.size(),
541  "modified_neighbor_solution_buffer->size() = "
542  << modified_neighbor_solution_buffer->size()
543  << "\nneighbor_data.size() = " << neighbor_data.size()
544  << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
545  "before calling hweno_impl.");
546 
548  neighbor_data, [&element, &mesh](const auto& neighbor_and_data) noexcept {
549  ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
550  element, neighbor_and_data.first.first),
551  "Found some amount of h-refinement; this is not supported");
552  ASSERT(neighbor_and_data.second.mesh == mesh,
553  "Found some amount of p-refinement; this is not supported");
554  });
555 
556  for (size_t tensor_index = 0; tensor_index < tensor->size(); ++tensor_index) {
557  const auto& tensor_component = (*tensor)[tensor_index];
558  for (const auto& neighbor_and_data : neighbor_data) {
559  const auto& primary_neighbor = neighbor_and_data.first;
560  const auto neighbors_to_exclude =
561  secondary_neighbors_to_exclude_from_fit<Tag>(
562  mean_value(tensor_component, mesh), tensor_index, neighbor_data,
563  primary_neighbor);
564 
565  DataVector& buffer =
566  modified_neighbor_solution_buffer->at(primary_neighbor);
567  solve_constrained_fit<Tag>(make_not_null(&buffer), tensor_component,
568  tensor_index, mesh, element, neighbor_data,
569  primary_neighbor, neighbors_to_exclude);
570  }
571 
572  // Sum local and modified neighbor polynomials for the WENO reconstruction
573  Weno_detail::reconstruct_from_weighted_sum(
574  make_not_null(&((*tensor)[tensor_index])), neighbor_linear_weight,
575  Weno_detail::DerivativeWeight::PowTwoEllOverEllFactorial, mesh,
576  *modified_neighbor_solution_buffer);
577  }
578 }
579 
580 } // namespace Limiters::Weno_detail
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
apply_matrices
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
DataBoxTag.hpp
utility
std::fabs
T fabs(T... args)
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:689
functional
std::pair
vector
Tags::Mean
Given the Tag holding a Tensor<DataVector, ...>, swap the DataVector with a double.
Definition: Tags.hpp:23
Index
Definition: Index.hpp:31
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
std::numeric_limits::lowest
T lowest(T... args)
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
std::reference_wrapper
cstddef
Assert.hpp
alg::found
bool found(const Container &c, const T &value)
Convenience wrapper around std::find, returns true if value is in c.
Definition: Algorithm.hpp:210
array
db::item_type
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
Index.hpp
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
alg::for_each
decltype(auto) for_each(const Container &c, UnaryFunction &&f)
Convenience wrapper around std::for_each, returns the result of std::for_each(begin(c),...
Definition: Algorithm.hpp:241
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:50
DirectionMap
Definition: DirectionMap.hpp:15
limits
Gsl.hpp
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
Matrix.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
equal_within_roundoff
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
Definition: EqualWithinRoundoff.hpp:15
Spectral::quadrature_weights
const DataVector & quadrature_weights(const Mesh< 1 > &mesh) noexcept
Quadrature weights for a one-dimensional mesh.
Spectral::interpolation_matrix
Matrix interpolation_matrix(const size_t num_points, const T &target_points) noexcept
Matrix used to interpolate to the target_points.
Definition: Spectral.cpp:325
unordered_map
MeanValue.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183