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>
8 : #include <cstddef>
9 : #include <limits>
10 : #include <string>
11 : #include <type_traits>
12 : #include <unordered_map>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "DataStructures/TaggedTuple.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Domain/SizeOfElement.hpp"
20 : #include "Domain/Structure/Element.hpp"
21 : #include "Domain/Structure/OrientationMapHelpers.hpp"
22 : #include "Domain/Tags.hpp"
23 : #include "Evolution/DiscontinuousGalerkin/Limiters/HwenoImpl.hpp"
24 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodTci.hpp"
25 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
26 : #include "Evolution/DiscontinuousGalerkin/Limiters/SimpleWenoImpl.hpp"
27 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
28 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoType.hpp"
30 : #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
31 : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
32 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
33 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
34 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
35 : #include "Options/String.hpp"
36 : #include "Utilities/Algorithm.hpp"
37 : #include "Utilities/ErrorHandling/Assert.hpp"
38 : #include "Utilities/Gsl.hpp"
39 : #include "Utilities/MakeArray.hpp"
40 : #include "Utilities/TMPL.hpp"
41 :
42 : /// \cond
43 : template <size_t VolumeDim>
44 : class Direction;
45 : template <size_t VolumeDim>
46 : class ElementId;
47 :
48 : namespace PUP {
49 : class er;
50 : } // namespace PUP
51 : /// \endcond
52 :
53 : namespace Limiters {
54 : /// \ingroup LimitersGroup
55 : /// \brief A compact-stencil WENO limiter for DG
56 : ///
57 : /// Implements the simple WENO limiter of \cite Zhong2013 and the Hermite WENO
58 : /// (HWENO) limiter of \cite Zhu2016. The implementation is system-agnostic and
59 : /// can act on an arbitrary set of tensors.
60 : ///
61 : /// #### Summary of the compact-stencil WENO algorithms:
62 : //
63 : /// The compact-stencil WENO limiters require communication only between
64 : /// nearest-neighbor elements, but aim to preserve the full order of the DG
65 : /// solution when the solution is smooth. To achieve this, full volume data is
66 : /// communicated between neighbors.
67 : //
68 : /// For each tensor component to limit, the new solution is obtained by a
69 : /// standard WENO procedure --- the new solution is a linear combination of
70 : /// different polynomials, with weights chosen so that the smoother (i.e., less
71 : /// oscillatory) polynomials contribute the most to the sum.
72 : ///
73 : /// For the simple WENO and HWENO limiters, the polynomials used are the local
74 : /// DG solution as well as a "modified" solution from each neighbor element. For
75 : /// the simple WENO limiter, the modified solution is obtained by simply
76 : /// extrapolating the neighbor solution onto the troubled element. For the HWENO
77 : /// limiter, the modified solution is obtained by a least-squares fit to the
78 : /// solution across multiple neighboring elements.
79 : ///
80 : /// #### Notes on the SpECTRE implemention of the WENO limiters:
81 : ///
82 : /// There are a few differences between the limiters as implemented in SpECTRE
83 : /// and as presented in the references. We list them here and discuss them
84 : /// further below.
85 : /// 1. The choice of basis to represent the DG solution
86 : /// 2. The system-agnostic implementation
87 : /// 3. The oscillation indicator
88 : ///
89 : /// Finally, in 4., we will discuss the geometric limitations of the
90 : /// implementation (which are not a deviation from the references).
91 : ///
92 : /// ##### 1. The choice of basis
93 : ///
94 : /// SpECTRE uses a Legendre basis, rather than the polynomial basis that we
95 : /// understand to be used in the references. Because the construction of the
96 : /// modified neighbor solutions and the WENO sum is geometrically motivated, the
97 : /// overall algorithm should work similarly. However, the precise numerics may
98 : /// differ.
99 : ///
100 : /// ##### 2. The system-agnostic implementation
101 : //
102 : /// This implementation can act on an arbitrary set of tensors. To reach this
103 : /// generality, our HWENO implementation uses a different troubled-cell
104 : /// indicator (TCI) than the reference, which instead specializes the TCI to the
105 : /// Newtonian Euler system of equations.
106 : ///
107 : /// This implementation uses the minmod-based TVB TCI of \cite Cockburn1999 to
108 : /// identify elements that need limiting. The simple WENO implementation follows
109 : /// its reference: it checks the TCI independently for each tensor component, so
110 : /// that only certain tensor components may be limited. The HWENO implementation
111 : /// checks the TVB TCI for all tensor components, and if any single component is
112 : /// troubled, then all components of all tensors are limited.
113 : ///
114 : /// When the evolution system has multiple evolved variables, the recommendation
115 : /// of the references is to apply the limiter to the system's characteristic
116 : /// variables to reduce spurious post-limiting oscillations. In SpECTRE,
117 : /// applying the limiter to the characteristic variables requires specializing
118 : /// the limiter to each evolution system. The system-specific limiter can also
119 : /// implement a system-specific TCI (as the HWENO reference does) to more
120 : /// precisely trigger the limiter.
121 : ///
122 : /// ##### 3. The oscillation indicator
123 : ///
124 : /// We use the oscillation indicator of \cite Dumbser2007, modified for use on
125 : /// the square/cube grids of SpECTRE. We favor this indicator because portions
126 : /// of the work can be precomputed, leading to an oscillation measure that is
127 : /// efficient to evaluate.
128 : ///
129 : /// ##### 4. The geometric limitations
130 : ///
131 : /// Does not support non-Legendre bases; this is checked in DEBUG mode. In
132 : /// principle other bases could be supported, but this would require
133 : /// generalizing the many internal algorithms that assume a Legendre basis.
134 : ///
135 : /// Does not support h- or p-refinement; this is checked always. In principle
136 : /// this could be supported. The modified neighbor solution algorithm would
137 : /// need to be generalized (reasonable for simple WENO, extremely tedious for
138 : /// HWENO), and the sum of neighbor solutions may need to be updated as well.
139 : ///
140 : /// Does not support curved elements; this is not enforced. The code will run
141 : /// but we make no guarantees about the results. Specifically, the limiter acts
142 : /// in the `Frame::ElementLogical` coordinates, because in these coordinates it
143 : /// is straightforward to formulate the algorithm. This means the limiter can
144 : /// operate on generic deformed grids --- however, some things can start to
145 : /// break down, especially on strongly deformed grids:
146 : /// 1. When the Jacobian (from `Frame::ElementLogical` to `Frame::Inertial`)
147 : /// varies across the element, then the limiter fails to be conservative.
148 : /// This is because the integral of a tensor `u` over the element will change
149 : /// after the limiter activates on `u`.
150 : /// 2. When computing the modified neighbor solution for the WENO sum, the
151 : /// extrapolation or fitting procedure may not properly account for the
152 : /// coordinates of the source data. If the coordinate map of the neighbor
153 : /// differs from that of the local element, then the logical-coordinate
154 : /// representation of the neighbor data may be incorrect. This may be a
155 : /// large error at Block boundaries with discontinuous map changes, and may
156 : /// be a small error from smoothly-varying maps that are not sufficiently
157 : /// resolved from one element to the next.
158 : template <size_t VolumeDim, typename TagsToLimit>
159 1 : class Weno;
160 :
161 : template <size_t VolumeDim, typename... Tags>
162 0 : class Weno<VolumeDim, tmpl::list<Tags...>> {
163 : public:
164 : /// \brief The WenoType
165 : ///
166 : /// One of `Limiters::WenoType`. See the `Limiters::Weno`
167 : /// documentation for details.
168 1 : struct Type {
169 0 : using type = WenoType;
170 0 : static constexpr Options::String help = {"Type of WENO limiter"};
171 : };
172 : /// \brief The linear weight given to each neighbor
173 : ///
174 : /// This linear weight gets combined with the oscillation indicator to
175 : /// compute the weight for each WENO estimated solution. The standard value
176 : /// in the literature is 0.001; larger values may be better suited for
177 : /// problems with strong shocks, and smaller values may be better suited to
178 : /// smooth problems.
179 1 : struct NeighborWeight {
180 0 : using type = double;
181 0 : static type lower_bound() { return 1e-6; }
182 0 : static type upper_bound() { return 0.1; }
183 0 : static constexpr Options::String help = {
184 : "Linear weight for each neighbor element's solution"};
185 : };
186 : /// \brief The TVB constant for the minmod TCI
187 : ///
188 : /// See `Limiters::Minmod` documentation for details.
189 1 : struct TvbConstant {
190 0 : using type = double;
191 0 : static type lower_bound() { return 0.0; }
192 0 : static constexpr Options::String help = {"TVB constant 'm'"};
193 : };
194 : /// \brief Turn the limiter off
195 : ///
196 : /// This option exists to temporarily disable the limiter for debugging
197 : /// purposes. For problems where limiting is not needed, the preferred
198 : /// approach is to not compile the limiter into the executable.
199 1 : struct DisableForDebugging {
200 0 : using type = bool;
201 0 : static type suggested_value() { return false; }
202 0 : static constexpr Options::String help = {"Disable the limiter"};
203 : };
204 0 : using options =
205 : tmpl::list<Type, NeighborWeight, TvbConstant, DisableForDebugging>;
206 0 : static constexpr Options::String help = {"A WENO limiter for DG"};
207 :
208 0 : Weno(WenoType weno_type, double neighbor_linear_weight, double tvb_constant,
209 : bool disable_for_debugging = false);
210 :
211 0 : Weno() = default;
212 0 : Weno(const Weno& /*rhs*/) = default;
213 0 : Weno& operator=(const Weno& /*rhs*/) = default;
214 0 : Weno(Weno&& /*rhs*/) = default;
215 0 : Weno& operator=(Weno&& /*rhs*/) = default;
216 0 : ~Weno() = default;
217 :
218 : // NOLINTNEXTLINE(google-runtime-references)
219 0 : void pup(PUP::er& p);
220 :
221 : /// \brief Data to send to neighbor elements
222 1 : struct PackagedData {
223 0 : Variables<tmpl::list<Tags...>> volume_data;
224 0 : tuples::TaggedTuple<::Tags::Mean<Tags>...> means;
225 0 : Mesh<VolumeDim> mesh;
226 0 : std::array<double, VolumeDim> element_size =
227 : make_array<VolumeDim>(std::numeric_limits<double>::signaling_NaN());
228 :
229 : // NOLINTNEXTLINE(google-runtime-references)
230 0 : void pup(PUP::er& p) {
231 : p | volume_data;
232 : p | means;
233 : p | mesh;
234 : p | element_size;
235 : }
236 : };
237 :
238 0 : using package_argument_tags =
239 : tmpl::list<Tags..., domain::Tags::Mesh<VolumeDim>,
240 : domain::Tags::SizeOfElement<VolumeDim>>;
241 :
242 : /// \brief Package data for sending to neighbor elements
243 1 : void package_data(gsl::not_null<PackagedData*> packaged_data,
244 : const typename Tags::type&... tensors,
245 : const Mesh<VolumeDim>& mesh,
246 : const std::array<double, VolumeDim>& element_size,
247 : const OrientationMap<VolumeDim>& orientation_map) const;
248 :
249 0 : using limit_tags = tmpl::list<Tags...>;
250 0 : using limit_argument_tags =
251 : tmpl::list<domain::Tags::Mesh<VolumeDim>,
252 : domain::Tags::Element<VolumeDim>,
253 : domain::Tags::SizeOfElement<VolumeDim>>;
254 :
255 : /// \brief Limit the solution on the element
256 1 : bool operator()(
257 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
258 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
259 : const std::array<double, VolumeDim>& element_size,
260 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
261 : boost::hash<DirectionalId<VolumeDim>>>&
262 : neighbor_data) const;
263 :
264 : private:
265 : template <size_t LocalDim, typename LocalTagList>
266 : // NOLINTNEXTLINE(readability-redundant-declaration) false positive
267 0 : friend bool operator==(const Weno<LocalDim, LocalTagList>& lhs,
268 : const Weno<LocalDim, LocalTagList>& rhs);
269 :
270 0 : WenoType weno_type_;
271 0 : double neighbor_linear_weight_;
272 0 : double tvb_constant_;
273 0 : bool disable_for_debugging_;
274 : };
275 :
276 : template <size_t VolumeDim, typename... Tags>
277 : Weno<VolumeDim, tmpl::list<Tags...>>::Weno(const WenoType weno_type,
278 : const double neighbor_linear_weight,
279 : const double tvb_constant,
280 : const bool disable_for_debugging)
281 : : weno_type_(weno_type),
282 : neighbor_linear_weight_(neighbor_linear_weight),
283 : tvb_constant_(tvb_constant),
284 : disable_for_debugging_(disable_for_debugging) {}
285 :
286 : template <size_t VolumeDim, typename... Tags>
287 : // NOLINTNEXTLINE(google-runtime-references)
288 : void Weno<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) {
289 : p | weno_type_;
290 : p | neighbor_linear_weight_;
291 : p | tvb_constant_;
292 : p | disable_for_debugging_;
293 : }
294 :
295 : template <size_t VolumeDim, typename... Tags>
296 : void Weno<VolumeDim, tmpl::list<Tags...>>::package_data(
297 : const gsl::not_null<PackagedData*> packaged_data,
298 : const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
299 : const std::array<double, VolumeDim>& element_size,
300 : const OrientationMap<VolumeDim>& orientation_map) const {
301 : // By always initializing the PackagedData Variables member, we avoid an
302 : // assertion that arises from having a default-constructed Variables in a
303 : // disabled limiter. There is a performance cost, because the package_data()
304 : // function does non-zero work even for a disabled limiter... but since the
305 : // limiter should never be disabled in a production simulation, this cost
306 : // should never matter.
307 : (packaged_data->volume_data).initialize(mesh.number_of_grid_points());
308 :
309 : if (UNLIKELY(disable_for_debugging_)) {
310 : // Do not initialize packaged_data
311 : // (except for the Variables member "volume_data", see above)
312 : return;
313 : }
314 :
315 : const auto wrap_compute_means = [&mesh, &packaged_data](auto tag,
316 : const auto tensor) {
317 : for (size_t i = 0; i < tensor.size(); ++i) {
318 : // Compute the mean using the local orientation of the tensor and mesh.
319 : get<::Tags::Mean<decltype(tag)>>(packaged_data->means)[i] =
320 : mean_value(tensor[i], mesh);
321 : }
322 : return '0';
323 : };
324 : expand_pack(wrap_compute_means(Tags{}, tensors)...);
325 :
326 : packaged_data->element_size =
327 : orientation_map.permute_from_neighbor(element_size);
328 :
329 : const auto wrap_copy_tensor = [&packaged_data](auto tag, const auto tensor) {
330 : get<decltype(tag)>(packaged_data->volume_data) = tensor;
331 : return '0';
332 : };
333 : expand_pack(wrap_copy_tensor(Tags{}, tensors)...);
334 : packaged_data->volume_data = orient_variables(
335 : packaged_data->volume_data, mesh.extents(), orientation_map);
336 :
337 : // Warning: the WENO limiter is currently only tested with aligned meshes.
338 : // The orientation of the mesh, the `element_size` computed above, and the
339 : // variables should be carefully tested when used with domains that involve
340 : // orientation maps
341 : packaged_data->mesh = orientation_map(mesh);
342 : }
343 :
344 : template <size_t VolumeDim, typename... Tags>
345 : bool Weno<VolumeDim, tmpl::list<Tags...>>::operator()(
346 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
347 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
348 : const std::array<double, VolumeDim>& element_size,
349 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
350 : boost::hash<DirectionalId<VolumeDim>>>&
351 : neighbor_data) const {
352 : if (UNLIKELY(disable_for_debugging_)) {
353 : // Do not modify input tensors
354 : return false;
355 : }
356 :
357 : // Check that basis is LGL or LG
358 : // A Legendre basis is assumed for the oscillation indicator (used in both
359 : // SimpleWeno and Hweno) and in the Hweno reconstruction.
360 : ASSERT(mesh.basis() == make_array<VolumeDim>(Spectral::Basis::Legendre),
361 : "Unsupported basis: " << mesh);
362 : ASSERT(mesh.quadrature() ==
363 : make_array<VolumeDim>(Spectral::Quadrature::GaussLobatto) or
364 : mesh.quadrature() ==
365 : make_array<VolumeDim>(Spectral::Quadrature::Gauss),
366 : "Unsupported quadrature: " << mesh);
367 :
368 : // Enforce restrictions on h-refinement, p-refinement
369 : if (UNLIKELY(
370 : alg::any_of(element.neighbors(), [](const auto& direction_neighbors) {
371 : return direction_neighbors.second.size() != 1;
372 : }))) {
373 : ERROR("The Weno limiter does not yet support h-refinement");
374 : // Removing this limitation will require:
375 : // - Generalizing the computation of the modified neighbor solutions.
376 : // - Generalizing the WENO weighted sum for multiple neighbors in each
377 : // direction.
378 : }
379 : alg::for_each(neighbor_data, [&mesh](const auto& neighbor_and_data) {
380 : if (UNLIKELY(neighbor_and_data.second.mesh != mesh)) {
381 : ERROR("The Weno limiter does not yet support p-refinement");
382 : // Removing this limitation will require generalizing the
383 : // computation of the modified neighbor solutions.
384 : }
385 : });
386 :
387 : if (weno_type_ == WenoType::Hweno) {
388 : // Troubled-cell detection for HWENO flags the element for limiting if any
389 : // component of any tensor needs limiting.
390 : const bool cell_is_troubled =
391 : Tci::tvb_minmod_indicator<VolumeDim, PackagedData, Tags...>(
392 : tvb_constant_, (*tensors)..., mesh, element, element_size,
393 : neighbor_data);
394 : if (not cell_is_troubled) {
395 : // No limiting is needed
396 : return false;
397 : }
398 :
399 : std::unordered_map<DirectionalId<VolumeDim>, DataVector,
400 : boost::hash<DirectionalId<VolumeDim>>>
401 : modified_neighbor_solution_buffer{};
402 : for (const auto& neighbor_and_data : neighbor_data) {
403 : const auto& neighbor = neighbor_and_data.first;
404 : modified_neighbor_solution_buffer.insert(
405 : std::make_pair(neighbor, DataVector(mesh.number_of_grid_points())));
406 : }
407 :
408 : EXPAND_PACK_LEFT_TO_RIGHT(Weno_detail::hweno_impl<Tags>(
409 : make_not_null(&modified_neighbor_solution_buffer), tensors,
410 : neighbor_linear_weight_, mesh, element, neighbor_data));
411 : return true; // cell_is_troubled
412 :
413 : } else if (weno_type_ == WenoType::SimpleWeno) {
414 : // Buffers and pre-computations for TCI
415 : Minmod_detail::BufferWrapper<VolumeDim> tci_buffer(mesh);
416 : const auto effective_neighbor_sizes =
417 : Minmod_detail::compute_effective_neighbor_sizes(element, neighbor_data);
418 :
419 : // Buffers for simple WENO implementation
420 : std::unordered_map<DirectionalId<VolumeDim>, intrp::RegularGrid<VolumeDim>,
421 : boost::hash<DirectionalId<VolumeDim>>>
422 : interpolator_buffer{};
423 : std::unordered_map<DirectionalId<VolumeDim>, DataVector,
424 : boost::hash<DirectionalId<VolumeDim>>>
425 : modified_neighbor_solution_buffer{};
426 :
427 : bool some_component_was_limited = false;
428 :
429 : const auto wrap_minmod_tci_and_simple_weno_impl =
430 : [this, &some_component_was_limited, &tci_buffer, &interpolator_buffer,
431 : &modified_neighbor_solution_buffer, &mesh, &element, &element_size,
432 : &neighbor_data,
433 : &effective_neighbor_sizes](auto tag, const auto tensor) {
434 : for (size_t tensor_storage_index = 0;
435 : tensor_storage_index < tensor->size(); ++tensor_storage_index) {
436 : // Check TCI
437 : const auto effective_neighbor_means =
438 : Minmod_detail::compute_effective_neighbor_means<decltype(tag)>(
439 : tensor_storage_index, element, neighbor_data);
440 : const bool component_needs_limiting = Tci::tvb_minmod_indicator(
441 : make_not_null(&tci_buffer), tvb_constant_,
442 : (*tensor)[tensor_storage_index], mesh, element, element_size,
443 : effective_neighbor_means, effective_neighbor_sizes);
444 :
445 : if (component_needs_limiting) {
446 : if (modified_neighbor_solution_buffer.empty()) {
447 : // Allocate the neighbor solution buffers only if the limiter is
448 : // triggered. This reduces allocation when no limiting occurs.
449 : for (const auto& neighbor_and_data : neighbor_data) {
450 : const auto& neighbor = neighbor_and_data.first;
451 : modified_neighbor_solution_buffer.insert(std::make_pair(
452 : neighbor, DataVector(mesh.number_of_grid_points())));
453 : }
454 : }
455 : Weno_detail::simple_weno_impl<decltype(tag)>(
456 : make_not_null(&interpolator_buffer),
457 : make_not_null(&modified_neighbor_solution_buffer), tensor,
458 : neighbor_linear_weight_, tensor_storage_index, mesh, element,
459 : neighbor_data);
460 : some_component_was_limited = true;
461 : }
462 : }
463 : return '0';
464 : };
465 : expand_pack(wrap_minmod_tci_and_simple_weno_impl(Tags{}, tensors)...);
466 : return some_component_was_limited; // cell_is_troubled
467 : } else {
468 : ERROR("WENO limiter not implemented for WenoType: " << weno_type_);
469 : }
470 :
471 : return false; // cell_is_troubled
472 : }
473 :
474 : template <size_t LocalDim, typename LocalTagList>
475 0 : bool operator==(const Weno<LocalDim, LocalTagList>& lhs,
476 : const Weno<LocalDim, LocalTagList>& rhs) {
477 : return lhs.weno_type_ == rhs.weno_type_ and
478 : lhs.neighbor_linear_weight_ == rhs.neighbor_linear_weight_ and
479 : lhs.tvb_constant_ == rhs.tvb_constant_ and
480 : lhs.disable_for_debugging_ == rhs.disable_for_debugging_;
481 : }
482 :
483 : template <size_t VolumeDim, typename TagList>
484 0 : bool operator!=(const Weno<VolumeDim, TagList>& lhs,
485 : const Weno<VolumeDim, TagList>& rhs) {
486 : return not(lhs == rhs);
487 : }
488 :
489 : } // namespace Limiters
|
