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