MinmodImpl.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>
8 #include <cstdlib>
9 #include <iterator>
10 #include <limits>
11 #include <memory>
12 #include <pup.h>
13 #include <type_traits>
14 #include <unordered_map>
15 #include <utility>
16 
17 #include "DataStructures/DataVector.hpp"
18 #include "DataStructures/Tags.hpp"
21 #include "Domain/Structure/DirectionMap.hpp"
24 #include "Domain/Structure/OrientationMap.hpp"
25 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodHelpers.hpp"
26 #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
30 #include "Utilities/Gsl.hpp"
31 #include "Utilities/Literals.hpp"
32 #include "Utilities/Numeric.hpp"
33 #include "Utilities/TMPL.hpp"
34 
35 namespace Limiters::Minmod_detail {
36 
37 // This function combines the evaluation of the troubled-cell indicator with the
38 // computation of the post-limiter reduced slopes. The returned bool indicates
39 // whether the slopes are to be reduced. The slopes themselves are returned by
40 // pointer.
41 //
42 // Note: This function is only made available in this header file to facilitate
43 // testing.
44 template <size_t VolumeDim>
45 bool minmod_limited_slopes(
46  gsl::not_null<DataVector*> u_lin_buffer,
47  gsl::not_null<BufferWrapper<VolumeDim>*> buffer,
50  Limiters::MinmodType minmod_type, double tvb_constant, const DataVector& u,
51  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
52  const std::array<double, VolumeDim>& element_size,
53  const DirectionMap<VolumeDim, double>& effective_neighbor_means,
54  const DirectionMap<VolumeDim, double>& effective_neighbor_sizes) noexcept;
55 
56 // Implements the minmod limiter for one Tensor<DataVector> at a time.
57 template <size_t VolumeDim, typename Tag, typename PackagedData>
58 bool minmod_impl(
59  const gsl::not_null<DataVector*> u_lin_buffer,
60  const gsl::not_null<BufferWrapper<VolumeDim>*> buffer,
62  const Limiters::MinmodType minmod_type, const double tvb_constant,
63  const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
64  const tnsr::I<DataVector, VolumeDim, Frame::Logical>& logical_coords,
65  const std::array<double, VolumeDim>& element_size,
66  const std::unordered_map<
69  neighbor_data) noexcept {
70  // True if the mesh is linear-order in every direction
71  const bool mesh_is_linear = (mesh.extents() == Index<VolumeDim>(2));
72  const bool minmod_type_is_linear =
73  (minmod_type != Limiters::MinmodType::LambdaPiN);
74  const bool using_linear_limiter_on_non_linear_mesh =
75  minmod_type_is_linear and not mesh_is_linear;
76 
77  // In each direction, average the size of all different neighbors in that
78  // direction. Note that only the component of neighor_size that is normal
79  // to the face is needed (and, therefore, computed). Note that this average
80  // does not depend on the solution on the neighboring elements, so could be
81  // precomputed outside of `limit_one_tensor`. Changing the code to
82  // precompute the average may or may not be a measurable optimization.
83  const auto effective_neighbor_sizes =
84  compute_effective_neighbor_sizes(element, neighbor_data);
85 
86  bool some_component_was_limited = false;
87  for (size_t i = 0; i < tensor->size(); ++i) {
88  // In each direction, average the mean of the i'th tensor component over
89  // all different neighbors in that direction. This produces one effective
90  // neighbor per direction.
91  const auto effective_neighbor_means =
92  compute_effective_neighbor_means<Tag>(i, element, neighbor_data);
93 
94  DataVector& u = (*tensor)[i];
96  std::array<double, VolumeDim> u_limited_slopes{};
97  const bool reduce_slopes = minmod_limited_slopes(
98  u_lin_buffer, buffer, make_not_null(&u_mean),
99  make_not_null(&u_limited_slopes), minmod_type, tvb_constant, u, mesh,
100  element, element_size, effective_neighbor_means,
101  effective_neighbor_sizes);
102 
103  if (reduce_slopes or using_linear_limiter_on_non_linear_mesh) {
104  u = u_mean;
105  for (size_t d = 0; d < VolumeDim; ++d) {
106  u += logical_coords.get(d) * gsl::at(u_limited_slopes, d);
107  }
108  some_component_was_limited = true;
109  }
110  } // end for loop over tensor components
111 
112  return some_component_was_limited;
113 }
114 
115 } // namespace Limiters::Minmod_detail
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
utility
Literals.hpp
std::pair
iterator
Index< VolumeDim >
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:50
ElementId.hpp
Assert.hpp
array
Limiters::MinmodType
MinmodType
Possible types of the minmod slope limiter and/or troubled-cell indicator.
Definition: MinmodType.hpp:20
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
memory
Element.hpp
Mesh< VolumeDim >
DirectionMap
Definition: DirectionMap.hpp:15
cstdlib
limits
Gsl.hpp
Tensor.hpp
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Direction.hpp
unordered_map
type_traits
TMPL.hpp
MeanValue.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13