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 <optional>
       8             : #include <pup.h>
       9             : 
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "Utilities/Serialization/CharmPupable.hpp"
      12             : 
      13             : namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
      14             : 
      15             : /*!
      16             :  * \brief Abstract base class for the transition function $G(r,\theta,\phi)$
      17             :  * used by the domain::CoordinateMaps::TimeDependent::Shape map.
      18             :  *
      19             :  * \details This base class defines the required methods of a transition
      20             :  * function that is defined as
      21             :  *
      22             :  * \begin{equation}
      23             :  * G(r,\theta,\phi) = \frac{f(r,\theta,\phi)}{r}
      24             :  * \end{equation}
      25             :  *
      26             :  * where $0 <= f(r,\theta,\phi) <= 1$. Different domains may require the shape
      27             :  * map to fall off towards the boundary in different ways. This behavior is
      28             :  * controlled by the transition function. It is also needed to find the inverse
      29             :  * of the shape map. Since the shape map preserves angles, the problem of
      30             :  * finding its inverse reduces to the 1-dimensional problem of finding the
      31             :  * original radius from the mapped radius. The mapped radius $\tilde{r}$ is
      32             :  * related to the original $r$ radius by:
      33             :  * \begin{equation}
      34             :  * \label{eq:shape_map_radius}
      35             :  * \tilde{r} = r \left(1 - G(r,\theta,\phi)\sum_{lm}
      36             :  *     \lambda_{lm}(t)Y_{lm}(\theta,\phi)\right),
      37             :  * \end{equation}
      38             :  * Depending on the format of the transition function, it should be possible to
      39             :  * analytically derive this map's inverse because it preserves angles and shifts
      40             :  * only the radius of each point. Otherwise the inverse has to be computed
      41             :  * numerically.
      42             :  *
      43             :  * The transition function must also be able to compute the gradient and the
      44             :  * value of the function divided by the radius. Care must be taken that this
      45             :  * does not divide by zero.
      46             :  *
      47             :  * All member functions with the exception of `original_radius_over_radius`
      48             :  * exist as overloads for types `double` and `DataVector` so that they work with
      49             :  * the templated shape map methods calling them. To avoid code duplication these
      50             :  * can be forwarded to templated implementation methods held by the derived
      51             :  * classes only.
      52             :  *
      53             :  * For an example, see SphereTransition.
      54             :  *
      55             :  * #### Design Decisions:
      56             :  *
      57             :  * It was decided to make the ShapeMapTransitionFunction an abstract base class
      58             :  * with overloads for types `double` and `DataVector` corresponding to the
      59             :  * template parameter `T` of the shape map's methods. The shape map holds a
      60             :  * `unique_ptr` to this abstract class using a common dynamic dispatch design
      61             :  * pattern. This approach avoids templating the shape map all together.
      62             :  *
      63             :  * An alternative approach would be to directly template the transition
      64             :  * functions onto the shape map so that no abstract base class is necessary.
      65             :  * These approaches can also be combined by making the transition function an
      66             :  * abstract base class but also templating it onto the shape map. In this way
      67             :  * the shape map does not need to hold a `unique_ptr` but can hold the
      68             :  * transition function directly as a member.
      69             :  */
      70           1 : class ShapeMapTransitionFunction : public PUP::able {
      71             :  public:
      72           0 :   ShapeMapTransitionFunction() = default;
      73             : 
      74             :   /// @{
      75             :   /*!
      76             :    * Evaluate the transition function $G(r,\theta,\phi)$ at the
      77             :    * Cartesian coordinates `source_coords` and possibly divide by the radius
      78             :    * $r$.
      79             :    *
      80             :    * \details If \p one_over_radius_power has a value, then divide $G$ by $r$ to
      81             :    * that power. This is done here because the transition function knows how to
      82             :    * handle points in various regions of the map.
      83             :    */
      84           1 :   virtual double operator()(
      85             :       const std::array<double, 3>& source_coords,
      86             :       const std::optional<size_t>& one_over_radius_power) const = 0;
      87           1 :   virtual DataVector operator()(
      88             :       const std::array<DataVector, 3>& source_coords,
      89             :       const std::optional<size_t>& one_over_radius_power) const = 0;
      90             :   /// @}
      91             : 
      92             :   /*!
      93             :    * \brief The inverse of the transition function
      94             :    *
      95             :    * This method returns $r/\tilde{r}$ given the mapped coordinates
      96             :    * $\tilde{x}^i$ (`target_coords`) and the spherical harmonic expansion
      97             :    * $\Sigma(t, \theta, \phi) = \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)$
      98             :    * (`radial_distortion`). See domain::CoordinateMaps::TimeDependent::Shape for
      99             :    * details on how this quantity is used to compute the inverse of the Shape
     100             :    * map.
     101             :    *
     102             :    * To derive the expression for this inverse, solve Eq.
     103             :    * ($\ref{eq:shape_map_radius}$) for $r$ after substituting
     104             :    * $G(r,\theta,\phi)$.
     105             :    *
     106             :    * \param target_coords The mapped Cartesian coordinates $\tilde{x}^i$.
     107             :    * \param radial_distortion The spherical harmonic expansion
     108             :    * $\Sigma(t, \theta, \phi)$.
     109             :    * \return The quantity $r/\tilde{r}$.
     110             :    */
     111           1 :   virtual std::optional<double> original_radius_over_radius(
     112             :       const std::array<double, 3>& target_coords,
     113             :       double radial_distortion) const = 0;
     114             : 
     115             :   /*!
     116             :    * Evaluate the gradient of the transition function with respect to the
     117             :    * Cartesian coordinates x, y and z at the Cartesian coordinates
     118             :    * `source_coords`.
     119             :    */
     120             :   /// @{
     121           1 :   virtual std::array<double, 3> gradient(
     122             :       const std::array<double, 3>& source_coords) const = 0;
     123           1 :   virtual std::array<DataVector, 3> gradient(
     124             :       const std::array<DataVector, 3>& source_coords) const = 0;
     125             :   /// @}
     126             : 
     127           0 :   virtual std::unique_ptr<ShapeMapTransitionFunction> get_clone() const = 0;
     128             : 
     129           0 :   virtual bool operator==(const ShapeMapTransitionFunction& other) const = 0;
     130           0 :   virtual bool operator!=(const ShapeMapTransitionFunction& other) const = 0;
     131             : 
     132           0 :   WRAPPED_PUPable_abstract(ShapeMapTransitionFunction);
     133           0 :   explicit ShapeMapTransitionFunction(CkMigrateMessage* m) : PUP::able(m) {}
     134             : };
     135             : }  // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions