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 <limits>
8 : #include <type_traits>
9 : #include <unordered_map>
10 : #include <utility>
11 :
12 : #include "DataStructures/Tags.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodImpl.hpp"
15 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
16 : #include "Options/String.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/MakeArray.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TaggedTuple.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : template <size_t VolumeDim>
25 : class Direction;
26 : template <size_t Dim, typename T>
27 : class DirectionMap;
28 : template <size_t VolumeDim>
29 : class Element;
30 : template <size_t VolumeDim>
31 : class ElementId;
32 : template <size_t VolumeDim>
33 : class Mesh;
34 : template <size_t VolumeDim>
35 : class OrientationMap;
36 :
37 : namespace boost {
38 : template <class T>
39 : struct hash;
40 : } // namespace boost
41 :
42 : namespace PUP {
43 : class er;
44 : } // namespace PUP
45 :
46 : namespace Limiters::Minmod_detail {
47 : template <size_t VolumeDim>
48 : class BufferWrapper;
49 : } // namespace Limiters::Minmod_detail
50 :
51 : namespace domain::Tags {
52 : template <size_t Dim, typename Frame>
53 : struct Coordinates;
54 : template <size_t VolumeDim>
55 : struct Element;
56 : template <size_t VolumeDim>
57 : struct Mesh;
58 : template <size_t VolumeDim>
59 : struct SizeOfElement;
60 : } // namespace domain::Tags
61 : /// \endcond
62 :
63 : namespace Limiters {
64 : /// \ingroup LimitersGroup
65 : /// \brief A minmod-based generalized slope limiter
66 : ///
67 : /// Implements the three minmod-based generalized slope limiters from
68 : /// \cite Cockburn1999 Sec. 2.4: \f$\Lambda\Pi^1\f$, \f$\Lambda\Pi^N\f$, and
69 : /// MUSCL. The implementation is system-agnostic and can act on an arbitrary
70 : /// set of tensors.
71 : ///
72 : /// #### Summary of the generalized slope limiter algorithms:
73 : ///
74 : /// The MUSCL and \f$\Lambda\Pi^1\f$ limiters are both intended for use on
75 : /// piecewise-linear solutions, i.e., on linear-order elements with two points
76 : /// per dimension. These limiters operate by reducing the spatial slope of the
77 : /// tensor components if the data look like they may contain oscillations.
78 : /// Between these two, MUSCL is more dissipative --- it more aggressively
79 : /// reduces the slopes of the data, so it may better handle strong shocks, but
80 : /// it correspondingly produces the most broadening of features in the solution.
81 : ///
82 : /// Note that we do not _require_ the MUSCL and \f$\Lambda\Pi^1\f$ limiters to
83 : /// be used on linear-order elements. However, when they are used on a
84 : /// higher-resolution grid, the limiters act to linearize the solution (by
85 : /// discarding higher-order mode content) whether or not the slopes must be
86 : /// reduced.
87 : ///
88 : /// The \f$\Lambda\Pi^N\f$ limiter is intended for use with higher-order
89 : /// elements (with more than two points per dimension), where the solution is a
90 : /// piecewise polynomial of higher-than-linear order. This limiter generalizes
91 : /// \f$\Lambda\Pi^1\f$: the post-limiter solution is the linearized solution of
92 : /// \f$\Lambda\Pi^1\f$ in the case that the slopes must be reduced, but is the
93 : /// original (higher-order) data in the case that the slopes are acceptable.
94 : ///
95 : /// For all three types of minmod limiter, the algorithm can be relaxed from
96 : /// TVD (total variation diminishing) in the means to TVB (total variation
97 : /// bound) in the means. This may avoid limiting away smooth extrema in the
98 : /// solution that would otherwise look like spurious oscillations. When this
99 : /// correction is enabled, the limiter will not reduce the slope (but may still
100 : /// linearize) on elements where the slope is less than \f$m h^2\f$, where
101 : /// \f$m\f$ is the TVB constant and \f$h\f$ is the size of the DG element.
102 : /// Note the "in the means" qualifier: the limiter controls the oscillation
103 : /// between the mean solution values across neighboring cells, but may not
104 : /// control oscillations within the cells.
105 : ///
106 : /// The choice of the TVB constant \f$m\f$ is difficult. Larger values result in
107 : /// fewer limiter activations, especially near smooth extrema in the solution
108 : /// --- this can help to avoid incorrectly limiting away these smooth extrema,
109 : /// but can also result in insufficient limiting of truly spurious oscillations.
110 : /// The reference uses a value of 50 when presenting the limiter with simple
111 : /// shock tests, but in general the value \f$m\f$ that optimizes between
112 : /// robustness and minimal loss of accuracy is problem dependent.
113 : ///
114 : /// #### Notes on the SpECTRE implementation of the generalized slope limiters:
115 : ///
116 : /// This implementation can act on an arbitrary set of tensors; the limiting
117 : /// algorithm is applied to each component of each tensor independently. This
118 : /// is a convenient and general interface. However, when the evolution system
119 : /// has multiple evolved variables, the recommendation of the reference is to
120 : /// apply the limiter to the system's characteristic variables to reduce
121 : /// spurious post-limiting oscillations. In SpECTRE, applying the limiter to
122 : /// the characteristic variables requires specializing the limiter to each
123 : /// evolution system.
124 : ///
125 : /// The limiter acts in the `Frame::ElementLogical` coordinates, because in
126 : /// these coordinates it is straightforward to formulate the algorithm. This
127 : /// means the limiter can operate on generic deformed grids --- however, some
128 : /// things can start to break down, especially on strongly deformed grids:
129 : /// 1. When the Jacobian (from `Frame::ElementLogical` to `Frame::Inertial`)
130 : /// varies across the element, then the limiter fails to be conservative.
131 : /// This is because the integral of a tensor `u` over the element will change
132 : /// after the limiter activates on `u`.
133 : /// 2. When there is a sudden change in the size of the elements (perhaps at an
134 : /// h-refinement boundary, or at the boundary between two blocks with very
135 : /// different mappings), a smooth solution in `Frame::Inertial` can appear
136 : /// to have a kink in `Frame::ElementLogical`. The Minmod implementation
137 : /// includes some (tested but unproven) corrections based on the size of the
138 : /// elements that try to reduce spurious limiter activations near these fake
139 : /// kinks.
140 : ///
141 : /// When an element has multiple neighbors in any direction, an effective mean
142 : /// and neighbor size in this direction are computed by averaging over the
143 : /// multiple neighbors. This simple generalization of the minmod limiter enables
144 : /// it to operate on h-refined grids.
145 : ///
146 : /// \tparam VolumeDim The number of spatial dimensions.
147 : /// \tparam TagsToLimit A typelist of tags specifying the tensors to limit.
148 : template <size_t VolumeDim, typename TagsToLimit>
149 1 : class Minmod;
150 :
151 : template <size_t VolumeDim, typename... Tags>
152 0 : class Minmod<VolumeDim, tmpl::list<Tags...>> {
153 : public:
154 : /// \brief The MinmodType
155 : ///
156 : /// One of `Limiters::MinmodType`. See `Limiters::Minmod`
157 : /// documentation for details. Note in particular that on grids with more than
158 : /// two points per dimension, the recommended type is `LambdaPiN`.
159 1 : struct Type {
160 0 : using type = MinmodType;
161 0 : static constexpr Options::String help = {"Type of minmod"};
162 : };
163 : /// \brief The TVB constant
164 : ///
165 : /// See `Limiters::Minmod` documentation for details. The optimal value of
166 : /// this parameter is unfortunately problem-dependent.
167 1 : struct TvbConstant {
168 0 : using type = double;
169 0 : static type lower_bound() { return 0.0; }
170 0 : static constexpr Options::String help = {"TVB constant 'm'"};
171 : };
172 : /// \brief Turn the limiter off
173 : ///
174 : /// This option exists to temporarily disable the limiter for debugging
175 : /// purposes. For problems where limiting is not needed, the preferred
176 : /// approach is to not compile the limiter into the executable.
177 1 : struct DisableForDebugging {
178 0 : using type = bool;
179 0 : static type suggested_value() { return false; }
180 0 : static constexpr Options::String help = {"Disable the limiter"};
181 : };
182 0 : using options = tmpl::list<Type, TvbConstant, DisableForDebugging>;
183 0 : static constexpr Options::String help = {
184 : "A minmod-based generalized slope limiter.\n"
185 : "The different types of minmod are more or less aggressive in trying\n"
186 : "to reduce slopes. The TVB correction allows the limiter to ignore\n"
187 : "'small' slopes, and helps to avoid limiting of smooth extrema in the\n"
188 : "solution.\n"};
189 :
190 : /// \brief Constuct a Minmod slope limiter
191 : ///
192 : /// \param minmod_type The type of Minmod slope limiter.
193 : /// \param tvb_constant The value of the TVB constant.
194 : /// \param disable_for_debugging Switch to turn the limiter off.
195 1 : explicit Minmod(MinmodType minmod_type, double tvb_constant,
196 : bool disable_for_debugging = false);
197 :
198 0 : Minmod() = default;
199 0 : Minmod(const Minmod& /*rhs*/) = default;
200 0 : Minmod& operator=(const Minmod& /*rhs*/) = default;
201 0 : Minmod(Minmod&& /*rhs*/) = default;
202 0 : Minmod& operator=(Minmod&& /*rhs*/) = default;
203 0 : ~Minmod() = default;
204 :
205 : // NOLINTNEXTLINE(google-runtime-references)
206 0 : void pup(PUP::er& p);
207 :
208 : // To facilitate testing
209 : /// \cond
210 : MinmodType minmod_type() const { return minmod_type_; }
211 : /// \endcond
212 :
213 : /// \brief Data to send to neighbor elements.
214 1 : struct PackagedData {
215 0 : tuples::TaggedTuple<::Tags::Mean<Tags>...> means;
216 0 : std::array<double, VolumeDim> element_size =
217 : make_array<VolumeDim>(std::numeric_limits<double>::signaling_NaN());
218 :
219 : // NOLINTNEXTLINE(google-runtime-references)
220 0 : void pup(PUP::er& p) {
221 : p | means;
222 : p | element_size;
223 : }
224 : };
225 :
226 0 : using package_argument_tags =
227 : tmpl::list<Tags..., domain::Tags::Mesh<VolumeDim>,
228 : domain::Tags::SizeOfElement<VolumeDim>>;
229 :
230 : /// \brief Package data for sending to neighbor elements.
231 : ///
232 : /// The following quantities are stored in `PackagedData` and communicated
233 : /// between neighboring elements:
234 : /// - the cell-averaged mean of each tensor component, and
235 : /// - the size of the cell along each logical coordinate direction.
236 : ///
237 : /// \param packaged_data The data package to fill with this element's values.
238 : /// \param tensors The tensors to be averaged and packaged.
239 : /// \param mesh The mesh on which the tensor values are measured.
240 : /// \param element_size The size of the element in inertial coordinates, along
241 : /// each dimension of logical coordinates.
242 : /// \param orientation_map The orientation of the neighbor
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::Coordinates<VolumeDim, Frame::ElementLogical>,
254 : domain::Tags::SizeOfElement<VolumeDim>>;
255 :
256 : /// \brief Limits the solution on the element.
257 : ///
258 : /// For each component of each tensor, the limiter will (in general) linearize
259 : /// the data, then possibly reduce its slope, dimension-by-dimension, until it
260 : /// no longer looks oscillatory.
261 : ///
262 : /// \param tensors The tensors to be limited.
263 : /// \param mesh The mesh on which the tensor values are measured.
264 : /// \param element The element on which the tensors to limit live.
265 : /// \param logical_coords The element logical coordinates of the mesh
266 : /// gridpoints.
267 : /// \param element_size The size of the element, in the inertial coordinates.
268 : /// \param neighbor_data The data from each neighbor.
269 : ///
270 : /// \return whether the limiter modified the solution or not.
271 : ///
272 : /// \note The return value is false if the limiter knows it has not modified
273 : /// the solution. True return values can indicate:
274 : /// - The solution was limited to reduce the slope, whether by a large factor
275 : /// or by a factor only roundoff away from unity.
276 : /// - The solution was linearized but not limited.
277 : /// - The solution is identical to the input, if the input was a linear
278 : /// function on a higher-order mesh, so that the limiter cannot know that
279 : /// the linearization step did not actually modify the data. This is
280 : /// somewhat contrived and is unlikely to occur outside of code tests or
281 : /// test cases with very clean initial data.
282 1 : bool operator()(
283 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
284 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
285 : const tnsr::I<DataVector, VolumeDim, Frame::ElementLogical>&
286 : logical_coords,
287 : const std::array<double, VolumeDim>& element_size,
288 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
289 : boost::hash<DirectionalId<VolumeDim>>>&
290 : neighbor_data) const;
291 :
292 : private:
293 : template <size_t LocalDim, typename LocalTagList>
294 : // NOLINTNEXTLINE(readability-redundant-declaration) false positive
295 0 : friend bool operator==(const Minmod<LocalDim, LocalTagList>& lhs,
296 : const Minmod<LocalDim, LocalTagList>& rhs);
297 :
298 0 : MinmodType minmod_type_;
299 0 : double tvb_constant_;
300 0 : bool disable_for_debugging_;
301 : };
302 :
303 : template <size_t VolumeDim, typename... Tags>
304 : Minmod<VolumeDim, tmpl::list<Tags...>>::Minmod(const MinmodType minmod_type,
305 : const double tvb_constant,
306 : const bool disable_for_debugging)
307 : : minmod_type_(minmod_type),
308 : tvb_constant_(tvb_constant),
309 : disable_for_debugging_(disable_for_debugging) {
310 : ASSERT(tvb_constant >= 0.0, "The TVB constant must be non-negative.");
311 : }
312 :
313 : template <size_t VolumeDim, typename... Tags>
314 : void Minmod<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) {
315 : p | minmod_type_;
316 : p | tvb_constant_;
317 : p | disable_for_debugging_;
318 : }
319 :
320 : template <size_t VolumeDim, typename... Tags>
321 : void Minmod<VolumeDim, tmpl::list<Tags...>>::package_data(
322 : const gsl::not_null<PackagedData*> packaged_data,
323 : const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
324 : const std::array<double, VolumeDim>& element_size,
325 : const OrientationMap<VolumeDim>& orientation_map) const {
326 : if (UNLIKELY(disable_for_debugging_)) {
327 : // Do not initialize packaged_data
328 : return;
329 : }
330 :
331 : const auto wrap_compute_means = [&mesh, &packaged_data](auto tag,
332 : const auto tensor) {
333 : for (size_t i = 0; i < tensor.size(); ++i) {
334 : // Compute the mean using the local orientation of the tensor and mesh:
335 : // this avoids the work of reorienting the tensor while giving the same
336 : // result.
337 : get<::Tags::Mean<decltype(tag)>>(packaged_data->means)[i] =
338 : mean_value(tensor[i], mesh);
339 : }
340 : return '0';
341 : };
342 : expand_pack(wrap_compute_means(Tags{}, tensors)...);
343 : packaged_data->element_size =
344 : orientation_map.permute_from_neighbor(element_size);
345 : }
346 :
347 : template <size_t VolumeDim, typename... Tags>
348 : bool Minmod<VolumeDim, tmpl::list<Tags...>>::operator()(
349 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
350 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
351 : const tnsr::I<DataVector, VolumeDim, Frame::ElementLogical>& logical_coords,
352 : const std::array<double, VolumeDim>& element_size,
353 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
354 : boost::hash<DirectionalId<VolumeDim>>>&
355 : neighbor_data) const {
356 : if (UNLIKELY(disable_for_debugging_)) {
357 : // Do not modify input tensors
358 : return false;
359 : }
360 :
361 : DataVector u_lin_buffer(mesh.number_of_grid_points());
362 : Minmod_detail::BufferWrapper<VolumeDim> buffer(mesh);
363 :
364 : bool limiter_activated = false;
365 : const auto wrap_minmod_impl = [this, &limiter_activated, &element, &mesh,
366 : &logical_coords, &element_size, &neighbor_data,
367 : &u_lin_buffer,
368 : &buffer](auto tag, const auto tensor) {
369 : limiter_activated =
370 : Minmod_detail::minmod_impl<VolumeDim, decltype(tag)>(
371 : &u_lin_buffer, &buffer, tensor, minmod_type_, tvb_constant_, mesh,
372 : element, logical_coords, element_size, neighbor_data) or
373 : limiter_activated;
374 : return '0';
375 : };
376 : expand_pack(wrap_minmod_impl(Tags{}, tensors)...);
377 : return limiter_activated;
378 : }
379 :
380 : template <size_t LocalDim, typename LocalTagList>
381 0 : bool operator==(const Minmod<LocalDim, LocalTagList>& lhs,
382 : const Minmod<LocalDim, LocalTagList>& rhs) {
383 : return lhs.minmod_type_ == rhs.minmod_type_ and
384 : lhs.tvb_constant_ == rhs.tvb_constant_ and
385 : lhs.disable_for_debugging_ == rhs.disable_for_debugging_;
386 : }
387 :
388 : template <size_t VolumeDim, typename TagList>
389 0 : bool operator!=(const Minmod<VolumeDim, TagList>& lhs,
390 : const Minmod<VolumeDim, TagList>& rhs) {
391 : return not(lhs == rhs);
392 : }
393 :
394 : } // namespace Limiters
|
