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/NewtonianMhd/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 0 : namespace NewtonianMhd::Limiters::Tci {
40 :
41 : /// \ingroup LimitersGroup
42 : /// \brief Implements the troubled-cell indicator from Krivodonova et al, 2004.
43 : ///
44 : /// The KXRCF (these are the author initials) TCI is described in
45 : /// \cite Krivodonova2004.
46 : ///
47 : /// In summary, this TCI uses the size of discontinuities between neighboring DG
48 : /// elements to determine the smoothness of the solution. This works because the
49 : /// discontinuities converge rapidly for smooth solutions, therefore a large
50 : /// discontinuity suggests a lack of smoothness and the need to apply a limiter.
51 : ///
52 : /// The reference sets the constant we call `kxrcf_constant` to 1. This should
53 : /// generally be a good threshold to use, though it might not be the optimal
54 : /// value (in balancing robustness vs accuracy) for any particular problem.
55 : ///
56 : /// This implementation
57 : /// - does not support h- or p-refinement; this is checked by assertion.
58 : /// - chooses not to check external boundaries, because this adds complexity.
59 : /// However, by not checking external boundaries, the implementation may not
60 : /// be robust for problems that feed in shocks through boundary conditions.
61 : template <size_t VolumeDim, typename PackagedData>
62 1 : bool kxrcf_indicator(
63 : const double kxrcf_constant, const Scalar<DataVector>& cons_mass_density,
64 : const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
65 : const Scalar<DataVector>& cons_energy_density, const Mesh<VolumeDim>& mesh,
66 : const Element<VolumeDim>& element,
67 : const std::array<double, VolumeDim>& element_size,
68 : const Scalar<DataVector>& det_logical_to_inertial_jacobian,
69 : const typename evolution::dg::Tags::NormalCovectorAndMagnitude<
70 : VolumeDim>::type& normals_and_magnitudes,
71 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
72 : boost::hash<DirectionalId<VolumeDim>>>&
73 : neighbor_data) {
74 : // Enforce restrictions on h-refinement, p-refinement
75 : if (UNLIKELY(
76 : alg::any_of(element.neighbors(), [](const auto& direction_neighbors) {
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) {
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.direction();
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<NewtonianMhd::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<NewtonianMhd::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 NewtonianMhd::Limiters::Tci
|
