HwenoModifiedSolution.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 "NumericalAlgorithms/LinearOperators/ApplyMatrices.hpp"
27 #include "Utilities/Algorithm.hpp"
28 #include "Utilities/EqualWithinRoundoff.hpp"
29 #include "Utilities/Gsl.hpp"
30 
31 /// \cond
32 template <size_t VolumeDim>
33 class Element;
34 template <size_t VolumeDim>
35 class ElementId;
36 template <size_t>
37 class Mesh;
38 /// \endcond
39 
40 namespace Limiters {
41 namespace Hweno_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 Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
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)) {
157  return std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>{};
158  }
159 
160  // Identify element with maximum mean difference
161  const auto mean_difference = [&tensor_index, &local_mean ](
163  Package>& neighbor_and_data) noexcept {
164  return fabs(
165  get<::Tags::Mean<Tag>>(neighbor_and_data.second.means)[tensor_index] -
166  local_mean);
167  };
168 
169  double max_difference = std::numeric_limits<double>::lowest();
170  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>
171  neighbor_max_difference{};
172  for (const auto& neighbor_and_data : neighbor_data) {
173  const auto& neighbor = neighbor_and_data.first;
174  if (neighbor == primary_neighbor) {
175  continue;
176  }
177  const double difference = mean_difference(neighbor_and_data);
178  if (difference > max_difference) {
179  max_difference = difference;
180  neighbor_max_difference = neighbor;
181  }
182  }
183 
184  std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>
185  neighbors_to_exclude{{neighbor_max_difference}};
186 
187  // See if other elements share this maximum mean difference. This loop should
188  // only rarely find other neighbors with the same maximal mean difference to
189  // add to the vector, so it will usually not change the vector.
190  for (const auto& neighbor_and_data : neighbor_data) {
191  const auto& neighbor = neighbor_and_data.first;
192  if (neighbor == primary_neighbor or neighbor == neighbor_max_difference) {
193  continue;
194  }
195  const double difference = mean_difference(neighbor_and_data);
196  if (UNLIKELY(equal_within_roundoff(difference, max_difference))) {
197  neighbors_to_exclude.push_back(neighbor);
198  }
199  }
200 
201  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
202  "Logical inconsistency: trying to exclude the primary neighbor.");
203 
204  return neighbors_to_exclude;
205 }
206 
207 // Compute the vector b_s for the constrained fit. For details, see the
208 // documentation of `hweno_modified_neighbor_solution` below.
209 template <typename Tag, size_t VolumeDim, typename Package>
210 DataVector b_vector(
211  const Mesh<VolumeDim>& mesh, const size_t tensor_index,
213  const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
215  quadrature_weights_dot_interpolation_matrices,
216  const std::unordered_map<
217  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
218  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
219  neighbor_data,
220  const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
221  primary_neighbor,
222  const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
223  neighbors_to_exclude) noexcept {
224  const size_t number_of_grid_points = mesh.number_of_grid_points();
225  DataVector b(number_of_grid_points, 0.);
226 
227  for (const auto& neighbor_and_data : neighbor_data) {
228  // Generally, neighbors_to_exclude contains just one neighbor to exclude,
229  // and the loop over neighbor_data will hit this (in 3D) roughly 1/6 times.
230  // So the condition is somewhat, but not overwhelmingly, unlikely:
231  if (UNLIKELY(alg::found(neighbors_to_exclude, neighbor_and_data.first))) {
232  continue;
233  }
234 
235  const auto& direction = neighbor_and_data.first.first;
236  ASSERT(interpolation_matrices.contains(direction),
237  "interpolation_matrices does not contain key: " << direction);
238  ASSERT(
239  quadrature_weights_dot_interpolation_matrices.contains(direction),
240  "quadrature_weights_dot_interpolation_matrices does not contain key: "
241  << direction);
242 
243  const auto& neighbor_mesh = mesh;
244  const auto& neighbor_quadrature_weights = quadrature_weights;
245  const auto& interpolation_matrix = interpolation_matrices.at(direction);
246  const auto& quadrature_weights_dot_interpolation_matrix =
247  quadrature_weights_dot_interpolation_matrices.at(direction);
248 
249  const auto& neighbor_tensor_component =
250  get<Tag>(neighbor_and_data.second.volume_data)[tensor_index];
251 
252  // Add terms from the primary neighbor
253  if (neighbor_and_data.first == primary_neighbor) {
254  for (size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
255  for (size_t s = 0; s < number_of_grid_points; ++s) {
256  b[s] += neighbor_tensor_component[r] *
257  neighbor_quadrature_weights[r] * interpolation_matrix(r, s);
258  }
259  }
260  }
261  // Add terms from the secondary neighbors
262  else {
263  const double quadrature_weights_dot_u = [&]() noexcept {
264  double result = 0.;
265  for (size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
266  result +=
267  neighbor_tensor_component[r] * neighbor_quadrature_weights[r];
268  }
269  return result;
270  }
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,
286  const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
287  const std::unordered_map<
288  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
289  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
290  neighbor_data,
291  const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
292  primary_neighbor,
293  const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
294  neighbors_to_exclude) noexcept {
295  // Get the cache of linear algebra quantities for this element
296  const ConstrainedFitCache<VolumeDim>& cache =
297  constrained_fit_cache(element, mesh);
298 
299  // Because we don't support h-refinement, the direction is the only piece
300  // of the neighbor information that we actually need.
301  const Direction<VolumeDim> primary_direction = primary_neighbor.first;
303  directions_to_exclude = [&neighbors_to_exclude]() noexcept {
304  std::vector<Direction<VolumeDim>> result(neighbors_to_exclude.size());
305  for (size_t i = 0; i < result.size(); ++i) {
306  result[i] = neighbors_to_exclude[i].first;
307  }
308  return result;
309  }
310  ();
311 
312  const DataVector& w = cache.quadrature_weights;
313  const DirectionMap<VolumeDim, Matrix>& interp_matrices =
314  cache.interpolation_matrices;
315  const DirectionMap<VolumeDim, DataVector>& w_dot_interp_matrices =
316  cache.quadrature_weights_dot_interpolation_matrices;
317 
318  // Use cache if possible, or compute matrix if we are in edge case
319  const Matrix& inverse_a =
320  LIKELY(directions_to_exclude.size() < 2)
321  ? cache.retrieve_inverse_a_matrix(primary_direction,
322  directions_to_exclude)
323  : Hweno_detail::inverse_a_matrix(
324  element, mesh, w, interp_matrices, w_dot_interp_matrices,
325  primary_direction, directions_to_exclude);
326 
327  const DataVector b = b_vector<Tag>(mesh, tensor_index, w, interp_matrices,
328  w_dot_interp_matrices, neighbor_data,
329  primary_neighbor, neighbors_to_exclude);
330 
331  const size_t number_of_points = b.size();
332  const DataVector inverse_a_times_b = apply_matrices(
334  Index<1>(number_of_points));
335  const DataVector inverse_a_times_w = apply_matrices(
337  Index<1>(number_of_points));
338 
339  // Compute Lagrange multiplier:
340  // Note: we take w as an argument (instead of as a lambda capture), because
341  // some versions of Clang incorrectly warn about capturing w.
342  const double lagrange_multiplier =
343  [&number_of_points, &inverse_a_times_b, &inverse_a_times_w, &
344  u ](const DataVector& local_w) noexcept {
345  double numerator = 0.;
346  double denominator = 0.;
347  for (size_t s = 0; s < number_of_points; ++s) {
348  numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
349  denominator += local_w[s] * inverse_a_times_w[s];
350  }
351  return -numerator / denominator;
352  }
353  (w);
354 
355  // Compute solution:
356  *constrained_fit_result =
357  inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
358 }
359 
360 } // namespace Hweno_detail
361 
362 template <size_t VolumeDim, typename Package>
363 using LimiterNeighborData = std::unordered_map<
364  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
365  boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>;
366 
367 /*!
368  * \ingroup LimitersGroup
369  * \brief Compute the HWENO modified solution for a particular tensor
370  * from a particular neighbor element
371  *
372  * The HWENO limiter reconstructs a new solution from the linear combination of
373  * the local DG solution and a "modified" solution from each neighbor element.
374  * This function computes the modified solution for a particular tensor and
375  * neighbor, following Section 3 of \cite Zhu2016.
376  *
377  * The modified solution associated with a particular neighbor (the "primary"
378  * neighbor) is obtained by solving a constrained fit over the local element,
379  * the primary neighbor, and the other ("secondary") neighbors of the local
380  * element. This fit seeks to minimize in a least-squared sense:
381  * 1. The distance between the modified solution and the original solution on
382  * the primary neighbor.
383  * 2. The distance between the cell average of the modified solution and the
384  * cell average of the original solution on each secondary neighbor. Note
385  * however that one secondary neighbor (or, rarely, several) is excluded from
386  * this minimization: the secondary neighbor(s) where the original solution
387  * has the most different cell average from the local element. This helps to
388  * prevent an outlier (e.g., near a shock) from biasing the fit.
389  *
390  * The constraint on the minimization is the following: the cell average of the
391  * modified solution on the local element must equal the cell average of the
392  * local element's original solution.
393  *
394  * Below we give the mathematical form of the constraints described above and
395  * show how these are translated into a numerical algorithm.
396  *
397  * Consider an element \f$I_0\f$ with neighbors \f$I_1, I_2, ...\f$. For a
398  * given tensor component \f$u\f$, the values on each of these elements are
399  * \f$u^{(0)}, u^{(1)}, u^{(2)}, ...\f$. Taking for the sake of example the
400  * primary neighbor to be \f$I_1\f$, the modified solution \f$\phi\f$ must
401  * minimize
402  *
403  * \f[
404  * \chi^2 = \int_{I_1} (\phi - u^{(1)})^2 dV
405  * + \sum_{\ell \in L} \left(
406  * \int_{I_{\ell}} ( \phi - u^{(\ell)} ) dV
407  * \right)^2,
408  * \f]
409  *
410  * subject to the constaint
411  *
412  * \f[
413  * C = \int_{I_0} ( \phi - u^{(0)} ) dV = 0.
414  * \f]
415  *
416  * where \f$\ell\f$ ranges over a subset \f$L\f$ of the secondary neighbors.
417  * \f$L\f$ excludes the one (or more) secondary neighbor(s) where the mean of
418  * \f$u\f$ is the most different from the mean of \f$u^{(0)}\f$. Typically, only
419  * one secondary neighbor is excluded, so \f$L\f$ contains two fewer neighbors
420  * than the total number of neighbors to the element. Note that in 1D, this
421  * implies that \f$L\f$ is the empty set; for each modified solution, one
422  * neighbor is the primary neighbor and the other is the excluded neighbor.
423  *
424  * The integrals are evaluated by quadrature. We denote the quadrature weights
425  * by \f$w_s\f$ and the values of some data \f$X\f$ at the quadrature
426  * nodes by \f$X_s\f$. We use subscripts \f$r,s\f$ to denote quadrature nodes
427  * on the neighbor and local elements, respectively. The minimization becomes
428  *
429  * \f[
430  * \chi^2 = \sum_r w^{(1)}_r ( \phi_r - u^{(1)}_r )^2
431  * + \sum_{\ell \in L} \left(
432  * \sum_r w^{(\ell)}_r ( \phi_r - u^{(\ell)}_r )
433  * \right)^2,
434  * \f]
435  *
436  * subject to the constraint
437  *
438  * \f[
439  * C = \sum_s w^{(0)}_s ( \phi_s - u^{(0)}_s ) = 0.
440  * \f]
441  *
442  * Note that \f$\phi\f$ is a function defined on the local element \f$I_0\f$,
443  * and so is fully represented by its values \f$\phi_s\f$ at the quadrature
444  * points on this element. When evaluating \f$\phi\f$ on element \f$I_{\ell}\f$,
445  * we obtain the function values \f$\phi_r\f$ by polynomial extrapolation,
446  * \f$\phi_r = \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s\f$, where
447  * \f$\mathcal{I}^{(\ell)}_{rs}\f$ is the interpolation/extrapolation matrix
448  * that interpolates data defined at grid points \f$x_s\f$ and evaluates it at
449  * grid points \f$x_r\f$ of \f$I_{\ell}\f$. Thus,
450  *
451  * \f[
452  * \chi^2 = \sum_r w^{(1)}_r
453  * \left(
454  * \sum_s \mathcal{I}^{(1)}_{rs} \phi_s - u^{(1)}_r
455  * \right)^2
456  * + \sum_{\ell \in L} \left(
457  * \sum_r w^{(\ell)}_r
458  * \left(
459  * \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s
460  * - u^{(\ell)}_r
461  * \right)
462  * \right)^2.
463  * \f]
464  *
465  * The solution to this optimization problem is found in the standard way,
466  * using a Lagrange multiplier \f$\lambda\f$ to impose the constraint:
467  *
468  * \f[
469  * 0 = \frac{d}{d \phi_s} \left( \chi^2 + \lambda C \right).
470  * \f]
471  *
472  * Working out the differentiation with respect to \f$\phi_s\f$ leads to the
473  * linear problem that must be inverted to obtain the solution,
474  *
475  * \f[
476  * 0 = A_{st} \phi_t - b_s - \lambda w^{(0)}_s,
477  * \f]
478  *
479  * where
480  *
481  * \f{align*}{
482  * A_{st} &= \sum_r \left( w^{(1)}_r
483  * \mathcal{I}^{(1)}_{rs} \mathcal{I}^{(1)}_{rt}
484  * \right)
485  * + \sum_{\ell \in L} \left(
486  * \sum_r \left(
487  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rt}
488  * \right)
489  * \cdot
490  * \sum_r \left(
491  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
492  * \right)
493  * \right)
494  * \\
495  * b_s &= \sum_r \left( w^{(1)}_r u^{(1)}_r \mathcal{I}^{(1)}_{rs} \right)
496  * + \sum_{\ell \in L} \left(
497  * \sum_r \left( w^{(\ell)}_r u^{(\ell)}_r \right)
498  * \cdot
499  * \sum_r \left(
500  * w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
501  * \right)
502  * \right).
503  * \f}
504  *
505  * Finally, the solution to the constrained fit is
506  *
507  * \f{align*}{
508  * \lambda &= - \frac{ \sum_s w^{(0)}_s
509  * \left( (A^{-1})_{st} b_t - u^{(0)}_s \right)
510  * }{ \sum_s w^{(0)}_s (A^{-1})_{st} w^{(0)}_t }
511  * \\
512  * \phi_s &= (A^{-1})_{st} ( b_t + \lambda w^{(0)}_t ).
513  * \f}
514  *
515  * Note that the matrix \f$A\f$ does not depend on the values of the tensor
516  * \f$u\f$, so its inverse \f$A^{-1}\f$ can be precomputed and stored.
517  *
518  * \warning
519  * Note also that the implementation currently does not support h- or
520  * p-refinement; this is checked by some assertions. The implementation is
521  * untested for grids where elements are curved, and it should not be expected
522  * to work in these cases.
523  */
524 template <typename Tag, size_t VolumeDim, typename Package>
526  const gsl::not_null<db::item_type<Tag>*> modified_tensor,
527  const db::item_type<Tag>& local_tensor, const Element<VolumeDim>& element,
528  const Mesh<VolumeDim>& mesh,
529  const LimiterNeighborData<VolumeDim, Package>& neighbor_data,
530  const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
531  primary_neighbor) noexcept {
532  ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
533  element, primary_neighbor.first),
534  "Found some amount of h-refinement; this is not supported");
535  alg::for_each(neighbor_data, [&mesh](const auto& neighbor_and_data) noexcept {
536  ASSERT(neighbor_and_data.second.mesh == mesh,
537  "Found some amount of p-refinement; this is not supported");
538  });
539 
540  for (size_t tensor_index = 0; tensor_index < local_tensor.size();
541  ++tensor_index) {
542  const auto& tensor_component = local_tensor[tensor_index];
543  const auto neighbors_to_exclude =
544  Hweno_detail::secondary_neighbors_to_exclude_from_fit<Tag>(
545  mean_value(tensor_component, mesh), tensor_index, neighbor_data,
546  primary_neighbor);
547  Hweno_detail::solve_constrained_fit<Tag>(
548  make_not_null(&(*modified_tensor)[tensor_index]),
549  local_tensor[tensor_index], tensor_index, element, mesh, neighbor_data,
550  primary_neighbor, neighbors_to_exclude);
551  }
552 }
553 
554 } // 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:49
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:62
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
void hweno_modified_neighbor_solution(const gsl::not_null< db::item_type< Tag > *> modified_tensor, const db::item_type< Tag > &local_tensor, const Element< VolumeDim > &element, const Mesh< VolumeDim > &mesh, const LimiterNeighborData< VolumeDim, Package > &neighbor_data, const std::pair< Direction< VolumeDim >, ElementId< VolumeDim >> &primary_neighbor) noexcept
Compute the HWENO modified solution for a particular tensor from a particular neighbor element...
Definition: HwenoModifiedSolution.hpp:525
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:290
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 grid-function over a manifold. .
Definition: MeanValue.hpp:38
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:30
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 is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:436
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:863
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 function 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: HwenoModifiedSolution.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:20
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12