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