Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <string>
8 :
9 : #include "DataStructures/DataBox/Tag.hpp"
10 : #include "DataStructures/Tensor/EagerMath/Determinant.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Domain/ElementMap.hpp"
13 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
14 : #include "Utilities/SetNumberOfGridPoints.hpp"
15 :
16 : /// \cond
17 : namespace Tags {
18 : struct Time;
19 : } // namespace Tags
20 : /// \endcond
21 :
22 1 : namespace evolution::dg::subcell::fd::Tags {
23 : /// \brief The inverse Jacobian from the element logical frame to the grid frame
24 : /// at the cell centers.
25 : ///
26 : /// Specifically, \f$\partial x^{\bar{i}} / \partial x^i\f$, where \f$\bar{i}\f$
27 : /// denotes the element logical frame and \f$i\f$ denotes the grid frame.
28 : template <size_t Dim>
29 1 : struct InverseJacobianLogicalToGrid : db::SimpleTag {
30 0 : static std::string name() { return "InverseJacobian(Logical,Grid)"; }
31 0 : using type =
32 : ::InverseJacobian<DataVector, Dim, Frame::ElementLogical, Frame::Grid>;
33 : };
34 :
35 : /// \brief The determinant of the inverse Jacobian from the element logical
36 : /// frame to the grid frame at the cell centers.
37 1 : struct DetInverseJacobianLogicalToGrid : db::SimpleTag {
38 0 : static std::string name() { return "Det(InverseJacobian(Logical,Grid))"; }
39 0 : using type = Scalar<DataVector>;
40 : };
41 :
42 : /// \brief The determinant of the inverse Jacobian from the element logical
43 : /// frame to the inertial frame at the cell centers.
44 1 : struct DetInverseJacobianLogicalToInertial : db::SimpleTag {
45 0 : static std::string name() { return "Det(InverseJacobian(Logical,Inertial))"; }
46 0 : using type = Scalar<DataVector>;
47 : };
48 :
49 : /// \brief The inverse Jacobian from the element logical frame to the inertial
50 : /// frame at the cell centers.
51 : ///
52 : /// Specifically, \f$\partial x^{\bar{i}} / \partial x^i\f$, where \f$\bar{i}\f$
53 : /// denotes the element logical frame and \f$i\f$ denotes the inertial frame.
54 : template <size_t Dim>
55 1 : struct InverseJacobianLogicalToInertial : db::SimpleTag {
56 0 : static std::string name() { return "InverseJacobian(Logical,Inertial)"; }
57 0 : using type = ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
58 : Frame::Inertial>;
59 : };
60 :
61 : /// \brief The inverse Hessian from the element logical frame to the grid frame
62 : /// at the cell centers.
63 : ///
64 : /// Specifically, \f$\partial^2 x^{\bar{i}} / {\partial x^i \partial x^j}\f$,
65 : /// where \f$\bar{i}\f$ denotes the element logical frame and \f$i\f$ denotes
66 : /// the grid frame.
67 : template <size_t Dim>
68 1 : struct InverseHessianLogicalToGrid : db::SimpleTag {
69 0 : static std::string name() { return "InverseHessian(Logical,Grid)"; }
70 0 : using type =
71 : ::InverseHessian<DataVector, Dim, Frame::ElementLogical, Frame::Grid>;
72 : };
73 :
74 : /// \brief The inverse Hessian from the element logical frame to the inertial
75 : /// frame at the cell centers.
76 : ///
77 : /// Specifically, \f$\partial^2 x^{\bar{i}} / {\partial x^i \partial x^j}\f$,
78 : /// where \f$\bar{i}\f$ denotes the element logical frame and \f$i\f$ denotes
79 : /// the inertial frame.
80 : template <size_t Dim>
81 1 : struct InverseHessianLogicalToInertial : db::SimpleTag {
82 0 : static std::string name() { return "InverseHessian(Logical,Inertial)"; }
83 0 : using type =
84 : ::InverseHessian<DataVector, Dim, Frame::ElementLogical, Frame::Inertial>;
85 : };
86 :
87 : /// Compute item for the inverse jacobian matrix from logical to
88 : /// grid coordinates
89 : template <typename MapTagLogicalToGrid, size_t Dim>
90 1 : struct InverseJacobianLogicalToGridCompute : InverseJacobianLogicalToGrid<Dim>,
91 : db::ComputeTag {
92 0 : static constexpr size_t dim = Dim;
93 0 : using base = InverseJacobianLogicalToGrid<Dim>;
94 0 : using return_type = typename base::type;
95 0 : using argument_tags = tmpl::list<
96 : MapTagLogicalToGrid,
97 : evolution::dg::subcell::Tags::Coordinates<dim, Frame::ElementLogical>>;
98 0 : static void function(
99 : const gsl::not_null<return_type*> inverse_jacobian_logical_to_grid,
100 : const ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
101 : const tnsr::I<DataVector, dim, Frame::ElementLogical>& logical_coords) {
102 : *inverse_jacobian_logical_to_grid =
103 : logical_to_grid_map.inv_jacobian(logical_coords);
104 : }
105 : };
106 :
107 : /// Compute item for the determinant of the inverse jacobian matrix
108 : /// from logical to grid coordinates
109 : template <size_t Dim>
110 1 : struct DetInverseJacobianLogicalToGridCompute
111 : : DetInverseJacobianLogicalToGrid,
112 : db::ComputeTag {
113 0 : static constexpr size_t dim = Dim;
114 0 : using base = DetInverseJacobianLogicalToGrid;
115 0 : using return_type = typename base::type;
116 0 : using argument_tags = tmpl::list<InverseJacobianLogicalToGrid<Dim>>;
117 0 : static void function(
118 : const gsl::not_null<return_type*> det_inverse_jacobian_logical_to_grid,
119 : const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
120 : Frame::Grid>& inverse_jacobian_logical_to_grid) {
121 : *det_inverse_jacobian_logical_to_grid =
122 : determinant(inverse_jacobian_logical_to_grid);
123 : }
124 : };
125 :
126 : /// Compute item for the inverse jacobian matrix from logical to
127 : /// inertial coordinates
128 : template <typename MapTagGridToInertial, size_t Dim>
129 1 : struct InverseJacobianLogicalToInertialCompute
130 : : InverseJacobianLogicalToInertial<Dim>,
131 : db::ComputeTag {
132 0 : static constexpr size_t dim = Dim;
133 0 : using base = InverseJacobianLogicalToInertial<Dim>;
134 0 : using return_type = typename base::type;
135 0 : using argument_tags =
136 : tmpl::list<MapTagGridToInertial,
137 : evolution::dg::subcell::Tags::Coordinates<dim, Frame::Grid>,
138 : ::Tags::Time, ::domain::Tags::FunctionsOfTime,
139 : InverseJacobianLogicalToGrid<Dim>>;
140 0 : static void function(
141 : const gsl::not_null<return_type*> inverse_jacobian_logical_to_inertial,
142 : const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
143 : grid_to_inertial_map,
144 : const tnsr::I<DataVector, dim, Frame::Grid>& grid_coords,
145 : const double time,
146 : const std::unordered_map<
147 : std::string,
148 : std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
149 : functions_of_time,
150 : const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
151 : Frame::Grid>& inverse_jacobian_logical_to_grid) {
152 : if (grid_to_inertial_map.is_identity()) {
153 : // Optimization for time-independent maps; we just point to the
154 : // logical-to-grid inverse jacobian.
155 : const size_t num_pts = inverse_jacobian_logical_to_grid[0].size();
156 : for (size_t storage_index = 0;
157 : storage_index < inverse_jacobian_logical_to_grid.size();
158 : ++storage_index) {
159 : make_const_view(
160 : make_not_null(&std::as_const(
161 : (*inverse_jacobian_logical_to_inertial)[storage_index])),
162 : inverse_jacobian_logical_to_grid[storage_index], 0, num_pts);
163 : }
164 : } else {
165 : set_number_of_grid_points(inverse_jacobian_logical_to_inertial,
166 : get<0>(grid_coords).size());
167 : // Get grid to inertial inverse jacobian
168 : const auto& inverse_jacobian_grid_to_inertial =
169 : grid_to_inertial_map.inv_jacobian(grid_coords, time,
170 : functions_of_time);
171 : for (size_t i = 0; i < Dim; i++) {
172 : for (size_t j = 0; j < Dim; j++) {
173 : auto& inv_jacobian_component =
174 : inverse_jacobian_logical_to_inertial->get(i, j);
175 : inv_jacobian_component = 0.;
176 : for (size_t k = 0; k < Dim; k++) {
177 : inv_jacobian_component +=
178 : inverse_jacobian_logical_to_grid.get(i, k) *
179 : inverse_jacobian_grid_to_inertial.get(k, j);
180 : }
181 : }
182 : }
183 : }
184 : }
185 : };
186 :
187 : /// Compute item for the determinant of the inverse jacobian matrix
188 : /// from logical to inertial coordinates
189 : template <typename MapTagGridToInertial, size_t Dim>
190 1 : struct DetInverseJacobianLogicalToInertialCompute
191 : : DetInverseJacobianLogicalToInertial,
192 : db::ComputeTag {
193 0 : static constexpr size_t dim = Dim;
194 0 : using base = DetInverseJacobianLogicalToInertial;
195 0 : using return_type = typename base::type;
196 0 : using argument_tags =
197 : tmpl::list<MapTagGridToInertial, InverseJacobianLogicalToInertial<Dim>,
198 : DetInverseJacobianLogicalToGrid>;
199 0 : static void function(
200 : const gsl::not_null<return_type*>
201 : det_inverse_jacobian_logical_to_inertial,
202 : const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
203 : grid_to_inertial_map,
204 : const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
205 : Frame::Inertial>&
206 : inverse_jacobian_logical_to_inertial,
207 : const Scalar<DataVector>& det_inverse_jacobian_logical_to_grid) {
208 : if (grid_to_inertial_map.is_identity()) {
209 : const size_t num_pts = get(det_inverse_jacobian_logical_to_grid).size();
210 : make_const_view(make_not_null(&std::as_const(
211 : get(*det_inverse_jacobian_logical_to_inertial))),
212 : get(det_inverse_jacobian_logical_to_grid), 0, num_pts);
213 : } else {
214 : *det_inverse_jacobian_logical_to_inertial =
215 : determinant(inverse_jacobian_logical_to_inertial);
216 : }
217 : }
218 : };
219 :
220 : #ifdef SPECTRE_AUTODIFF
221 : /// Compute item for the inverse hessian matrix from logical to
222 : /// grid coordinates
223 : template <typename MapTagLogicalToGrid, size_t Dim>
224 1 : struct InverseHessianLogicalToGridCompute : InverseHessianLogicalToGrid<Dim>,
225 : db::ComputeTag {
226 0 : static constexpr size_t dim = Dim;
227 0 : using base = InverseHessianLogicalToGrid<Dim>;
228 0 : using return_type = typename base::type;
229 0 : using argument_tags = tmpl::list<
230 : MapTagLogicalToGrid,
231 : evolution::dg::subcell::Tags::Coordinates<dim, Frame::ElementLogical>>;
232 0 : static void function(
233 : const gsl::not_null<return_type*> inverse_hessian_logical_to_grid,
234 : const ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
235 : const tnsr::I<DataVector, dim, Frame::ElementLogical>& logical_coords) {
236 : *inverse_hessian_logical_to_grid =
237 : logical_to_grid_map.inv_hessian(logical_coords);
238 : }
239 : };
240 :
241 : /// Compute item for the inverse hessian matrix from logical to
242 : /// inertial coordinates
243 : template <typename MapTagGridToInertial, size_t Dim>
244 1 : struct InverseHessianLogicalToInertialCompute
245 : : InverseHessianLogicalToInertial<Dim>,
246 : db::ComputeTag {
247 0 : static constexpr size_t dim = Dim;
248 0 : using base = InverseHessianLogicalToInertial<Dim>;
249 0 : using return_type = typename base::type;
250 0 : using argument_tags =
251 : tmpl::list<MapTagGridToInertial,
252 : evolution::dg::subcell::Tags::Coordinates<dim, Frame::Grid>,
253 : ::Tags::Time, ::domain::Tags::FunctionsOfTime,
254 : InverseJacobianLogicalToGrid<Dim>,
255 : InverseHessianLogicalToGrid<Dim>>;
256 0 : static void function(
257 : const gsl::not_null<return_type*> inverse_hessian_logical_to_inertial,
258 : const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
259 : grid_to_inertial_map,
260 : const tnsr::I<DataVector, dim, Frame::Grid>& grid_coords,
261 : const double time,
262 : const std::unordered_map<
263 : std::string,
264 : std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
265 : functions_of_time,
266 : const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
267 : Frame::Grid>& inverse_jacobian_logical_to_grid,
268 : const ::InverseHessian<DataVector, Dim, Frame::ElementLogical,
269 : Frame::Grid>& inverse_hessian_logical_to_grid) {
270 : if (grid_to_inertial_map.is_identity()) {
271 : // Optimization for time-independent maps; we just point to the
272 : // logical-to-grid inverse hessian.
273 : const size_t num_pts = inverse_hessian_logical_to_grid[0].size();
274 : for (size_t storage_index = 0;
275 : storage_index < inverse_hessian_logical_to_grid.size();
276 : ++storage_index) {
277 : make_const_view(
278 : make_not_null(&std::as_const(
279 : (*inverse_hessian_logical_to_inertial)[storage_index])),
280 : inverse_hessian_logical_to_grid[storage_index], 0, num_pts);
281 : }
282 : } else {
283 : set_number_of_grid_points(inverse_hessian_logical_to_inertial,
284 : get<0>(grid_coords).size());
285 : // Get grid to inertial inverse jacobian
286 : const auto& inverse_jacobian_grid_to_inertial =
287 : grid_to_inertial_map.inv_jacobian(grid_coords, time,
288 : functions_of_time);
289 : // Get grid to inertial inverse hessian
290 : const auto& inverse_hessian_grid_to_inertial =
291 : grid_to_inertial_map.inv_hessian(grid_coords,
292 : inverse_jacobian_grid_to_inertial,
293 : time, functions_of_time);
294 : for (size_t i = 0; i < Dim; i++) {
295 : for (size_t j = 0; j < Dim; j++) {
296 : for (size_t k = j; k < Dim; k++) {
297 : auto& inv_hessian_component =
298 : inverse_hessian_logical_to_inertial->get(i, j, k);
299 : inv_hessian_component = 0.;
300 : for (size_t n = 0; n < Dim; n++) {
301 : inv_hessian_component +=
302 : inverse_hessian_grid_to_inertial.get(n, j, k) *
303 : inverse_jacobian_logical_to_grid.get(i, n);
304 : for (size_t m = 0; m < Dim; m++) {
305 : inv_hessian_component +=
306 : inverse_jacobian_grid_to_inertial.get(n, k) *
307 : inverse_hessian_logical_to_grid.get(i, m, n) *
308 : inverse_jacobian_grid_to_inertial.get(m, j);
309 : }
310 : }
311 : }
312 : }
313 : }
314 : }
315 : }
316 : };
317 : #endif // SPECTRE_AUTODIFF
318 : } // namespace evolution::dg::subcell::fd::Tags
|
