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 "ErrorHandling/Assert.hpp"
25 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
29 
30 /// \cond
31 template <size_t VolumeDim>
32 class ElementId;
33 /// \endcond
34 
35 namespace NewtonianEuler {
36 namespace Limiters {
37 namespace Tci {
38 
39 /// \ingroup LimitersGroup
40 /// \brief Implements the troubled-cell indicator from Krivodonova et al, 2004.
41 ///
42 /// The KXRCF (these are the author initials) TCI is described in
43 /// \cite Krivodonova2004.
44 ///
45 /// In summary, this TCI uses the size of discontinuities between neighboring DG
46 /// elements to determine the smoothness of the solution. This works because the
47 /// discontinuities converge rapidly for smooth solutions, therefore a large
48 /// discontinuity suggests a lack of smoothness and the need to apply a limiter.
49 ///
50 /// The reference sets the constant we call `kxrcf_constant` to 1. This should
51 /// generally be a good threshold to use, though it might not be the optimal
52 /// value (in balancing robustness vs accuracy) for any particular problem.
53 ///
54 /// This implementation
55 /// - does not support h- or p-refinement; this is checked by assertion.
56 /// - chooses not to check external boundaries, because this adds complexity.
57 /// However, by not checking external boundaries, the implementation may not
58 /// be robust for problems that feed in shocks through boundary conditions.
59 template <size_t VolumeDim, typename PackagedData>
61  const double kxrcf_constant, const Scalar<DataVector>& cons_mass_density,
62  const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
63  const Scalar<DataVector>& cons_energy_density, const Mesh<VolumeDim>& mesh,
64  const Element<VolumeDim>& element,
65  const std::array<double, VolumeDim>& element_size,
66  const Scalar<DataVector>& det_logical_to_inertial_jacobian,
68  tnsr::i<DataVector, VolumeDim>>&
69  unnormalized_normals,
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  for (const auto& [neighbor, data] : neighbor_data) {
110  // Skip computations on boundary if the boundary is external. This choice
111  // might be problematic for evolutions (likely only simple test cases) that
112  // feed in shocks through the boundary condition: the limiter might fail to
113  // activate in the cell that touches the boundary.
114  //
115  // Note that to do this properly we would need the limiter to know about the
116  // boundary condition in general, which may be difficult. Furthermore, such
117  // a change may also require changing the tags at the call site to ensure
118  // that ALL face normals are grabbed from the databox (and not only the
119  // internal face normals, as occurs when grabbing
120  // `Tags::Interface<Tags::InternalDirections, ...>`).
121  const auto& dir = neighbor.first;
122  if (unnormalized_normals.find(dir) == unnormalized_normals.end()) {
123  continue;
124  }
125 
126  const size_t sliced_dim = dir.dimension();
127  const size_t index_of_slice =
128  (dir.side() == Side::Lower ? 0 : mesh.extents()[sliced_dim] - 1);
129  const auto momentum_on_slice = data_on_slice(
130  cons_momentum_density, mesh.extents(), sliced_dim, index_of_slice);
131  const auto momentum_dot_normal =
132  dot_product(momentum_on_slice, unnormalized_normals.at(dir));
133 
134  // Skip boundaries with no significant inflow
135  // Note: the cutoff value here is small but arbitrarily chosen.
136  if (min(get(momentum_dot_normal)) > -1e-12) {
137  continue;
138  }
139  inflow_boundaries_present = true;
140 
141  // This mask has value 1. for momentum_dot_normal < 0.
142  // 0. for momentum_dot_normal >= 0.
143  const DataVector inflow_mask = 1. - step_function(get(momentum_dot_normal));
144  // Mask is then weighted pointwise by the Jacobian determinant giving
145  // surface integrals in inertial coordinates. This Jacobian determinant is
146  // given by the product of the volume Jacobian determinant with the norm of
147  // the unnormalized face normals.
148  const DataVector weighted_inflow_mask =
149  inflow_mask *
150  get(data_on_slice(det_logical_to_inertial_jacobian, mesh.extents(),
151  sliced_dim, index_of_slice)) *
152  get(magnitude(unnormalized_normals.at(dir)));
153 
154  inflow_area +=
155  definite_integral(weighted_inflow_mask, mesh.slice_away(sliced_dim));
156 
157  // This is the step that is incompatible with h/p refinement. For use with
158  // h/p refinement, would need to correctly obtain the neighbor solution on
159  // the local grid points.
160  const auto neighbor_vars_on_slice = data_on_slice(
161  data.volume_data, mesh.extents(), sliced_dim, index_of_slice);
162 
163  const auto density_on_slice = data_on_slice(
164  cons_mass_density, mesh.extents(), sliced_dim, index_of_slice);
165  const auto& neighbor_density_on_slice =
166  get<NewtonianEuler::Tags::MassDensityCons>(neighbor_vars_on_slice);
167  inflow_delta_density += definite_integral(
168  (get(density_on_slice) - get(neighbor_density_on_slice)) *
169  weighted_inflow_mask,
170  mesh.slice_away(sliced_dim));
171 
172  const auto energy_on_slice = data_on_slice(
173  cons_energy_density, mesh.extents(), sliced_dim, index_of_slice);
174  const auto& neighbor_energy_on_slice =
175  get<NewtonianEuler::Tags::EnergyDensity>(neighbor_vars_on_slice);
176  inflow_delta_energy += definite_integral(
177  (get(energy_on_slice) - get(neighbor_energy_on_slice)) *
178  weighted_inflow_mask,
179  mesh.slice_away(sliced_dim));
180  }
181 
182  if (not inflow_boundaries_present) {
183  // No boundaries had inflow, so not a troubled cell
184  return false;
185  }
186 
187  // KXRCF take h to be the radius of the circumscribed circle
188  const double h = 0.5 * magnitude(element_size);
189  const double h_pow = pow(h, 0.5 * mesh.extents(0));
190 
191  ASSERT(inflow_area > 0.,
192  "Sanity check failed: negative area of inflow boundaries");
193 
194  const double norm_squared_density = mean_value(
195  get(det_logical_to_inertial_jacobian) * square(get(cons_mass_density)),
196  mesh);
197  ASSERT(norm_squared_density > 0.,
198  "Sanity check failed: negative density norm over element");
199  const double norm_density = sqrt(norm_squared_density);
200  const double ratio_for_density =
201  abs(inflow_delta_density) / (h_pow * inflow_area * norm_density);
202 
203  const double norm_squared_energy = mean_value(
204  get(det_logical_to_inertial_jacobian) * square(get(cons_energy_density)),
205  mesh);
206  ASSERT(norm_squared_energy > 0.,
207  "Sanity check failed: negative energy norm over element");
208  const double norm_energy = sqrt(norm_squared_energy);
209  const double ratio_for_energy =
210  abs(inflow_delta_energy) / (h_pow * inflow_area * norm_energy);
211 
212  return (ratio_for_density > kxrcf_constant or
213  ratio_for_energy > kxrcf_constant);
214 }
215 
216 } // namespace Tci
217 } // namespace Limiters
218 } // namespace NewtonianEuler
Limiters
Things relating to limiting.
Definition: HwenoImpl.cpp:111
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
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
std::pair
Tags.hpp
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
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:24
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
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:42
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
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
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
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 std::unordered_map< Direction< VolumeDim >, tnsr::i< DataVector, VolumeDim >> &unnormalized_normals, 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:60
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
Mesh.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
string