SpECTRE Documentation Coverage Report
Current view: top level - Domain/CoordinateMaps/TimeDependent - Translation.hpp Hit Total Coverage
Commit: 9a905b0737f373631c1b8e8389b8f26e67fa5313 Lines: 2 31 6.5 %
Date: 2024-03-28 09:03:18
Legend: Lines: hit not hit

          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 <memory>
       9             : #include <optional>
      10             : #include <string>
      11             : #include <unordered_map>
      12             : #include <unordered_set>
      13             : 
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "PointwiseFunctions/MathFunctions/MathFunction.hpp"
      16             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      17             : 
      18             : /// \cond
      19             : namespace domain::FunctionsOfTime {
      20             : class FunctionOfTime;
      21             : }  // namespace domain::FunctionsOfTime
      22             : namespace PUP {
      23             : class er;
      24             : }  // namespace PUP
      25             : /// \endcond
      26             : 
      27             : namespace domain::CoordinateMaps::TimeDependent {
      28             : /*!
      29             :  * \ingroup CoordMapsTimeDependentGroup
      30             :  * \brief Translation map defined by \f$\vec{x} = \vec{\xi}+F(r)\vec{T}(t)\f$
      31             :  * where $F(r)$ takes on different forms based on which constructor is used.
      32             :  *
      33             :  * \details The map adds a translation to the coordinates $\vec{\xi}$ based on
      34             :  * what type of translation is needed. For the piecewise translation, a
      35             :  * translation $F(r)\vec{T}(t)$ is added to $\vec{\xi}$ based on what region
      36             :  * $|\vec{\xi}|$ is in. For coordinates within the inner radius, $F(r) = 1$
      37             :  * causing a uniform translation. Coordinates in between the inner and outer
      38             :  * radius have a linear radial falloff applied to them. Coordinates beyond the
      39             :  * outer radius have no translation applied to them $F(r) = 0$. The piecewise
      40             :  * translation assumes that the center of your map is at (0., 0., 0.). For the
      41             :  * radial MathFunction translation, a radial translation \f$F(r)\vec{T}(t)\f$ is
      42             :  * added to the coordinates \f$\vec{\xi}\f$, where \f$\vec{T}(t)\f$ is a
      43             :  * FunctionOfTime and\f$F(r)\f$ is a 1D radial MathFunction. The radius of each
      44             :  * point is found by subtracting the center map argument from the coordinates
      45             :  * \f$\vec{\xi}\f$ or the target coordinates \f$\vec{\bar{\xi}}\f$. The
      46             :  * Translation Map class is overloaded so that the user can choose between a
      47             :  * piecewise translation, radial translation or a uniform translation based on
      48             :  * their problem. If a radial dependence is not specified, this sets \f$F(r) =
      49             :  * 1\f$.
      50             :  *
      51             :  * ### Mapped Coordinates
      52             :  * The piecewise translation translates the coordinates $\vec{\xi}$
      53             :  * to the target coordinates $\vec{\bar{\xi}}$ based on the region $\vec{\xi}$
      54             :  * is in.
      55             :  * \f{equation}{
      56             :  * \vec{\bar{\xi}} = \left\{\begin{array}{ll}\vec{\xi} + \vec{T}(t), &
      57             :  * |\vec{\xi}| \leq R_{in}, \\ \vec{\xi} + wT(t), &  R_{in} < |\vec{\xi}| <
      58             :  * R_{out}, \\ \vec{\xi}, & |\vec{\xi}| \geq R_{out} \end{array}\right.
      59             :  * \f}
      60             :  *
      61             :  * Where $R_{in}$ is the inner radius, $R_{out}$ is the outer radius, and $w$ is
      62             :  * the radial falloff factor found through
      63             :  * \f{equation}{
      64             :  * w = \frac{R_{out} - |\vec{\xi}|}{R_{out} - R_{in}}
      65             :  * \f}
      66             :  *
      67             :  * The radial MathFunction translation translates the coordinates
      68             :  * \f$\vec{\xi}\f$ to the target coordinates \f{equation}{\vec{\bar{\xi}} =
      69             :  * \vec{\xi} + F(r)\vec{T}(t) \f}
      70             :  *
      71             :  * If you only supply a FunctionOfTime to the constructor of this class, the
      72             :  * radial function will be set to 1.0 causing a uniform translation for your
      73             :  * coordinates. If a FunctionOfTime, MathFunction, and map center are passed in,
      74             :  * the radius will be found through
      75             :  * \f{equation}{
      76             :  * r = |\vec{\xi} - \vec{c}|
      77             :  * \f}
      78             :  * where r is the radius and \f$\vec{c}\f$ is the center argument.
      79             :  *
      80             :  * ### Inverse Translation
      81             :  * The piecewise inverse translates the coordinates
      82             :  * \f$\vec{\bar{\xi}}\f$ to the original coordinates based on what region
      83             :  * $\vec{\bar{\xi}}$ is in.
      84             :  * \f{equation}{
      85             :  * \vec{\xi} = \left\{\begin{array}{ll}\vec{\bar{\xi}} -
      86             :  * \vec{T}(t), & |\vec{\bar{\xi}}| \leq R_{in}, or, |\vec{\bar{\xi}} - T(t)|
      87             :  * \leq R_{in}, \\
      88             :  * \vec{\bar{\xi}} - wT(t), &  R_{in} < |\vec{\bar{\xi}}| < R_{out}, \\
      89             :  * \vec{\bar{\xi}}, & |\vec{\bar{\xi}}| \geq R_{out}\end{array}\right.
      90             :  * \f}
      91             :  * Where $w$ is the radial falloff factor found through a quadratic solve of the
      92             :  * form
      93             :  * \f{equation}{
      94             :  * w^2(\vec{T}(t)^2 - (R_{out} - R_{in})^2) - 2w(\vec{T}(t)\vec{\bar{\xi}} -
      95             :  * R_{out}(R_{out} - R_{in})) + \vec{\bar{\xi}}^2 - R_{out}^2
      96             :  * \f}
      97             :  * The inverse map also assumes that if $\vec{\bar{\xi}}
      98             :  * - \vec{T}(t) \leq R_{in}$ then the translated point originally came from
      99             :  * within the inner radius so it'll be translated back without a quadratic
     100             :  * solve.
     101             :  *
     102             :  * The radial MathFunction inverse translates the coordinates
     103             :  * \f$\vec{\bar{\xi}}\f$ to the original coordinates using
     104             :  * \f{equation}{
     105             :  * \vec{\xi} = \vec{\bar{\xi}} - F(r)\vec{T}(t)
     106             :  * \f}
     107             :  * where \f$r^2\f$ is found as the root of
     108             :  * \f{equation}{
     109             :  *   r^2 = \Big(\vec{\bar{\xi}} - \vec{c} - F(r) \vec{T}(t)\Big)^2.
     110             :  * \f}
     111             :  *
     112             :  * ### Frame Velocity
     113             :  * For the piecewise translation, the frame velocity is found through
     114             :  * \f{equation}{
     115             :  * \vec{v} = \left\{\begin{array}{ll}\frac{\vec{dT}(t)}{dt}, & |\vec{\xi}| \leq
     116             :  * R_{in}, \\ w\frac{\vec{dT}(t)}{dt}, &  R_{in} < |\vec{\xi}| < R_{out}, \\ 0,
     117             :  * & |\vec{\xi}| \geq R_{out} \end{array}\right.
     118             :  * \f}
     119             :  *
     120             :  * For the radial MathFunction translation, the frame velocity is found through
     121             :  * \f{equation}{
     122             :  * \vec{v} = \frac{\vec{dT}(t)}{dt} F(r)
     123             :  * \f}
     124             :  * where \f$\frac{\vec{dT}(t)}{dt}\f$ is the first derivative of the
     125             :  * FunctionOfTime.
     126             :  *
     127             :  * ### Jacobian
     128             :  * For the piecewise translation, the jacobian is computed based on what region
     129             :  * the coordinates $\vec{\xi}$ is in.
     130             :  * \f{equation}{
     131             :  * {J^{i}}_{j} = \frac{dw}{dr} T(t)^i \frac{\xi_j}{r}, R_{in} <
     132             :  * |\vec{\bar{\xi}}| < R_{out}
     133             :  * \f}
     134             :  * otherwise, it will return the identity matrix.
     135             :  *
     136             :  * For the radial MathFunction translation, the jacobian is computed through the
     137             :  * first derivative when the radius is bigger than 1.e-13:
     138             :  * \f{equation}{
     139             :  * {J^{i}}_{j} = \frac{dF(r)}{dr} T(t)^i \frac{(\xi_j - c_j)}{r}
     140             :  * \f}
     141             :  * Where \f$\frac{dF(r)}{dr}\f$ is the first derivative of the MathFunction,
     142             :  * \f$\vec{\xi_j}\f$ is the source coordinates, \f$\vec{c}\f$ is the center of
     143             :  * your map, and r is the radius.
     144             :  *
     145             :  * At a radius smaller than 1e-13, we ASSERT that the radial MathFunction is
     146             :  * smooth $\frac{dF(r)}{dr} \approx 0$, so return the identity matrix.
     147             :  *
     148             :  *
     149             :  * ### Inverse Jacobian
     150             :  * The inverse jacobian is computed numerically by inverting the jacobian.
     151             :  */
     152             : template <size_t Dim>
     153           1 : class Translation {
     154             :  public:
     155           0 :   static constexpr size_t dim = Dim;
     156             : 
     157           0 :   Translation() = default;
     158           0 :   explicit Translation(std::string function_of_time_name);
     159             : 
     160           0 :   explicit Translation(std::string function_of_time_name, double inner_radius,
     161             :                        double outer_radius);
     162             : 
     163           0 :   explicit Translation(
     164             :       std::string function_of_time_name,
     165             :       std::unique_ptr<MathFunction<1, Frame::Inertial>> radial_function,
     166             :       std::array<double, Dim>& center);
     167             : 
     168           0 :   Translation(const Translation<Dim>& Translation_Map);
     169             : 
     170           0 :   ~Translation() = default;
     171           0 :   Translation(Translation&&) = default;
     172           0 :   Translation& operator=(Translation&&) = default;
     173           0 :   Translation& operator=(const Translation& Translation_Map);
     174             : 
     175             :   template <typename T>
     176           0 :   std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
     177             :       const std::array<T, Dim>& source_coords, double time,
     178             :       const std::unordered_map<
     179             :           std::string,
     180             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     181             :           functions_of_time) const;
     182             : 
     183             :   /// The inverse function is only callable with doubles because the inverse
     184             :   /// might fail if called for a point out of range, and it is unclear
     185             :   /// what should happen if the inverse were to succeed for some points in a
     186             :   /// DataVector but fail for other points.
     187           1 :   std::optional<std::array<double, Dim>> inverse(
     188             :       const std::array<double, Dim>& target_coords, double time,
     189             :       const std::unordered_map<
     190             :           std::string,
     191             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     192             :           functions_of_time) const;
     193             : 
     194             :   template <typename T>
     195           0 :   std::array<tt::remove_cvref_wrap_t<T>, Dim> frame_velocity(
     196             :       const std::array<T, Dim>& source_coords, double time,
     197             :       const std::unordered_map<
     198             :           std::string,
     199             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     200             :           functions_of_time) const;
     201             : 
     202             :   template <typename T>
     203           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
     204             :       const std::array<T, Dim>& source_coords, double time,
     205             :       const std::unordered_map<
     206             :           std::string,
     207             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     208             :           functions_of_time) const;
     209             : 
     210             :   template <typename T>
     211           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
     212             :       const std::array<T, Dim>& source_coords, double time,
     213             :       const std::unordered_map<
     214             :           std::string,
     215             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     216             :           functions_of_time) const;
     217             : 
     218             :   // NOLINTNEXTLINE(google-runtime-references)
     219           0 :   void pup(PUP::er& p);
     220             : 
     221           0 :   static bool is_identity() { return false; }
     222             : 
     223           0 :   const std::unordered_set<std::string>& function_of_time_names() const {
     224             :     return f_of_t_names_;
     225             :   }
     226             : 
     227             :  private:
     228             :   template <size_t LocalDim>
     229           0 :   friend bool operator==(  // NOLINT(readability-redundant-declaration)
     230             :       const Translation<LocalDim>& lhs, const Translation<LocalDim>& rhs);
     231             : 
     232             :   // These 2 helper functions compute the translated coordinates or frame
     233             :   // velocity based on the option passed in, 0 for translated coordinates, and
     234             :   // frame velocity for any other number.
     235             :   template <typename T>
     236           0 :   std::array<tt::remove_cvref_wrap_t<T>, Dim> math_function_helper(
     237             :       const std::array<T, Dim>& source_coords, double time,
     238             :       const std::unordered_map<
     239             :           std::string,
     240             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     241             :           functions_of_time,
     242             :       size_t function_or_deriv_index) const;
     243             : 
     244             :   template <typename T>
     245           0 :   std::array<tt::remove_cvref_wrap_t<T>, Dim> piecewise_helper(
     246             :       const std::array<T, Dim>& source_coords, double time,
     247             :       const std::unordered_map<
     248             :           std::string,
     249             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
     250             :           functions_of_time,
     251             :       size_t function_or_deriv_index) const;
     252             : 
     253           0 :   double root_finder(const std::array<double, Dim>& distance_to_center,
     254             :                      const DataVector& function_of_time) const;
     255             : 
     256           0 :   std::string f_of_t_name_{};
     257           0 :   std::unordered_set<std::string> f_of_t_names_;
     258           0 :   std::optional<double> inner_radius_;
     259           0 :   std::optional<double> outer_radius_;
     260           0 :   std::unique_ptr<MathFunction<1, Frame::Inertial>> f_of_r_{};
     261           0 :   std::array<double, Dim> center_{};
     262             : };
     263             : 
     264             : template <size_t Dim>
     265           0 : inline bool operator!=(const Translation<Dim>& lhs,
     266             :                        const Translation<Dim>& rhs) {
     267             :   return not(lhs == rhs);
     268             : }
     269             : 
     270             : }  // namespace domain::CoordinateMaps::TimeDependent

Generated by: LCOV version 1.14