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 <cstddef>
8 : #include <limits>
9 : #include <string>
10 : #include <type_traits>
11 : #include <unordered_map>
12 : #include <utility>
13 :
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tags.hpp"
16 : #include "DataStructures/Tags/TempTensor.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Evolution/Systems/NewtonianEuler/Characteristics.hpp"
20 : #include "Evolution/Systems/NewtonianEuler/SoundSpeedSquared.hpp"
21 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
22 : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
23 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
24 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
25 : #include "Utilities/Gsl.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 0 : namespace NewtonianEuler::Limiters {
29 :
30 : /*!
31 : * \brief Compute the transform matrices between the conserved variables and
32 : * the characteristic variables of the NewtonianEuler system.
33 : *
34 : * Wraps calls to `NewtonianEuler::right_eigenvectors` and
35 : * `NewtonianEuler::left_eigenvectors`.
36 : */
37 : template <size_t VolumeDim>
38 1 : std::pair<Matrix, Matrix> right_and_left_eigenvectors(
39 : const Scalar<double>& mean_density,
40 : const tnsr::I<double, VolumeDim>& mean_momentum,
41 : const Scalar<double>& mean_energy,
42 : const EquationsOfState::EquationOfState<false, 2>& equation_of_state,
43 : const tnsr::i<double, VolumeDim>& unit_vector,
44 : bool compute_char_transformation_numerically = false);
45 :
46 : /// @{
47 : /// \brief Compute characteristic fields from conserved fields
48 : ///
49 : /// Note that these functions apply the same transformation to every grid point
50 : /// in the element, using the same matrix from `left_eigenvectors`.
51 : /// This is in contrast to the characteristic transformation used in the
52 : /// GeneralizedHarmonic upwind flux, which is computed pointwise.
53 : ///
54 : /// By using a fixed matrix of eigenvectors computed from the cell-averaged
55 : /// fields, we ensure we can consistently transform back from the characteristic
56 : /// variables even after applying the limiter (the limiter changes the pointwise
57 : /// values but preserves the cell averages).
58 : template <size_t VolumeDim>
59 1 : void characteristic_fields(
60 : gsl::not_null<tuples::TaggedTuple<
61 : ::Tags::Mean<NewtonianEuler::Tags::VMinus>,
62 : ::Tags::Mean<NewtonianEuler::Tags::VMomentum<VolumeDim>>,
63 : ::Tags::Mean<NewtonianEuler::Tags::VPlus>>*>
64 : char_means,
65 : const tuples::TaggedTuple<
66 : ::Tags::Mean<NewtonianEuler::Tags::MassDensityCons>,
67 : ::Tags::Mean<NewtonianEuler::Tags::MomentumDensity<VolumeDim>>,
68 : ::Tags::Mean<NewtonianEuler::Tags::EnergyDensity>>& cons_means,
69 : const Matrix& left);
70 :
71 : template <size_t VolumeDim>
72 1 : void characteristic_fields(
73 : gsl::not_null<Scalar<DataVector>*> char_v_minus,
74 : gsl::not_null<tnsr::I<DataVector, VolumeDim>*> char_v_momentum,
75 : gsl::not_null<Scalar<DataVector>*> char_v_plus,
76 : const Scalar<DataVector>& cons_mass_density,
77 : const tnsr::I<DataVector, VolumeDim>& cons_momentum_density,
78 : const Scalar<DataVector>& cons_energy_density, const Matrix& left);
79 :
80 : template <size_t VolumeDim>
81 1 : void characteristic_fields(
82 : gsl::not_null<
83 : Variables<tmpl::list<NewtonianEuler::Tags::VMinus,
84 : NewtonianEuler::Tags::VMomentum<VolumeDim>,
85 : NewtonianEuler::Tags::VPlus>>*>
86 : char_vars,
87 : const Variables<tmpl::list<NewtonianEuler::Tags::MassDensityCons,
88 : NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
89 : NewtonianEuler::Tags::EnergyDensity>>& cons_vars,
90 : const Matrix& left);
91 : /// @}
92 :
93 : /// \brief Compute conserved fields from characteristic fields
94 : ///
95 : /// Note that this function applies the same transformation to every grid point
96 : /// in the element, using the same matrix from `right_eigenvectors`.
97 : /// This is in contrast to the characteristic transformation used in the
98 : /// GeneralizedHarmonic upwind flux, which is computed pointwise.
99 : ///
100 : /// By using a fixed matrix of eigenvectors computed from the cell-averaged
101 : /// fields, we ensure we can consistently transform back from the characteristic
102 : /// variables even after applying the limiter (the limiter changes the pointwise
103 : /// values but preserves the cell averages).
104 : template <size_t VolumeDim>
105 1 : void conserved_fields_from_characteristic_fields(
106 : gsl::not_null<Scalar<DataVector>*> cons_mass_density,
107 : gsl::not_null<tnsr::I<DataVector, VolumeDim>*> cons_momentum_density,
108 : gsl::not_null<Scalar<DataVector>*> cons_energy_density,
109 : const Scalar<DataVector>& char_v_minus,
110 : const tnsr::I<DataVector, VolumeDim>& char_v_momentum,
111 : const Scalar<DataVector>& char_v_plus, const Matrix& right);
112 :
113 : /// \brief Apply a limiter to the characteristic fields computed with respect
114 : /// to each direction in the volume, then take average of results
115 : ///
116 : /// When computing the characteristic fields in the volume to pass as inputs to
117 : /// the limiter, it is not necessarily clear (in more than one dimension) which
118 : /// unit vector should be used in the characteristic decomposition. This is in
119 : /// contrast to uses of characteristic fields for, e.g., computing boundary
120 : /// conditions where the boundary normal provides a clear choice of unit vector.
121 : //
122 : /// The common solution to this challenge is to average the results of applying
123 : /// the limiter with different choices of characteristic decomposition:
124 : /// - For each direction (e.g. in 3D for each of
125 : /// \f$(\hat{x}, \hat{y}, \hat{z})\f$), compute the characteristic fields with
126 : /// respect to this direction. This gives `VolumeDim` sets of characteristic
127 : /// fields.
128 : /// - Apply the limiter to each set of characteristic fields, then convert back
129 : /// to conserved fields. This gives `VolumeDim` different sets of limited
130 : /// conserved fields.
131 : /// - Average the new conserved fields; this is the result of the limiting
132 : /// process.
133 : ///
134 : /// This function handles the logic of computing the different characteristic
135 : /// fields, limiting them, converting back to conserved fields, and averaging.
136 : ///
137 : /// The limiter to apply is passed in via a lambda which is responsible for
138 : /// converting any neighbor data to the characteristic representation, and then
139 : /// applying the limiter to all NewtonianEuler characteristic fields.
140 : template <size_t VolumeDim, typename LimiterLambda>
141 1 : bool apply_limiter_to_characteristic_fields_in_all_directions(
142 : const gsl::not_null<Scalar<DataVector>*> mass_density_cons,
143 : const gsl::not_null<tnsr::I<DataVector, VolumeDim>*> momentum_density,
144 : const gsl::not_null<Scalar<DataVector>*> energy_density,
145 : const Mesh<VolumeDim>& mesh,
146 : const EquationsOfState::EquationOfState<false, 2>& equation_of_state,
147 : const LimiterLambda& prepare_and_apply_limiter,
148 : const bool compute_char_transformation_numerically = false) {
149 : // Temp variables for calculations
150 : // There are quite a few tensors in this allocation because in general we need
151 : // to preserve the input (cons field) tensors until they are overwritten at
152 : // the very end of the computation... then we need additional buffers for the
153 : // char fields, the limited cons fields, and the accumulated cons fields for
154 : // averaging.
155 : //
156 : // Possible optimization: specialize the 1D case which doesn't do average so
157 : // doesn't need as many buffers
158 : Variables<tmpl::list<
159 : NewtonianEuler::Tags::VMinus, NewtonianEuler::Tags::VMomentum<VolumeDim>,
160 : NewtonianEuler::Tags::VPlus, ::Tags::TempScalar<0>,
161 : ::Tags::TempI<0, VolumeDim>, ::Tags::TempScalar<1>, ::Tags::TempScalar<2>,
162 : ::Tags::TempI<1, VolumeDim>, ::Tags::TempScalar<3>>>
163 : temp_buffer(mesh.number_of_grid_points());
164 : auto& char_v_minus = get<NewtonianEuler::Tags::VMinus>(temp_buffer);
165 : auto& char_v_momentum =
166 : get<NewtonianEuler::Tags::VMomentum<VolumeDim>>(temp_buffer);
167 : auto& char_v_plus = get<NewtonianEuler::Tags::VPlus>(temp_buffer);
168 : auto& temp_mass_density_cons = get<::Tags::TempScalar<0>>(temp_buffer);
169 : auto& temp_momentum_density = get<::Tags::TempI<0, VolumeDim>>(temp_buffer);
170 : auto& temp_energy_density = get<::Tags::TempScalar<1>>(temp_buffer);
171 : auto& accumulate_mass_density_cons = get<::Tags::TempScalar<2>>(temp_buffer);
172 : auto& accumulate_momentum_density =
173 : get<::Tags::TempI<1, VolumeDim>>(temp_buffer);
174 : auto& accumulate_energy_density = get<::Tags::TempScalar<3>>(temp_buffer);
175 :
176 : // Initialize the accumulating tensors
177 : get(accumulate_mass_density_cons) = 0.;
178 : for (size_t i = 0; i < VolumeDim; ++i) {
179 : accumulate_momentum_density.get(i) = 0.;
180 : }
181 : get(accumulate_energy_density) = 0.;
182 :
183 : // Cellwise means, used in computing the cons/char transformations
184 : const auto mean_density =
185 : Scalar<double>{mean_value(get(*mass_density_cons), mesh)};
186 : const auto mean_momentum = [&momentum_density, &mesh]() {
187 : tnsr::I<double, VolumeDim> result{};
188 : for (size_t i = 0; i < VolumeDim; ++i) {
189 : result.get(i) = mean_value(momentum_density->get(i), mesh);
190 : }
191 : return result;
192 : }();
193 : const auto mean_energy =
194 : Scalar<double>{mean_value(get(*energy_density), mesh)};
195 :
196 : bool some_component_was_limited = false;
197 :
198 : // Loop over directions, then compute chars w.r.t. this direction and limit
199 : for (size_t d = 0; d < VolumeDim; ++d) {
200 : const auto unit_vector = [&d]() {
201 : auto components = make_array<VolumeDim>(0.);
202 : components[d] = 1.;
203 : return tnsr::i<double, VolumeDim>(components);
204 : }();
205 : const auto right_and_left =
206 : NewtonianEuler::Limiters::right_and_left_eigenvectors(
207 : mean_density, mean_momentum, mean_energy, equation_of_state,
208 : unit_vector, compute_char_transformation_numerically);
209 : const auto& right = right_and_left.first;
210 : const auto& left = right_and_left.second;
211 :
212 : // Transform tensors to characteristics
213 : NewtonianEuler::Limiters::characteristic_fields(
214 : make_not_null(&char_v_minus), make_not_null(&char_v_momentum),
215 : make_not_null(&char_v_plus), *mass_density_cons, *momentum_density,
216 : *energy_density, left);
217 :
218 : // Transform neighbor data and apply limiter
219 : const bool some_component_was_limited_with_this_unit_vector =
220 : prepare_and_apply_limiter(make_not_null(&char_v_minus),
221 : make_not_null(&char_v_momentum),
222 : make_not_null(&char_v_plus), left);
223 :
224 : some_component_was_limited =
225 : some_component_was_limited_with_this_unit_vector or
226 : some_component_was_limited;
227 :
228 : // Transform back to conserved variables. But skip the transformation if no
229 : // limiting occured with this unit vector.
230 : if (some_component_was_limited_with_this_unit_vector) {
231 : NewtonianEuler::Limiters::conserved_fields_from_characteristic_fields(
232 : make_not_null(&temp_mass_density_cons),
233 : make_not_null(&temp_momentum_density),
234 : make_not_null(&temp_energy_density), char_v_minus, char_v_momentum,
235 : char_v_plus, right);
236 : } else {
237 : temp_mass_density_cons = *mass_density_cons;
238 : temp_momentum_density = *momentum_density;
239 : temp_energy_density = *energy_density;
240 : }
241 :
242 : // Add to running sum for averaging
243 : const double one_over_dim = 1.0 / static_cast<double>(VolumeDim);
244 : get(accumulate_mass_density_cons) +=
245 : one_over_dim * get(temp_mass_density_cons);
246 : for (size_t i = 0; i < VolumeDim; ++i) {
247 : accumulate_momentum_density.get(i) +=
248 : one_over_dim * temp_momentum_density.get(i);
249 : }
250 : get(accumulate_energy_density) += one_over_dim * get(temp_energy_density);
251 : } // for loop over dimensions
252 :
253 : *mass_density_cons = accumulate_mass_density_cons;
254 : *momentum_density = accumulate_momentum_density;
255 : *energy_density = accumulate_energy_density;
256 : return some_component_was_limited;
257 : }
258 :
259 : } // namespace NewtonianEuler::Limiters
|
