KxrcfTci.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 <cstddef>
9 #include <string>
10 #include <unordered_map>
11 #include <utility>
12 
13 #include "DataStructures/SliceVariables.hpp"
14 #include "DataStructures/Tags.hpp"
16 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
17 #include "DataStructures/Tensor/Slice.hpp"
20 #include "Domain/SizeOfElement.hpp"
23 #include "Domain/Tags.hpp"
24 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
25 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
30 
31 /// \cond
32 template <size_t VolumeDim>
33 class ElementId;
34 /// \endcond
35 
36 namespace NewtonianEuler {
37 namespace Limiters {
38 namespace Tci {
39 
40 /// \ingroup LimitersGroup
41 /// \brief Implements the troubled-cell indicator from Krivodonova et al, 2004.
42 ///
43 /// The KXRCF (these are the author initials) TCI is described in
44 /// \cite Krivodonova2004.
45 ///
46 /// In summary, this TCI uses the size of discontinuities between neighboring DG
47 /// elements to determine the smoothness of the solution. This works because the
48 /// discontinuities converge rapidly for smooth solutions, therefore a large
49 /// discontinuity suggests a lack of smoothness and the need to apply a limiter.
50 ///
51 /// The reference sets the constant we call `kxrcf_constant` to 1. This should
52 /// generally be a good threshold to use, though it might not be the optimal
53 /// value (in balancing robustness vs accuracy) for any particular problem.
54 ///
55 /// This implementation
56 /// - does not support h- or p-refinement; this is checked by assertion.
57 /// - chooses not to check external boundaries, because this adds complexity.
58 /// However, by not checking external boundaries, the implementation may not
59 /// be robust for problems that feed in shocks through boundary conditions.
60 template <size_t VolumeDim, typename PackagedData>
62  const double kxrcf_constant, const Scalar<DataVector>& cons_mass_density,
63  const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
64  const Scalar<DataVector>& cons_energy_density, const Mesh<VolumeDim>& mesh,
65  const Element<VolumeDim>& element,
66  const std::array<double, VolumeDim>& element_size,
67  const Scalar<DataVector>& det_logical_to_inertial_jacobian,
69  VolumeDim>::type& normals_and_magnitudes,
70  const std::unordered_map<
73  neighbor_data) noexcept {
74  // Enforce restrictions on h-refinement, p-refinement
75  if (UNLIKELY(alg::any_of(element.neighbors(),
76  [](const auto& direction_neighbors) noexcept {
77  return direction_neighbors.second.size() != 1;
78  }))) {
79  ERROR("The Kxrcf TCI does not yet support h-refinement");
80  // Removing this limitation will require adapting the surface integrals to
81  // correctly acount for,
82  // - multiple (smaller) neighbors contributing to the integral
83  // - only a portion of a (larger) neighbor contributing to the integral
84  }
85  alg::for_each(neighbor_data, [&mesh](const auto& neighbor_and_data) noexcept {
86  if (UNLIKELY(neighbor_and_data.second.mesh != mesh)) {
87  ERROR("The Kxrcf TCI does not yet support p-refinement");
88  // Removing this limitation will require generalizing the surface
89  // integrals to make sure the meshes are consistent.
90  }
91  });
92  // Check the mesh matches expectations:
93  // - the extents must be uniform, because the TCI expects a unique "order" for
94  // the DG scheme. This shows up when computing the threshold parameter,
95  // h^(degree+1)/2
96  // - the quadrature must be GL, because the current implementation doesn't
97  // handle extrapolation to the boundary, though this could be changed.
98  // - the basis in principle could be changed, but until we need to change it
99  // we just check it's the expected Legendre for simplicity
100  ASSERT(mesh == Mesh<VolumeDim>(mesh.extents(0), Spectral::Basis::Legendre,
101  Spectral::Quadrature::GaussLobatto),
102  "The Kxrcf TCI expects a uniform LGL mesh, but got mesh = " << mesh);
103 
104  bool inflow_boundaries_present = false;
105  double inflow_area = 0.;
106  double inflow_delta_density = 0.;
107  double inflow_delta_energy = 0.;
108 
109  // Skip boundary integrations on external boundaries. This choice might be
110  // problematic for evolutions (likely only simple test cases) that feed in
111  // shocks through the boundary condition: the limiter might fail to activate
112  // in the cell that touches the boundary.
113  //
114  // To properly compute the limiter at external boundaries we would need the
115  // limiter to know about the boundary condition, which may be difficult to
116  // do in a general way.
117  for (const auto& [neighbor, data] : neighbor_data) {
118  const auto& dir = neighbor.first;
119 
120  // Check consistency of neighbor_data with element and normals
121  ASSERT(element.neighbors().contains(dir),
122  "Received neighbor data from dir = "
123  << dir << ", but element has no neighbor in this dir");
124  ASSERT(normals_and_magnitudes.contains(dir),
125  "Received neighbor data from dir = "
126  << dir
127  << ", but normals_and_magnitudes has no normal in this dir");
128  ASSERT(normals_and_magnitudes.at(dir).has_value(),
129  "The normals_and_magnitudes are not up-to-date in dir = " << dir);
130  const auto& normal = get<evolution::dg::Tags::NormalCovector<VolumeDim>>(
131  normals_and_magnitudes.at(dir).value());
132  const auto& magnitude_of_normal =
133  get<evolution::dg::Tags::MagnitudeOfNormal>(
134  normals_and_magnitudes.at(dir).value());
135 
136  const size_t sliced_dim = dir.dimension();
137  const size_t index_of_slice =
138  (dir.side() == Side::Lower ? 0 : mesh.extents()[sliced_dim] - 1);
139  const auto momentum_on_slice = data_on_slice(
140  cons_momentum_density, mesh.extents(), sliced_dim, index_of_slice);
141  const auto momentum_dot_normal = dot_product(momentum_on_slice, normal);
142 
143  // Skip boundaries with no significant inflow
144  // Note: the cutoff value here is small but arbitrarily chosen.
145  if (min(get(momentum_dot_normal)) > -1e-12) {
146  continue;
147  }
148  inflow_boundaries_present = true;
149 
150  // This mask has value 1. for momentum_dot_normal < 0.
151  // 0. for momentum_dot_normal >= 0.
152  const DataVector inflow_mask = 1. - step_function(get(momentum_dot_normal));
153  // Mask is then weighted pointwise by the Jacobian determinant giving
154  // surface integrals in inertial coordinates. This Jacobian determinant is
155  // given by the product of the volume Jacobian determinant with the
156  // magnitude of the unnormalized normal covectors.
157  const DataVector weighted_inflow_mask =
158  inflow_mask *
159  get(data_on_slice(det_logical_to_inertial_jacobian, mesh.extents(),
160  sliced_dim, index_of_slice)) *
161  get(magnitude_of_normal);
162 
163  inflow_area +=
164  definite_integral(weighted_inflow_mask, mesh.slice_away(sliced_dim));
165 
166  // This is the step that is incompatible with h/p refinement. For use with
167  // h/p refinement, would need to correctly obtain the neighbor solution on
168  // the local grid points.
169  const auto neighbor_vars_on_slice = data_on_slice(
170  data.volume_data, mesh.extents(), sliced_dim, index_of_slice);
171 
172  const auto density_on_slice = data_on_slice(
173  cons_mass_density, mesh.extents(), sliced_dim, index_of_slice);
174  const auto& neighbor_density_on_slice =
175  get<NewtonianEuler::Tags::MassDensityCons>(neighbor_vars_on_slice);
176  inflow_delta_density += definite_integral(
177  (get(density_on_slice) - get(neighbor_density_on_slice)) *
178  weighted_inflow_mask,
179  mesh.slice_away(sliced_dim));
180 
181  const auto energy_on_slice = data_on_slice(
182  cons_energy_density, mesh.extents(), sliced_dim, index_of_slice);
183  const auto& neighbor_energy_on_slice =
184  get<NewtonianEuler::Tags::EnergyDensity>(neighbor_vars_on_slice);
185  inflow_delta_energy += definite_integral(
186  (get(energy_on_slice) - get(neighbor_energy_on_slice)) *
187  weighted_inflow_mask,
188  mesh.slice_away(sliced_dim));
189  }
190 
191  if (not inflow_boundaries_present) {
192  // No boundaries had inflow, so not a troubled cell
193  return false;
194  }
195 
196  // KXRCF take h to be the radius of the circumscribed circle
197  const double h = 0.5 * magnitude(element_size);
198  const double h_pow = pow(h, 0.5 * mesh.extents(0));
199 
200  ASSERT(inflow_area > 0.,
201  "Sanity check failed: negative area of inflow boundaries");
202 
203  const double norm_squared_density = mean_value(
204  get(det_logical_to_inertial_jacobian) * square(get(cons_mass_density)),
205  mesh);
206  ASSERT(norm_squared_density > 0.,
207  "Sanity check failed: negative density norm over element");
208  const double norm_density = sqrt(norm_squared_density);
209  const double ratio_for_density =
210  abs(inflow_delta_density) / (h_pow * inflow_area * norm_density);
211 
212  const double norm_squared_energy = mean_value(
213  get(det_logical_to_inertial_jacobian) * square(get(cons_energy_density)),
214  mesh);
215  ASSERT(norm_squared_energy > 0.,
216  "Sanity check failed: negative energy norm over element");
217  const double norm_energy = sqrt(norm_squared_energy);
218  const double ratio_for_energy =
219  abs(inflow_delta_energy) / (h_pow * inflow_area * norm_energy);
220 
221  return (ratio_for_density > kxrcf_constant or
222  ratio_for_energy > kxrcf_constant);
223 }
224 
225 } // namespace Tci
226 } // namespace Limiters
227 } // namespace NewtonianEuler
Limiters
Things relating to limiting.
Definition: HwenoImpl.hpp:42
alg::any_of
decltype(auto) any_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::any_of.
Definition: Algorithm.hpp:190
pow
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:83
std::pair
evolution::dg::Tags::NormalCovectorAndMagnitude
The normal covector and its magnitude for all internal faces of an element.
Definition: NormalVectorTags.hpp:35
Tags.hpp
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
NewtonianEuler
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
dot_product
void dot_product(const gsl::not_null< Scalar< DataType > * > dot_product, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean dot product of two vectors or one forms.
Definition: DotProduct.hpp:25
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
DotProduct.hpp
db::data_on_slice
Variables< tmpl::list< TagsToSlice... > > data_on_slice(const db::DataBox< TagsList > &box, const Index< VolumeDim > &element_extents, const size_t sliced_dim, const size_t fixed_index, tmpl::list< TagsToSlice... >) noexcept
Slices volume Tensors from a DataBox into a Variables
Definition: DataOnSlice.hpp:33
cstddef
Assert.hpp
array
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
mean_value
double mean_value(const DataVector &f, const Mesh< Dim > &mesh) noexcept
Compute the mean value of a function over a manifold.
Definition: MeanValue.hpp:47
alg::for_each
decltype(auto) for_each(const Container &c, UnaryFunction &&f)
Convenience wrapper around std::for_each, returns the result of std::for_each(begin(c),...
Definition: Algorithm.hpp:291
NewtonianEuler::Limiters::Tci::kxrcf_indicator
bool kxrcf_indicator(const double kxrcf_constant, const Scalar< DataVector > &cons_mass_density, const tnsr::I< DataVector, VolumeDim > &cons_momentum_density, const Scalar< DataVector > &cons_energy_density, const Mesh< VolumeDim > &mesh, const Element< VolumeDim > &element, const std::array< double, VolumeDim > &element_size, const Scalar< DataVector > &det_logical_to_inertial_jacobian, const typename evolution::dg::Tags::NormalCovectorAndMagnitude< VolumeDim >::type &normals_and_magnitudes, const std::unordered_map< std::pair< Direction< VolumeDim >, ElementId< VolumeDim >>, PackagedData, boost::hash< std::pair< Direction< VolumeDim >, ElementId< VolumeDim >>>> &neighbor_data) noexcept
Implements the troubled-cell indicator from Krivodonova et al, 2004.
Definition: KxrcfTci.hpp:61
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
definite_integral
double definite_integral(const DataVector &integrand, const Mesh< Dim > &mesh) noexcept
Compute the definite integral of a function over a manifold.
Variables.hpp
Element.hpp
Mesh< VolumeDim >
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Tensor.hpp
magnitude
Scalar< DataType > magnitude(const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector) noexcept
Compute the Euclidean magnitude of a rank-1 tensor.
Definition: Magnitude.hpp:27
DefiniteIntegral.hpp
Direction.hpp
unordered_map
step_function
constexpr T step_function(const T &arg) noexcept
Defines the Heaviside step function for arithmetic types. .
Definition: Math.hpp:71
MeanValue.hpp
Mesh.hpp
string