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