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 <memory>
10 : #include <optional>
11 : #include <string>
12 : #include <unordered_map>
13 : #include <unordered_set>
14 :
15 : #include "DataStructures/Tensor/TypeAliases.hpp"
16 :
17 : /// \cond
18 : namespace domain {
19 : namespace FunctionsOfTime {
20 : class FunctionOfTime;
21 : } // namespace FunctionsOfTime
22 : } // namespace domain
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace domain {
29 : namespace CoordinateMaps {
30 : namespace TimeDependent {
31 : /*!
32 : * \ingroup CoordMapsTimeDependentGroup
33 : * \brief Maps the radius as \f$r(t) = a(t)\rho + \left(b(t) - a(t)\right)
34 : * \frac{\rho^3} {R^2}\f$ where \f$\rho\f$ is the radius of the source
35 : * coordinates.
36 : *
37 : * The map scales the radius \f$\rho\f$ in the source coordinates
38 : * \f$\xi^{\hat{i}}\f$ by a factor \f$a(t)\f$, while the coordinates near the
39 : * outer boundary \f$R\f$, are scaled by a factor \f$b(t)\f$. Here \f$a(t)\f$
40 : * and \f$b(t)\f$ are FunctionsOfTime. The target/mapped coordinates are denoted
41 : * by \f$x^i\f$.
42 : *
43 : * The mapped coordinates are given by:
44 : *
45 : * \f{align}{
46 : * x^i = \left[a + (b-a) \frac{\rho^2}{R^2}\right] \xi^{\hat{i}}
47 : * \delta^i_{\hat{i}},
48 : * \f}
49 : *
50 : * where \f$\xi^{\hat{i}}\f$ are the source coordinates, \f$a\f$ and \f$b\f$ are
51 : * functions of time, \f$\rho\f$ is the radius in the source coordinates, and
52 : * \f$R\f$ is the outer boundary.
53 : *
54 : * The inverse map is computed by solving the cubic equation:
55 : *
56 : * \f{align}{
57 : * (b-a)\frac{\rho^3}{R^2} + a \rho - r = 0,
58 : * \f}
59 : *
60 : * which is done by defining \f$q=\rho/R\f$, and solving
61 : *
62 : * \f{align}{
63 : * q \left[(b-a) q^2 + a\right] - \frac{r}{R} = 0.
64 : * \f}
65 : *
66 : * The source coordinates are obtained using:
67 : *
68 : * \f{align}{
69 : * \xi^{\hat{i}} = \frac{qR}{r} x^i(t) \delta^{\hat{i}}_i
70 : * \f}
71 : *
72 : * The Jacobian is given by:
73 : *
74 : * \f{align}{
75 : * \frac{\partial x^i}{\partial \xi^{\hat{i}}}=
76 : * \left[a + (b-a) \frac{\rho^2}{R^2}\right] \delta^i_{\hat{i}}
77 : * + \frac{2 (b-a)}{R^2} \xi^{\hat{j}} \delta^i_{\hat{j}} \xi^{\hat{k}}
78 : * \delta_{\hat{k}\hat{i}}
79 : * \f}
80 : *
81 : * The inverse Jacobian is given by:
82 : *
83 : * \f{align}{
84 : * \frac{\partial \xi^{\hat{i}}}{\partial x^i}=
85 : * \frac{1}{\left[a + (b-a)\rho^2/R^2\right]}
86 : * \left[\delta^{\hat{i}}_i -
87 : * \frac{2 (b-a)}{\left[a R^2 + 3(b-a)\rho^2\right]}
88 : * \xi^{\hat{i}}\xi^{\hat{j}}\delta_{\hat{j}i}\right]
89 : * \f}
90 : *
91 : * The mesh velocity \f$v_g^i\f$ is given by:
92 : *
93 : * \f{align}{
94 : * v_g^i = \left[\frac{da}{dt} + \left(\frac{db}{dt}-\frac{da}{dt}\right)
95 : * \frac{\rho^2}{R^2}\right] \xi^{\hat{i}} \delta^i_{\hat{i}}.
96 : * \f}
97 : */
98 : template <size_t Dim>
99 1 : class CubicScale {
100 : public:
101 0 : static constexpr size_t dim = Dim;
102 :
103 0 : explicit CubicScale(double outer_boundary,
104 : std::string function_of_time_name_a,
105 : std::string function_of_time_name_b);
106 0 : CubicScale() = default;
107 :
108 : template <typename T>
109 0 : std::array<T, Dim> operator()(
110 : const std::array<T, Dim>& source_coords, double time,
111 : const std::unordered_map<
112 : std::string,
113 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
114 : functions_of_time) const;
115 :
116 : /// Returns std::nullopt if the point is outside the range of the map.
117 : /// The inverse function is only callable with doubles because the inverse
118 : /// might fail if called for a point out of range, and it is unclear
119 : /// what should happen if the inverse were to succeed for some points in a
120 : /// DataVector but fail for other points.
121 1 : std::optional<std::array<double, Dim>> inverse(
122 : const std::array<double, Dim>& target_coords, double time,
123 : const std::unordered_map<
124 : std::string,
125 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
126 : functions_of_time) const;
127 :
128 : template <typename T>
129 0 : std::array<T, Dim> frame_velocity(
130 : const std::array<T, Dim>& source_coords, double time,
131 : const std::unordered_map<
132 : std::string,
133 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
134 : functions_of_time) const;
135 :
136 : template <typename T>
137 0 : tnsr::Ij<T, Dim, Frame::NoFrame> inv_jacobian(
138 : const std::array<T, Dim>& source_coords, double time,
139 : const std::unordered_map<
140 : std::string,
141 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
142 : functions_of_time) const;
143 :
144 : template <typename T>
145 0 : tnsr::Ij<T, Dim, Frame::NoFrame> jacobian(
146 : const std::array<T, Dim>& source_coords, double time,
147 : const std::unordered_map<
148 : std::string,
149 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
150 : functions_of_time) const;
151 :
152 : // NOLINTNEXTLINE(google-runtime-references)
153 0 : void pup(PUP::er& p);
154 :
155 0 : static bool is_identity() { return false; }
156 :
157 0 : static constexpr bool supports_hessian{false};
158 :
159 0 : const std::unordered_set<std::string>& function_of_time_names() const {
160 : return f_of_t_names_;
161 : }
162 :
163 : private:
164 : template <size_t LocalDim>
165 : // NOLINTNEXTLINE(readability-redundant-declaration)
166 0 : friend bool operator==(const CubicScale<LocalDim>& lhs,
167 : const CubicScale<LocalDim>& rhs);
168 :
169 0 : std::string f_of_t_a_{};
170 0 : std::string f_of_t_b_{};
171 0 : std::unordered_set<std::string> f_of_t_names_;
172 0 : double one_over_outer_boundary_{std::numeric_limits<double>::signaling_NaN()};
173 0 : bool functions_of_time_equal_{false};
174 : };
175 :
176 : template <size_t Dim>
177 0 : bool operator!=(const CubicScale<Dim>& lhs, const CubicScale<Dim>& rhs) {
178 : return not(lhs == rhs);
179 : }
180 :
181 : } // namespace TimeDependent
182 : } // namespace CoordinateMaps
183 : } // namespace domain
|
