7 #include <boost/functional/hash.hpp>
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"
24 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
31 template <
size_t VolumeDim>
59 template <
size_t VolumeDim,
typename PackagedData>
62 const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
68 tnsr::i<DataVector, VolumeDim>>&
73 neighbor_data) noexcept {
76 [](
const auto& direction_neighbors) noexcept {
77 return direction_neighbors.second.size() != 1;
79 ERROR(
"The Kxrcf TCI does not yet support h-refinement");
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");
101 Spectral::Quadrature::GaussLobatto),
102 "The Kxrcf TCI expects a uniform LGL mesh, but got mesh = " << mesh);
104 bool inflow_boundaries_present =
false;
105 double inflow_area = 0.;
106 double inflow_delta_density = 0.;
107 double inflow_delta_energy = 0.;
109 for (
const auto& [neighbor, data] : neighbor_data) {
121 const auto& dir = neighbor.first;
122 if (unnormalized_normals.find(dir) == unnormalized_normals.end()) {
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);
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));
136 if (min(
get(momentum_dot_normal)) > -1e-12) {
139 inflow_boundaries_present =
true;
151 sliced_dim, index_of_slice)) *
161 data.volume_data, mesh.extents(), sliced_dim, index_of_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);
168 (
get(density_on_slice) -
get(neighbor_density_on_slice)) *
169 weighted_inflow_mask,
170 mesh.slice_away(sliced_dim));
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);
177 (
get(energy_on_slice) -
get(neighbor_energy_on_slice)) *
178 weighted_inflow_mask,
179 mesh.slice_away(sliced_dim));
182 if (not inflow_boundaries_present) {
188 const double h = 0.5 *
magnitude(element_size);
189 const double h_pow =
pow(h, 0.5 * mesh.extents(0));
192 "Sanity check failed: negative area of inflow boundaries");
194 const double norm_squared_density =
mean_value(
195 get(det_logical_to_inertial_jacobian) *
square(
get(cons_mass_density)),
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);
203 const double norm_squared_energy =
mean_value(
204 get(det_logical_to_inertial_jacobian) *
square(
get(cons_energy_density)),
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);
212 return (ratio_for_density > kxrcf_constant or
213 ratio_for_energy > kxrcf_constant);
Things relating to limiting.
Definition: HwenoImpl.cpp:111
decltype(auto) any_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::any_of.
Definition: Algorithm.hpp:190
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
#define UNLIKELY(x)
Definition: Gsl.hpp:73
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
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:88
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8
Definition: Direction.hpp:23
Definition: Element.hpp:29
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
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
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
Stores a collection of function values.
Definition: DataVector.hpp:46
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
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
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
double definite_integral(const DataVector &integrand, const Mesh< Dim > &mesh) noexcept
Compute the definite integral of a function over a manifold.
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
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
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
constexpr T step_function(const T &arg) noexcept
Defines the Heaviside step function for arithmetic types. .
Definition: Math.hpp:71