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