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