CubicScale.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <array>
7 #include <boost/optional.hpp>
8 #include <cstddef>
9 #include <limits>
10 #include <memory>
11 #include <string>
12 #include <unordered_map>
13
15 #include "Utilities/TypeTraits/RemoveReferenceWrapper.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 class CubicScale {
100  public:
101  static constexpr size_t dim = Dim;
102
103  explicit CubicScale(double outer_boundary,
104  std::string function_of_time_name_a,
105  std::string function_of_time_name_b) noexcept;
106  CubicScale() = default;
107
108  template <typename T>
109  std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
110  const std::array<T, Dim>& source_coords, double time,
111  const std::unordered_map<
112  std::string,
114  functions_of_time) const noexcept;
115
116  /// Returns boost::none if the point is outside the range of the map.
117  template <typename T>
118  boost::optional<std::array<tt::remove_cvref_wrap_t<T>, Dim>> inverse(
119  const std::array<T, Dim>& target_coords, double time,
120  const std::unordered_map<
121  std::string,
123  functions_of_time) const noexcept;
124
125  template <typename T>
126  std::array<tt::remove_cvref_wrap_t<T>, Dim> frame_velocity(
127  const std::array<T, Dim>& source_coords, double time,
128  const std::unordered_map<
129  std::string,
131  functions_of_time) const noexcept;
132
133  template <typename T>
134  tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
135  const std::array<T, Dim>& source_coords, double time,
136  const std::unordered_map<
137  std::string,
139  functions_of_time) const noexcept;
140
141  template <typename T>
142  tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
143  const std::array<T, Dim>& source_coords, double time,
144  const std::unordered_map<
145  std::string,
147  functions_of_time) const noexcept;
148
150  void pup(PUP::er& p) noexcept; // NOLINT
151
152  static bool is_identity() noexcept { return false; }
153
154  private:
155  template <size_t LocalDim>
157  friend bool operator==(const CubicScale<LocalDim>& lhs,
158  const CubicScale<LocalDim>& rhs) noexcept;
159
160  std::string f_of_t_a_{};
161  std::string f_of_t_b_{};
162  double one_over_outer_boundary_{std::numeric_limits<double>::signaling_NaN()};
163  bool functions_of_time_equal_{false};
164 };
165
166 template <size_t Dim>
167 bool operator!=(const CubicScale<Dim>& lhs,
168  const CubicScale<Dim>& rhs) noexcept {
169  return not(lhs == rhs);
170 }
171
172 } // namespace TimeDependent
173 } // namespace CoordinateMaps
174 } // namespace domain
std::string
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
cstddef
domain::CoordinateMaps::TimeDependent::CubicScale
Maps the radius as where is the radius of the source coordinates.
Definition: CubicScale.hpp:99
array
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
memory
limits
TypeAliases.hpp
std::unique_ptr
domain::CoordinateMaps::TimeDependent::CubicScale::inverse
boost::optional< std::array< tt::remove_cvref_wrap_t< T >, Dim > > inverse(const std::array< T, Dim > &target_coords, double time, const std::unordered_map< std::string, std::unique_ptr< domain::FunctionsOfTime::FunctionOfTime >> &functions_of_time) const noexcept
Returns boost::none if the point is outside the range of the map.
Definition: CubicScale.cpp:101
unordered_map
string