SpECTRE Documentation Coverage Report
Current view: top level - Domain/CoordinateMaps/TimeDependent - Shape.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 1 33 3.0 %
Date: 2024-03-29 00:33:31
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 <cstddef>
       7             : #include <memory>
       8             : #include <optional>
       9             : #include <string>
      10             : #include <unordered_set>
      11             : 
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
      14             : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
      15             : #include "Utilities/Gsl.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             : /*!
      30             :  * \ingroup CoordMapsTimeDependentGroup
      31             :  * \brief Distorts a distribution of points radially according to a spherical
      32             :  * harmonic expansion while preserving angles.
      33             :  *
      34             :  * \details The shape map distorts the distance \f$r\f$ between a point and
      35             :  * the center while leaving the angles \f$\theta\f$, \f$\phi\f$ between them
      36             :  * preserved by applying a spherical harmonic expansion with time-dependent
      37             :  * coefficients \f$\lambda_{lm}(t)\f$. There are two ways to specify the
      38             :  * time-dependent coefficients \f$\lambda_{lm}(t)\f$:
      39             :  *
      40             :  * 1. A single FunctionOfTime which specifies all coefficients. This
      41             :  *    FunctionOfTime should have `ylm::Spherepack::spectral_size()` number of
      42             :  *    components. These are in Spherepack order and should be the Spherepack
      43             :  *    coefficients, *not* the spherical harmonic coefficients. See the note
      44             :  *    below. To use this, set the `size_function_of_time_name` argument of the
      45             :  *    constructor to `std::nullopt`.
      46             :  * 2. Two different FunctionOfTime%s. The first is similar to 1.) in that it
      47             :  *    should have the same number of components, be in Spherepack order, and be
      48             :  *    the Spherepack coefficients. The only difference is that the \f$l = 0\f$
      49             :  *    coefficient should be identically 0. The second FunctionOfTime should have
      50             :  *    a single component which will be the \f$l = 0\f$ coefficient. This
      51             :  *    component should be stored as the spherical harmonic coefficient and *not*
      52             :  *    a Spherepack coefficient. See the note below. To use this method, set the
      53             :  *    `size_function_of_time_name` argument of the constructor to the name of
      54             :  *    the FunctionOfTime that's in the cache. This method is useful if we have
      55             :  *    control systems because we have a separate control system controlling a
      56             :  *    separate function of time for the \f$l = 0\f$ coefficient than we do for
      57             :  *    the other coefficients.
      58             :  *
      59             :  * \note The quantities stored in the "shape" FunctionOfTime (the
      60             :  * `shape_function_of_time_name` argument in the constructor that must always be
      61             :  * specified) are ***not*** the complex spherical-harmonic coefficients
      62             :  * \f$\lambda_{lm}(t)\f$, but instead are the real-valued SPHEREPACK
      63             :  * coefficients \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$ used by Spherepack. This
      64             :  * is the same for both methods of specifying FunctionOfTime%s above. The
      65             :  * relationship between these two sets of coefficients is
      66             :  * \f{align}
      67             :  * a_{l0} & = \sqrt{\frac{2}{\pi}}\lambda_{l0}&\qquad l\geq 0,\\
      68             :  * a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(\lambda_{lm})
      69             :  * &\qquad l\geq 1, m\geq 1, \\
      70             :  * b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(\lambda_{lm})
      71             :  * &\qquad l\geq 1, m\geq 1.
      72             :  * \f}
      73             :  * The "shape" FunctionOfTime stores coefficients only for non-negative \f$m\f$;
      74             :  * this is because the function we are expanding is real, so the
      75             :  * coefficients for \f$m<0\f$ can be obtained from \f$m>0\f$ coefficients by
      76             :  * complex conjugation.
      77             :  * If the `size_function_of_time_name` argument is given to the constructor,
      78             :  * then it is asserted that the \f$l=0\f$ coefficient of the "shape" function of
      79             :  * time is exactly 0. The \f$l=0\f$ coefficient is then controlled by the "size"
      80             :  * FunctionOfTime. Unlike the "shape" FunctionOfTime, the quantity in the
      81             :  * "size" FunctionOfTime ***is*** the "complex" spherical harmonic coefficient
      82             :  * \f$\lambda_{00}(t)\f$, and not the SPHEREPACK coefficient \f$a_{00}(t)\f$
      83             :  * ("complex" is in quotes because all \f$m=0\f$ coefficients are always real.)
      84             :  * Here and below we write the equations in terms of \f$\lambda_{lm}(t)\f$
      85             :  * instead of \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$, regardless of which
      86             :  * FunctionOfTime representation we are using, because the resulting expressions
      87             :  * are much shorter.
      88             :  *
      89             :  * An additional dimensionless domain-dependent transition function
      90             :  * \f$f(r, \theta, \phi)\f$ ensures that the distortion falls off correctly to
      91             :  * zero at a certain boundary (must be a block boundary). The function \f$f(r,
      92             :  * \theta, \phi)\f$ is restricted such that
      93             :  *
      94             :  * \f{equation}{
      95             :  * 0 \leq f(r, \theta, \phi) \leq 1
      96             :  * \f}
      97             :  *
      98             :  * ### Mapped coordinates
      99             :  *
     100             :  * Given a point with cartesian coordinates \f$\xi^i\f$, let the polar
     101             :  * coordinates \f$(r, \theta, \phi)\f$ with respect to a center \f$x_c^i\f$ be
     102             :  * defined in the usual way:
     103             :  *
     104             :  * \f{align}{
     105             :  * \xi^0 - x_c^0 &= r \sin(\theta) \cos(\phi)\\
     106             :  * \xi^1 - x_c^1 &= r \sin(\theta) \sin(\phi)\\
     107             :  * \xi^2 - x_c^2 &= r \cos(\theta)
     108             :  * \f}
     109             :  *
     110             :  * The shape map maps the unmapped
     111             :  * coordinates \f$\xi^i\f$ to coordinates \f$x^i\f$:
     112             :  *
     113             :  * \f{equation}{
     114             :  * x^i = \xi^i - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm}
     115             :  * \lambda_{lm}(t)Y_{lm}(\theta, \phi).
     116             :  * \f}
     117             :  *
     118             :  * Or written another way
     119             :  *
     120             :  * \f{equation}{\label{eq:map_center_plus_distortion}
     121             :  * x^i = x_c^i + (\xi^i - x_c^i) \left(1 - \frac{f(r, \theta, \phi)}{r}
     122             :  * \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right).
     123             :  * \f}
     124             :  *
     125             :  * The form in Eq. \f$\ref{eq:map_center_plus_distortion}\f$ makes two things
     126             :  * clearer
     127             :  *
     128             :  * 1. This shape map is just a radial distortion about \f$x_c^i\f$
     129             :  * 2. The coefficients \f$\lambda_{lm}\f$ have units of distance because
     130             :  *    \f$\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r\f$ must be dimensionless
     131             :  *    (because \f$f\f$ is dimensionless).
     132             :  *
     133             :  * ### Inverse map
     134             :  *
     135             :  * The inverse map is given by:
     136             :  * \f{equation}{
     137             :  * \xi^i = x_c^i + (x^i-x_c^i)*(r/\tilde{r}),
     138             :  * \f}
     139             :  * where \f$\tilde{r}\f$ is the radius of \f$\xi\f$, calculated by the
     140             :  * transition map. In order to compute \f$r/\tilde{r}\f$, the following equation
     141             :  * must be solved
     142             :  *
     143             :  * \f{equation}{
     144             :  * \frac{r}{\tilde{r}} =
     145             :  * \frac{1}{1-f(r,\theta,\phi)\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r}
     146             :  * \f}
     147             :  *
     148             :  * For more details, see
     149             :  * \link domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
     150             :  * ShapeMapTransitionFunction::original_radius_over_radius \endlink.
     151             :  *
     152             :  * ### Frame velocity
     153             :  *
     154             :  * The frame velocity \f$v^i\ = dx^i / dt\f$ is calculated trivially:
     155             :  * \f{equation}{
     156             :  * v^i = - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm}
     157             :  * \dot{\lambda}_{lm}(t)Y_{lm}(\theta, \phi).
     158             :  * \f}
     159             :  *
     160             :  * ### Jacobian
     161             :  *
     162             :  * The Jacobian is given by:
     163             :  * \f{align}{
     164             :  * \frac{\partial x^i}{\partial \xi^j} = \delta_j^i &\left( 1 - \frac{f(r,
     165             :  * \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right)
     166             :  * \nonumber \\
     167             :  * &- (\xi^i - x_c^i)
     168             :  * \left[\left(\frac{1}{r}\frac{\partial}{\partial \xi^j} f(r, \theta, \phi) -
     169             :  * \frac{\xi_j}{r^3} f(r, \theta, \phi)\right) \sum_{lm}
     170             :  * \lambda_{lm}(t)Y_{lm}(\theta, \phi) + \frac{f(r, \theta, \phi)}{r}
     171             :  * \sum_{lm} \lambda_{lm}(t) \frac{\partial}{\partial \xi^j} Y_{lm}(\theta,
     172             :  * \phi) \right].
     173             :  * \f}
     174             :  *
     175             :  * where \f$\xi_j = \xi^j\f$.
     176             :  *
     177             :  * ### Inverse Jacobian
     178             :  *
     179             :  * The inverse Jacobian is computed by numerically inverting the Jacobian.
     180             :  *
     181             :  * For future optimization, the `interpolation_info` objects calculated in all
     182             :  * functions of this class could be cached. Since every element should evaluate
     183             :  * the same grid coordinates most time steps, this might greatly decrease
     184             :  * computation. Every element has their own clone of the shape map so the
     185             :  * caching could be done with member variables. Care must be taken that
     186             :  * `jacobian` currently calculates the `interpolation_info` with an order
     187             :  * higher.
     188             :  */
     189           1 : class Shape {
     190             :  public:
     191           0 :   using FunctionsOfTimeMap = std::unordered_map<
     192             :       std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
     193             : 
     194           0 :   explicit Shape(
     195             :       const std::array<double, 3>& center, size_t l_max, size_t m_max,
     196             :       std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
     197             :           transition_func,
     198             :       std::string shape_function_of_time_name,
     199             :       std::optional<std::string> size_function_of_time_name = std::nullopt);
     200             : 
     201           0 :   Shape() = default;
     202           0 :   ~Shape() = default;
     203           0 :   Shape(Shape&&) = default;
     204           0 :   Shape& operator=(Shape&&) = default;
     205           0 :   Shape(const Shape& rhs);
     206           0 :   Shape& operator=(const Shape& rhs);
     207             : 
     208             :   template <typename T>
     209           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
     210             :       const std::array<T, 3>& source_coords, double time,
     211             :       const FunctionsOfTimeMap& functions_of_time) const;
     212             : 
     213           0 :   std::optional<std::array<double, 3>> inverse(
     214             :       const std::array<double, 3>& target_coords, double time,
     215             :       const FunctionsOfTimeMap& functions_of_time) const;
     216             : 
     217             :   template <typename T>
     218           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> frame_velocity(
     219             :       const std::array<T, 3>& source_coords, double time,
     220             :       const FunctionsOfTimeMap& functions_of_time) const;
     221             : 
     222             :   template <typename T>
     223           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
     224             :       const std::array<T, 3>& source_coords, double time,
     225             :       const FunctionsOfTimeMap& functions_of_time) const;
     226             : 
     227             :   template <typename T>
     228           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
     229             :       const std::array<T, 3>& source_coords, double time,
     230             :       const FunctionsOfTimeMap& functions_of_time) const;
     231             : 
     232             :   // NOLINTNEXTLINE(google-runtime-references)
     233           0 :   void pup(PUP::er& p);
     234           0 :   static bool is_identity() { return false; }
     235           0 :   static constexpr size_t dim = 3;
     236             : 
     237           0 :   const std::unordered_set<std::string>& function_of_time_names() const {
     238             :     return f_of_t_names_;
     239             :   }
     240             : 
     241             :  private:
     242           0 :   std::string shape_f_of_t_name_;
     243           0 :   std::optional<std::string> size_f_of_t_name_;
     244           0 :   std::unordered_set<std::string> f_of_t_names_;
     245           0 :   std::array<double, 3> center_{};
     246           0 :   size_t l_max_ = 2;
     247           0 :   size_t m_max_ = 2;
     248           0 :   ylm::Spherepack ylm_{2, 2};
     249             :   std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
     250           0 :       transition_func_;
     251             : 
     252             :   template <typename T>
     253           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> center_coordinates(
     254             :       const std::array<T, 3>& coords) const {
     255             :     return {coords[0] - center_[0], coords[1] - center_[1],
     256             :             coords[2] - center_[2]};
     257             :   }
     258             : 
     259           0 :   void check_size(const gsl::not_null<DataVector*>& coefs,
     260             :                   const FunctionsOfTimeMap& functions_of_time, double time,
     261             :                   bool use_deriv) const;
     262             : 
     263             :   // Checks that the vector of coefficients has the right size and that the
     264             :   // monopole and dipole coefficients are zero.
     265           0 :   void check_coefficients(const DataVector& coefs) const;
     266             : 
     267             :   template <typename T>
     268           0 :   T check_and_compute_one_over_radius(
     269             :       const std::array<T, 3>& centered_coords) const;
     270             : 
     271           0 :   friend bool operator==(const Shape& lhs, const Shape& rhs);
     272             : };
     273           0 : bool operator!=(const Shape& lhs, const Shape& rhs);
     274             : 
     275             : }  // namespace domain::CoordinateMaps::TimeDependent

Generated by: LCOV version 1.14