Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines class CoordinateMap
6 :
7 : #pragma once
8 :
9 : #include <cstddef>
10 : #include <limits>
11 : #include <memory>
12 : #include <optional>
13 : #include <pup.h>
14 : #include <string>
15 : #include <tuple>
16 : #include <type_traits>
17 : #include <typeinfo>
18 : #include <unordered_map>
19 : #include <unordered_set>
20 : #include <utility>
21 : #include <vector>
22 :
23 : #include "DataStructures/Tensor/Tensor.hpp"
24 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
25 : #include "Utilities/Autodiff/Autodiff.hpp"
26 : #include "Utilities/PrettyType.hpp"
27 : #include "Utilities/Serialization/CharmPupable.hpp"
28 : #include "Utilities/Simd/Simd.hpp"
29 : #include "Utilities/TMPL.hpp"
30 : #include "Utilities/TypeTraits/CreateIsCallable.hpp"
31 :
32 : /// \cond
33 : class DataVector;
34 : /// \endcond
35 :
36 : namespace domain {
37 : /// Contains all coordinate maps.
38 : namespace CoordinateMaps {
39 : /// Contains the time-dependent coordinate maps
40 1 : namespace TimeDependent {}
41 : template <typename FirstMap, typename... Maps>
42 0 : constexpr size_t map_dim = FirstMap::dim;
43 : } // namespace CoordinateMaps
44 :
45 : namespace CoordinateMap_detail {
46 : CREATE_IS_CALLABLE(function_of_time_names)
47 : CREATE_IS_CALLABLE_V(function_of_time_names)
48 :
49 : template <typename T>
50 : struct map_type {
51 : using type = T;
52 : };
53 :
54 : template <typename T>
55 : struct map_type<std::unique_ptr<T>> {
56 : using type = T;
57 : };
58 :
59 : template <typename T>
60 : using map_type_t = typename map_type<T>::type;
61 :
62 : template <typename... Maps, size_t... Is>
63 : std::unordered_set<std::string> initialize_names(
64 : const std::tuple<Maps...>& maps, std::index_sequence<Is...> /*meta*/) {
65 : std::unordered_set<std::string> function_of_time_names{};
66 : const auto add_names = [&function_of_time_names, &maps](auto index) {
67 : const auto& map = std::get<decltype(index)::value>(maps);
68 : using TupleMap = std::decay_t<decltype(map)>;
69 : constexpr bool map_is_unique_ptr = tt::is_a_v<std::unique_ptr, TupleMap>;
70 : using Map = map_type_t<TupleMap>;
71 : if constexpr (is_function_of_time_names_callable_v<Map>) {
72 : if constexpr (map_is_unique_ptr) {
73 : const auto& names = map->function_of_time_names();
74 : function_of_time_names.insert(names.begin(), names.end());
75 : } else {
76 : const auto& names = map.function_of_time_names();
77 : function_of_time_names.insert(names.begin(), names.end());
78 : }
79 : } else {
80 : (void)function_of_time_names;
81 : }
82 : };
83 : EXPAND_PACK_LEFT_TO_RIGHT(add_names(std::integral_constant<size_t, Is>{}));
84 :
85 : return function_of_time_names;
86 : }
87 :
88 : template <typename... Maps>
89 : std::string get_unsupported_autodiff_maps_error() {
90 : std::string unsupported_maps;
91 : size_t index = 0;
92 : (
93 : [&]() {
94 : if constexpr (not Maps::supports_hessian) {
95 : if (not unsupported_maps.empty()) {
96 : unsupported_maps += ", ";
97 : }
98 : unsupported_maps += "Map " + std::to_string(index) + " (" +
99 : pretty_type::get_name<Maps>() + ")";
100 : }
101 : ++index;
102 : }(),
103 : ...);
104 : return unsupported_maps;
105 : }
106 :
107 : } // namespace CoordinateMap_detail
108 :
109 : /*!
110 : * \ingroup CoordinateMapsGroup
111 : * \brief Abstract base class for CoordinateMap
112 : */
113 : template <typename SourceFrame, typename TargetFrame, size_t Dim>
114 1 : class CoordinateMapBase : public PUP::able {
115 : public:
116 0 : static constexpr size_t dim = Dim;
117 0 : using source_frame = SourceFrame;
118 0 : using target_frame = TargetFrame;
119 :
120 0 : WRAPPED_PUPable_abstract(CoordinateMapBase); // NOLINT
121 :
122 0 : CoordinateMapBase() = default;
123 0 : CoordinateMapBase(const CoordinateMapBase& /*rhs*/) = default;
124 0 : CoordinateMapBase& operator=(const CoordinateMapBase& /*rhs*/) = default;
125 0 : CoordinateMapBase(CoordinateMapBase&& /*rhs*/) = default;
126 0 : CoordinateMapBase& operator=(CoordinateMapBase&& /*rhs*/) = default;
127 0 : ~CoordinateMapBase() override = default;
128 :
129 : virtual std::unique_ptr<CoordinateMapBase<SourceFrame, TargetFrame, Dim>>
130 0 : get_clone() const = 0;
131 :
132 : /// \brief Retrieve the same map but going from `SourceFrame` to
133 : /// `Frame::Grid`.
134 : ///
135 : /// This functionality is needed when composing time-dependent maps with
136 : /// time-independent maps, where the target frame of the time-independent map
137 : /// is `Frame::Grid`.
138 : virtual std::unique_ptr<CoordinateMapBase<SourceFrame, Frame::Grid, Dim>>
139 1 : get_to_grid_frame() const = 0;
140 :
141 : /// Returns `true` if the map is the identity
142 1 : virtual bool is_identity() const = 0;
143 :
144 : /// Returns `true` if the inverse Jacobian depends on time.
145 1 : virtual bool inv_jacobian_is_time_dependent() const = 0;
146 :
147 : /// Returns `true` if the Jacobian depends on time.
148 1 : virtual bool jacobian_is_time_dependent() const = 0;
149 :
150 : /// Get a set of all FunctionOfTime names used in this mapping
151 1 : virtual const std::unordered_set<std::string>& function_of_time_names()
152 : const = 0;
153 :
154 : /// Returns `true` if this CoordinateMap supports autodiff
155 1 : virtual bool supports_hessian() const = 0;
156 :
157 : /// @{
158 : /// Apply the `Maps` to the point(s) `source_point`
159 1 : virtual tnsr::I<double, Dim, TargetFrame> operator()(
160 : tnsr::I<double, Dim, SourceFrame> source_point,
161 : double time = std::numeric_limits<double>::signaling_NaN(),
162 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
163 1 : virtual tnsr::I<DataVector, Dim, TargetFrame> operator()(
164 : tnsr::I<DataVector, Dim, SourceFrame> source_point,
165 : double time = std::numeric_limits<double>::signaling_NaN(),
166 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
167 : /// @}
168 :
169 : #ifdef SPECTRE_AUTODIFF
170 : /// @{
171 : /// Apply the `Maps` to the point(s) `source_point`
172 : ///
173 : /// \note Require SPECTRE_AUTODIFF=ON.
174 : virtual tnsr::I<autodiff::HigherOrderDual<2, double>, Dim, TargetFrame>
175 1 : operator()(tnsr::I<autodiff::HigherOrderDual<2, double>, Dim,
176 : SourceFrame> /*source_point*/,
177 : double /*time*/ = std::numeric_limits<double>::signaling_NaN(),
178 : const FunctionsOfTimeMap& /*functions_of_time*/ = {}) const {
179 : ERROR("Call operator for autodiff types must be overriden.");
180 : }
181 : virtual tnsr::I<autodiff::HigherOrderDual<2, simd::batch<double>>, Dim,
182 : TargetFrame>
183 1 : operator()(tnsr::I<autodiff::HigherOrderDual<2, simd::batch<double>>, Dim,
184 : SourceFrame> /*source_point*/,
185 : double /*time*/ = std::numeric_limits<double>::signaling_NaN(),
186 : const FunctionsOfTimeMap& /*functions_of_time*/ = {}) const {
187 : ERROR("Call operator for autodiff types must be overriden.");
188 : }
189 : /// @}
190 : #endif // SPECTRE_AUTODIFF
191 :
192 : /// @{
193 : /// Apply the inverse `Maps` to the point(s) `target_point`.
194 : /// The returned std::optional is invalid if the map is not invertible
195 : /// at `target_point`, or if `target_point` can be easily determined to not
196 : /// make sense for the map. An example of the latter is passing a
197 : /// point with a negative value of z into a positive-z Wedge<3> inverse map.
198 : /// The inverse function is only callable with doubles because the inverse
199 : /// might fail if called for a point out of range, and it is unclear
200 : /// what should happen if the inverse were to succeed for some points in a
201 : /// DataVector but fail for other points.
202 1 : virtual std::optional<tnsr::I<double, Dim, SourceFrame>> inverse(
203 : tnsr::I<double, Dim, TargetFrame> target_point,
204 : double time = std::numeric_limits<double>::signaling_NaN(),
205 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
206 : /// @}
207 :
208 : /// @{
209 : /// Compute the inverse Jacobian of the `Maps` at the point(s)
210 : /// `source_point`
211 1 : virtual InverseJacobian<double, Dim, SourceFrame, TargetFrame> inv_jacobian(
212 : tnsr::I<double, Dim, SourceFrame> source_point,
213 : double time = std::numeric_limits<double>::signaling_NaN(),
214 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
215 : virtual InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>
216 1 : inv_jacobian(tnsr::I<DataVector, Dim, SourceFrame> source_point,
217 : double time = std::numeric_limits<double>::signaling_NaN(),
218 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
219 : #ifdef SPECTRE_AUTODIFF
220 : virtual InverseJacobian<autodiff::HigherOrderDual<2, double>, Dim,
221 : SourceFrame, TargetFrame>
222 1 : inv_jacobian(tnsr::I<autodiff::HigherOrderDual<2, double>, Dim, SourceFrame>
223 : source_point,
224 : double time = std::numeric_limits<double>::signaling_NaN(),
225 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
226 : #endif // SPECTRE_AUTODIFF
227 : /// @}
228 :
229 : #ifdef SPECTRE_AUTODIFF
230 : /// @{
231 : /// Compute the inverse Hessian of the `Maps` at the point(s)
232 : /// `source_point`
233 : ///
234 : /// \note Require SPECTRE_AUTODIFF=ON.
235 1 : virtual InverseHessian<double, Dim, SourceFrame, TargetFrame> inv_hessian(
236 : tnsr::I<double, Dim, SourceFrame> source_point,
237 : const InverseJacobian<double, Dim, SourceFrame, TargetFrame>& inverse_jac,
238 : double time = std::numeric_limits<double>::signaling_NaN(),
239 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
240 1 : virtual InverseHessian<DataVector, Dim, SourceFrame, TargetFrame> inv_hessian(
241 : tnsr::I<DataVector, Dim, SourceFrame> source_point,
242 : const InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>&
243 : inverse_jac,
244 : double time = std::numeric_limits<double>::signaling_NaN(),
245 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
246 : /// @}
247 : #endif // SPECTRE_AUTODIFF
248 :
249 : /// @{
250 : /// Compute the Jacobian of the `Maps` at the point(s) `source_point`
251 1 : virtual Jacobian<double, Dim, SourceFrame, TargetFrame> jacobian(
252 : tnsr::I<double, Dim, SourceFrame> source_point,
253 : double time = std::numeric_limits<double>::signaling_NaN(),
254 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
255 1 : virtual Jacobian<DataVector, Dim, SourceFrame, TargetFrame> jacobian(
256 : tnsr::I<DataVector, Dim, SourceFrame> source_point,
257 : double time = std::numeric_limits<double>::signaling_NaN(),
258 : const FunctionsOfTimeMap& functions_of_time = {}) const = 0;
259 : /// @}
260 :
261 : /// @{
262 : /// Compute the mapped coordinates, frame velocity, Jacobian, and inverse
263 : /// Jacobian
264 1 : virtual auto coords_frame_velocity_jacobians(
265 : tnsr::I<double, Dim, SourceFrame> source_point,
266 : double time = std::numeric_limits<double>::signaling_NaN(),
267 : const FunctionsOfTimeMap& functions_of_time = {}) const
268 : -> std::tuple<tnsr::I<double, Dim, TargetFrame>,
269 : InverseJacobian<double, Dim, SourceFrame, TargetFrame>,
270 : Jacobian<double, Dim, SourceFrame, TargetFrame>,
271 : tnsr::I<double, Dim, TargetFrame>> = 0;
272 1 : virtual auto coords_frame_velocity_jacobians(
273 : tnsr::I<DataVector, Dim, SourceFrame> source_point,
274 : double time = std::numeric_limits<double>::signaling_NaN(),
275 : const FunctionsOfTimeMap& functions_of_time = {}) const
276 : -> std::tuple<tnsr::I<DataVector, Dim, TargetFrame>,
277 : InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>,
278 : Jacobian<DataVector, Dim, SourceFrame, TargetFrame>,
279 : tnsr::I<DataVector, Dim, TargetFrame>> = 0;
280 : /// @}
281 :
282 : private:
283 0 : virtual bool is_equal_to(const CoordinateMapBase& other) const = 0;
284 0 : friend bool operator==(const CoordinateMapBase& lhs,
285 : const CoordinateMapBase& rhs) {
286 : return typeid(lhs) == typeid(rhs) and lhs.is_equal_to(rhs);
287 : }
288 0 : friend bool operator!=(const CoordinateMapBase& lhs,
289 : const CoordinateMapBase& rhs) {
290 : return not(lhs == rhs);
291 : }
292 : };
293 :
294 : /*!
295 : * \ingroup CoordinateMapsGroup
296 : * \brief A coordinate map or composition of coordinate maps
297 : *
298 : * Maps coordinates from the `SourceFrame` to the `TargetFrame` using the
299 : * coordinate maps `Maps...`. The individual maps are applied left to right
300 : * from the source to the target Frame. The inverse map, as well as Jacobian
301 : * and inverse Jacobian are also provided. The `CoordinateMap` class must
302 : * be used even if just wrapping a single coordinate map. It is designed to
303 : * be an extremely minimal interface to the underlying coordinate maps. For
304 : * a list of all coordinate maps see the CoordinateMaps group or namespace.
305 : *
306 : * Each coordinate map must contain a `static constexpr size_t dim` variable
307 : * that is equal to the dimensionality of the map. The Coordinatemap class
308 : * contains a member `static constexpr size_t dim`, a type alias `source_frame`,
309 : * a type alias `target_frame` and typelist of the `Maps...`.
310 : */
311 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
312 1 : class CoordinateMap
313 : : public CoordinateMapBase<SourceFrame, TargetFrame,
314 : CoordinateMaps::map_dim<Maps...>> {
315 : static_assert(sizeof...(Maps) > 0, "Must have at least one map");
316 : static_assert(
317 : tmpl::all<tmpl::integral_list<size_t, Maps::dim...>,
318 : std::is_same<tmpl::integral_constant<
319 : size_t, CoordinateMaps::map_dim<Maps...>>,
320 : tmpl::_1>>::value,
321 : "All Maps passed to CoordinateMap must be of the same dimensionality.");
322 :
323 : public:
324 0 : static constexpr size_t dim = CoordinateMaps::map_dim<Maps...>;
325 0 : using source_frame = SourceFrame;
326 0 : using target_frame = TargetFrame;
327 0 : using maps_list = tmpl::list<Maps...>;
328 :
329 : /// Used for Charm++ serialization
330 1 : CoordinateMap() = default;
331 :
332 0 : CoordinateMap(const CoordinateMap& /*rhs*/) = default;
333 0 : CoordinateMap& operator=(const CoordinateMap& /*rhs*/) = default;
334 0 : CoordinateMap(CoordinateMap&& /*rhs*/) = default;
335 0 : CoordinateMap& operator=(CoordinateMap&& /*rhs*/) = default;
336 0 : ~CoordinateMap() override = default;
337 :
338 0 : explicit CoordinateMap(Maps... maps);
339 :
340 0 : std::unique_ptr<CoordinateMapBase<SourceFrame, TargetFrame, dim>> get_clone()
341 : const override {
342 : return std::make_unique<CoordinateMap>(*this);
343 : }
344 :
345 : std::unique_ptr<CoordinateMapBase<SourceFrame, Frame::Grid, dim>>
346 1 : get_to_grid_frame() const override {
347 : return get_to_grid_frame_impl(std::make_index_sequence<sizeof...(Maps)>{});
348 : }
349 :
350 : /// Returns `true` if the map is the identity
351 1 : bool is_identity() const override;
352 :
353 : /// Returns `true` if the inverse Jacobian depends on time.
354 1 : bool inv_jacobian_is_time_dependent() const override;
355 :
356 : /// Returns `true` if the Jacobian depends on time.
357 1 : bool jacobian_is_time_dependent() const override;
358 :
359 : /// Get a set of all FunctionOfTime names from `Maps`
360 1 : const std::unordered_set<std::string>& function_of_time_names()
361 : const override {
362 : return function_of_time_names_;
363 : }
364 :
365 : /// Returns `true` if this coordinate map supports hessian
366 1 : bool supports_hessian() const override {
367 : return (Maps::supports_hessian && ...);
368 : }
369 :
370 : /// @{
371 : /// Apply the `Maps...` to the point(s) `source_point`
372 1 : tnsr::I<double, dim, TargetFrame> operator()(
373 : tnsr::I<double, dim, SourceFrame> source_point,
374 : const double time = std::numeric_limits<double>::signaling_NaN(),
375 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
376 : return call_impl(std::move(source_point), time, functions_of_time,
377 : std::make_index_sequence<sizeof...(Maps)>{});
378 : }
379 1 : tnsr::I<DataVector, dim, TargetFrame> operator()(
380 : tnsr::I<DataVector, dim, SourceFrame> source_point,
381 : const double time = std::numeric_limits<double>::signaling_NaN(),
382 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
383 : return call_impl(std::move(source_point), time, functions_of_time,
384 : std::make_index_sequence<sizeof...(Maps)>{});
385 : }
386 : /// @}
387 :
388 : #ifdef SPECTRE_AUTODIFF
389 : /// @{
390 : /// Apply the `Maps...` to the point(s) `source_point`
391 : ///
392 : /// \note Require SPECTRE_AUTODIFF=ON.
393 1 : tnsr::I<autodiff::HigherOrderDual<2, double>, dim, TargetFrame> operator()(
394 : tnsr::I<autodiff::HigherOrderDual<2, double>, dim, SourceFrame>
395 : source_point,
396 : const double time = std::numeric_limits<double>::signaling_NaN(),
397 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
398 : if constexpr ((Maps::supports_hessian && ...)) {
399 : return call_impl(std::move(source_point), time, functions_of_time,
400 : std::make_index_sequence<sizeof...(Maps)>{});
401 : } else {
402 : ERROR("At least one of the Maps does not support autodiff: "
403 : << CoordinateMap_detail::get_unsupported_autodiff_maps_error<
404 : Maps...>());
405 : }
406 : }
407 : tnsr::I<autodiff::HigherOrderDual<2, simd::batch<double>>, dim, TargetFrame>
408 1 : operator()(tnsr::I<autodiff::HigherOrderDual<2, simd::batch<double>>, dim,
409 : SourceFrame>
410 : source_point,
411 : const double time = std::numeric_limits<double>::signaling_NaN(),
412 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
413 : if constexpr ((Maps::supports_hessian && ...)) {
414 : return call_impl(std::move(source_point), time, functions_of_time,
415 : std::make_index_sequence<sizeof...(Maps)>{});
416 : } else {
417 : ERROR("At least one of the Maps does not support autodiff: "
418 : << CoordinateMap_detail::get_unsupported_autodiff_maps_error<
419 : Maps...>());
420 : }
421 : }
422 : /// @}
423 : #endif // SPECTRE_AUTODIFF
424 :
425 : /// @{
426 : /// Apply the inverse `Maps...` to the point(s) `target_point`
427 1 : std::optional<tnsr::I<double, dim, SourceFrame>> inverse(
428 : tnsr::I<double, dim, TargetFrame> target_point,
429 : const double time = std::numeric_limits<double>::signaling_NaN(),
430 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
431 : return inverse_impl(std::move(target_point), time, functions_of_time,
432 : std::make_index_sequence<sizeof...(Maps)>{});
433 : }
434 : /// @}
435 :
436 : /// @{
437 : /// Compute the inverse Jacobian of the `Maps...` at the point(s)
438 : /// `source_point`
439 1 : InverseJacobian<double, dim, SourceFrame, TargetFrame> inv_jacobian(
440 : tnsr::I<double, dim, SourceFrame> source_point,
441 : const double time = std::numeric_limits<double>::signaling_NaN(),
442 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
443 : return inv_jacobian_impl(std::move(source_point), time, functions_of_time);
444 : }
445 1 : InverseJacobian<DataVector, dim, SourceFrame, TargetFrame> inv_jacobian(
446 : tnsr::I<DataVector, dim, SourceFrame> source_point,
447 : const double time = std::numeric_limits<double>::signaling_NaN(),
448 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
449 : return inv_jacobian_impl(std::move(source_point), time, functions_of_time);
450 : }
451 : #ifdef SPECTRE_AUTODIFF
452 : InverseJacobian<autodiff::HigherOrderDual<2, double>, dim, SourceFrame,
453 : TargetFrame>
454 1 : inv_jacobian(
455 : tnsr::I<autodiff::HigherOrderDual<2, double>, dim, SourceFrame>
456 : source_point,
457 : const double time = std::numeric_limits<double>::signaling_NaN(),
458 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
459 : if constexpr ((Maps::supports_hessian && ...)) {
460 : return inv_jacobian_impl(std::move(source_point), time,
461 : functions_of_time);
462 : } else {
463 : ERROR("At least one of the Maps does not support autodiff: "
464 : << CoordinateMap_detail::get_unsupported_autodiff_maps_error<
465 : Maps...>());
466 : }
467 : }
468 : #endif // SPECTRE_AUTODIFF
469 : /// @}
470 :
471 : #ifdef SPECTRE_AUTODIFF
472 : /// @{
473 : /*!
474 : * Compute the inverse Hessian of the `Maps...` at the point(s)
475 : * `source_point` by computing the Hessian of the `Maps...` and
476 : * compose it with inverse Jacobians of the `Maps...`.
477 : *
478 : * This function propagates autodiff dual types
479 : * through the `call_impl` function to automatically get the
480 : * Hessian \f$ \frac{\partial^2 x^i}{\partial\xi^j \partial\xi^k} \f$
481 : * where \f$ x^i \f$ are the target coordinates and \f$ \xi^i \f$ are the
482 : * source coordinates. Then the inverse Hessian is computed by the
483 : * identity
484 : * \f[\frac{\partial^2 \xi^i}{\partial x^m \partial x^n} =
485 : * -\frac{\partial \xi^i}{\partial x^j} \frac{\partial \xi^k}{\partial x^m}
486 : * \frac{\partial \xi^l}{\partial x^n} \frac{\partial^2 x^j}{\partial \xi^k
487 : * \partial \xi^l}, \f] where the inverse Jacobian is passed in as a function
488 : * argument.
489 : *
490 : * See Test_CoordinateMap.cpp for a different implementation. The current
491 : * implementation is chosen in production as it is faster when we have
492 : * the inverse Jacobian already, and in most cases we do.
493 : *
494 : * We use forward mode autodiff here as it is simpler to implement and
495 : * has better optimization than the reverse mode. Reverse mode in
496 : * the [Autodiff](https://github.com/autodiff/autodiff) library
497 : * has higher cost per propagation and does not support taping.
498 : * Also see this
499 : * [github issue](https://github.com/autodiff/autodiff/issues/332).
500 : *
501 : * \note Require SPECTRE_AUTODIFF=ON.
502 : */
503 1 : InverseHessian<double, dim, SourceFrame, TargetFrame> inv_hessian(
504 : tnsr::I<double, dim, SourceFrame> source_point,
505 : const ::InverseJacobian<double, dim, SourceFrame, TargetFrame>&
506 : inverse_jac,
507 : const double time = std::numeric_limits<double>::signaling_NaN(),
508 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
509 : if constexpr ((Maps::supports_hessian && ...)) {
510 : return inv_hessian_impl(std::move(source_point), inverse_jac, time,
511 : functions_of_time);
512 : } else {
513 : ERROR("At least one of the Maps does not support autodiff: "
514 : << CoordinateMap_detail::get_unsupported_autodiff_maps_error<
515 : Maps...>());
516 : }
517 : }
518 1 : InverseHessian<DataVector, dim, SourceFrame, TargetFrame> inv_hessian(
519 : tnsr::I<DataVector, dim, SourceFrame> source_point,
520 : const ::InverseJacobian<DataVector, dim, SourceFrame, TargetFrame>&
521 : inverse_jac,
522 : const double time = std::numeric_limits<double>::signaling_NaN(),
523 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
524 : if constexpr ((Maps::supports_hessian && ...)) {
525 : return inv_hessian_impl(std::move(source_point), inverse_jac, time,
526 : functions_of_time);
527 : } else {
528 : ERROR("At least one of the Maps does not support autodiff: "
529 : << CoordinateMap_detail::get_unsupported_autodiff_maps_error<
530 : Maps...>());
531 : }
532 : }
533 : /// @}
534 : #endif // SPECTRE_AUTODIFF
535 :
536 : /// @{
537 : /// Compute the Jacobian of the `Maps...` at the point(s) `source_point`
538 1 : Jacobian<double, dim, SourceFrame, TargetFrame> jacobian(
539 : tnsr::I<double, dim, SourceFrame> source_point,
540 : const double time = std::numeric_limits<double>::signaling_NaN(),
541 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
542 : return jacobian_impl(std::move(source_point), time, functions_of_time);
543 : }
544 1 : Jacobian<DataVector, dim, SourceFrame, TargetFrame> jacobian(
545 : tnsr::I<DataVector, dim, SourceFrame> source_point,
546 : const double time = std::numeric_limits<double>::signaling_NaN(),
547 : const FunctionsOfTimeMap& functions_of_time = {}) const override {
548 : return jacobian_impl(std::move(source_point), time, functions_of_time);
549 : }
550 : /// @}
551 :
552 : /// @{
553 : /// Compute the mapped coordinates, frame velocity, Jacobian, and inverse
554 : /// Jacobian. The inverse Jacobian is computed by numerically inverting the
555 : /// Jacobian as this was measured to be quicker than computing it directly for
556 : /// more complex map concatenations.
557 1 : auto coords_frame_velocity_jacobians(
558 : tnsr::I<double, dim, SourceFrame> source_point,
559 : const double time = std::numeric_limits<double>::signaling_NaN(),
560 : const FunctionsOfTimeMap& functions_of_time = {}) const
561 : -> std::tuple<tnsr::I<double, dim, TargetFrame>,
562 : InverseJacobian<double, dim, SourceFrame, TargetFrame>,
563 : Jacobian<double, dim, SourceFrame, TargetFrame>,
564 : tnsr::I<double, dim, TargetFrame>> override {
565 : return coords_frame_velocity_jacobians_impl(std::move(source_point), time,
566 : functions_of_time);
567 : }
568 1 : auto coords_frame_velocity_jacobians(
569 : tnsr::I<DataVector, dim, SourceFrame> source_point,
570 : const double time = std::numeric_limits<double>::signaling_NaN(),
571 : const FunctionsOfTimeMap& functions_of_time = {}) const
572 : -> std::tuple<tnsr::I<DataVector, dim, TargetFrame>,
573 : InverseJacobian<DataVector, dim, SourceFrame, TargetFrame>,
574 : Jacobian<DataVector, dim, SourceFrame, TargetFrame>,
575 : tnsr::I<DataVector, dim, TargetFrame>> override {
576 : return coords_frame_velocity_jacobians_impl(std::move(source_point), time,
577 : functions_of_time);
578 : }
579 : /// @}
580 :
581 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
582 : SINGLE_ARG(CoordinateMapBase<SourceFrame, TargetFrame, dim>),
583 : CoordinateMap);
584 :
585 0 : explicit CoordinateMap(CkMigrateMessage* /*unused*/) {}
586 :
587 : // NOLINTNEXTLINE(google-runtime-references)
588 0 : void pup(PUP::er& p) override {
589 : size_t version = 0;
590 : p | version;
591 : // Remember to increment the version number when making changes to this
592 : // function. Retain support for unpacking data written by previous versions
593 : // whenever possible. See `Domain` docs for details.
594 : if (version >= 0) {
595 : CoordinateMapBase<SourceFrame, TargetFrame, dim>::pup(p);
596 : p | maps_;
597 : }
598 :
599 : // No need to pup this because it is uniquely determined by the maps
600 : if (p.isUnpacking()) {
601 : function_of_time_names_ = CoordinateMap_detail::initialize_names(
602 : maps_, std::make_index_sequence<sizeof...(Maps)>{});
603 : }
604 : }
605 :
606 : private:
607 0 : friend bool operator==(const CoordinateMap& lhs, const CoordinateMap& rhs) {
608 : return lhs.maps_ == rhs.maps_;
609 : }
610 :
611 : template <typename NewMap, typename LocalSourceFrame,
612 : typename LocalTargetFrame, typename... LocalMaps, size_t... Is>
613 : friend CoordinateMap<LocalSourceFrame, LocalTargetFrame, LocalMaps..., NewMap>
614 : // NOLINTNEXTLINE(readability-redundant-declaration,-warnings-as-errors)
615 0 : push_back_impl(
616 : CoordinateMap<LocalSourceFrame, LocalTargetFrame, LocalMaps...>&& old_map,
617 : NewMap new_map, std::index_sequence<Is...> /*meta*/);
618 :
619 : template <typename... NewMaps, typename LocalSourceFrame,
620 : typename LocalTargetFrame, typename... LocalMaps, size_t... Is,
621 : size_t... Js>
622 : friend CoordinateMap<LocalSourceFrame, LocalTargetFrame, LocalMaps...,
623 : NewMaps...>
624 : // NOLINTNEXTLINE(readability-redundant-declaration,-warnings-as-errors)
625 0 : push_back_impl(
626 : CoordinateMap<LocalSourceFrame, LocalTargetFrame, LocalMaps...>&& old_map,
627 : CoordinateMap<LocalSourceFrame, LocalTargetFrame, NewMaps...> new_map,
628 : std::index_sequence<Is...> /*meta*/, std::index_sequence<Js...> /*meta*/);
629 :
630 : template <typename NewMap, typename LocalSourceFrame,
631 : typename LocalTargetFrame, typename... LocalMaps, size_t... Is>
632 : friend CoordinateMap<LocalSourceFrame, LocalTargetFrame, NewMap, LocalMaps...>
633 : // NOLINTNEXTLINE(readability-redundant-declaration,-warnings-as-errors)
634 0 : push_front_impl(
635 : CoordinateMap<LocalSourceFrame, LocalTargetFrame, LocalMaps...>&& old_map,
636 : NewMap new_map, std::index_sequence<Is...> /*meta*/);
637 :
638 : template <size_t... Is>
639 : std::unique_ptr<CoordinateMapBase<SourceFrame, Frame::Grid, dim>>
640 0 : get_to_grid_frame_impl(std::index_sequence<Is...> /*meta*/) const;
641 :
642 0 : bool is_equal_to(const CoordinateMapBase<SourceFrame, TargetFrame, dim>&
643 : other) const override {
644 : const auto& cast_of_other = dynamic_cast<const CoordinateMap&>(other);
645 : return *this == cast_of_other;
646 : }
647 :
648 0 : void check_functions_of_time(
649 : const FunctionsOfTimeMap& functions_of_time) const;
650 :
651 : template <typename T, size_t... Is>
652 0 : tnsr::I<T, dim, TargetFrame> call_impl(
653 : tnsr::I<T, dim, SourceFrame>&& source_point, double time,
654 : const FunctionsOfTimeMap& functions_of_time,
655 : std::index_sequence<Is...> /*meta*/) const;
656 :
657 : template <typename T, size_t... Is>
658 0 : std::optional<tnsr::I<T, dim, SourceFrame>> inverse_impl(
659 : tnsr::I<T, dim, TargetFrame>&& target_point, double time,
660 : const FunctionsOfTimeMap& functions_of_time,
661 : std::index_sequence<Is...> /*meta*/) const;
662 :
663 : template <typename T>
664 0 : InverseJacobian<T, dim, SourceFrame, TargetFrame> inv_jacobian_impl(
665 : tnsr::I<T, dim, SourceFrame>&& source_point, double time,
666 : const FunctionsOfTimeMap& functions_of_time) const;
667 :
668 : template <typename T>
669 0 : InverseHessian<T, dim, SourceFrame, TargetFrame> inv_hessian_impl(
670 : tnsr::I<T, dim, SourceFrame>&& source_point,
671 : const ::InverseJacobian<T, dim, SourceFrame, TargetFrame>& inverse_jac,
672 : double time, const FunctionsOfTimeMap& functions_of_time) const;
673 :
674 : template <typename T>
675 0 : Jacobian<T, dim, SourceFrame, TargetFrame> jacobian_impl(
676 : tnsr::I<T, dim, SourceFrame>&& source_point, double time,
677 : const FunctionsOfTimeMap& functions_of_time) const;
678 :
679 : template <typename T>
680 : std::tuple<tnsr::I<T, dim, TargetFrame>,
681 : InverseJacobian<T, dim, SourceFrame, TargetFrame>,
682 : Jacobian<T, dim, SourceFrame, TargetFrame>,
683 : tnsr::I<T, dim, TargetFrame>>
684 0 : coords_frame_velocity_jacobians_impl(
685 : tnsr::I<T, dim, SourceFrame> source_point, double time,
686 : const FunctionsOfTimeMap& functions_of_time) const;
687 :
688 0 : std::tuple<Maps...> maps_;
689 0 : std::unordered_set<std::string> function_of_time_names_;
690 : };
691 :
692 : /// \ingroup ComputationalDomainGroup
693 : /// \brief Creates a `CoordinateMap` of `maps...`
694 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
695 1 : auto make_coordinate_map(Maps&&... maps)
696 : -> CoordinateMap<SourceFrame, TargetFrame, std::decay_t<Maps>...>;
697 :
698 : /// \ingroup ComputationalDomainGroup
699 : /// \brief Creates a `std::unique_ptr<CoordinateMapBase>` of `maps...`
700 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
701 1 : auto make_coordinate_map_base(Maps&&... maps)
702 : -> std::unique_ptr<CoordinateMapBase<
703 : SourceFrame, TargetFrame,
704 : CoordinateMap<SourceFrame, TargetFrame, std::decay_t<Maps>...>::dim>>;
705 :
706 : /// \ingroup ComputationalDomainGroup
707 : /// \brief Creates a `std::vector<std::unique_ptr<CoordinateMapBase>>`
708 : /// containing the result of `make_coordinate_map_base` applied to each
709 : /// argument passed in.
710 : template <typename SourceFrame, typename TargetFrame, typename Arg0,
711 : typename... Args>
712 1 : auto make_vector_coordinate_map_base(Arg0&& arg_0, Args&&... remaining_args)
713 : -> std::vector<std::unique_ptr<
714 : CoordinateMapBase<SourceFrame, TargetFrame, std::decay_t<Arg0>::dim>>>;
715 :
716 : /// \ingroup ComputationalDomainGroup
717 : /// \brief Creates a `std::vector<std::unique_ptr<CoordinateMapBase>>`
718 : /// containing the result of `make_coordinate_map_base` applied to each
719 : /// element of the vector of maps composed with the rest of the arguments
720 : /// passed in.
721 : template <typename SourceFrame, typename TargetFrame, size_t Dim, typename Map,
722 : typename... Maps>
723 1 : auto make_vector_coordinate_map_base(std::vector<Map> maps,
724 : const Maps&... remaining_maps)
725 : -> std::vector<
726 : std::unique_ptr<CoordinateMapBase<SourceFrame, TargetFrame, Dim>>>;
727 :
728 : /// \ingroup ComputationalDomainGroup
729 : /// \brief Creates a `CoordinateMap` by appending the new map to the end of the
730 : /// old maps
731 : template <typename SourceFrame, typename TargetFrame, typename... Maps,
732 : typename NewMap>
733 1 : CoordinateMap<SourceFrame, TargetFrame, Maps..., NewMap> push_back(
734 : CoordinateMap<SourceFrame, TargetFrame, Maps...> old_map, NewMap new_map);
735 :
736 : /// \ingroup ComputationalDomainGroup
737 : /// \brief Creates a `CoordinateMap` by prepending the new map to the beginning
738 : /// of the old maps
739 : template <typename SourceFrame, typename TargetFrame, typename... Maps,
740 : typename NewMap>
741 1 : CoordinateMap<SourceFrame, TargetFrame, NewMap, Maps...> push_front(
742 : CoordinateMap<SourceFrame, TargetFrame, Maps...> old_map, NewMap new_map);
743 :
744 : /// \cond
745 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
746 : PUP::able::PUP_ID
747 : CoordinateMap<SourceFrame, TargetFrame, Maps...>::my_PUP_ID = // NOLINT
748 : 0;
749 : /// \endcond
750 : } // namespace domain
|
