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 functions used by the
      17             :  * domain::CoordinateMaps::TimeDependent::Shape map.
      18             :  *
      19             :  * \details This base class defines the required methods of a transition
      20             :  * function used by the shape map. Different domains require the shape map to
      21             :  * fall off towards the boundary in different ways. This behavior is controlled
      22             :  * by the transition function. It is also needed to find the inverse of the
      23             :  * shape map. Since the shape map preserves angles, the problem of finding its
      24             :  * inverse reduces to the 1-dimensional problem of finding the original radius
      25             :  * from the mapped radius. The mapped radius \f$\tilde{r}\f$ is related to the
      26             :  * original \f$r\f$ radius by:
      27             :  * \f{equation}{\label{eq:shape_map_radius}
      28             :  *  \tilde{r} = r (1 - f(r,\theta,\phi) \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta,
      29             :  *  \phi)),
      30             :  *  \f}
      31             :  * where \f$f(r,\theta,\phi)\f$ is the transition function. Depending
      32             :  * on the format of the transition function, it should be possible to
      33             :  * analytically derive this map's inverse because it preserves angles and shifts
      34             :  * only the radius of each point. Otherwise the inverse has to be computed
      35             :  * numerically.
      36             :  *
      37             :  * The transition function must also be able to compute the gradient and the
      38             :  * value of the function divided by the radius. Care must be taken that this
      39             :  * does not divide by zero.
      40             :  *
      41             :  * All member functions with the exception of `original_radius_over_radius`
      42             :  * exist as overloads for types `double` and `DataVector` so that they work with
      43             :  * the templated shape map methods calling them. To avoid code duplication these
      44             :  * can be forwarded to templated implementation methods held by the derived
      45             :  * classes only.
      46             :  *
      47             :  * For an example, see SphereTransition.
      48             :  *
      49             :  * #### Design Decisions:
      50             :  *
      51             :  * It was decided to make the ShapeMapTransitionFunction an abstract base class
      52             :  * with overloads for types `double` and `DataVector` corresponding to the
      53             :  * template parameter `T` of the shape map's methods. The shape map holds a
      54             :  * `unique_ptr` to this abstract class using a common dynamic dispatch design
      55             :  * pattern. This approach avoids templating the shape map all together.
      56             :  *
      57             :  * An alternative approach would be to directly template the transition
      58             :  * functions onto the shape map so that no abstract base class is necessary.
      59             :  * These approaches can also be combined by making the transition function an
      60             :  * abstract base class but also templating it onto the shape map. In this way
      61             :  * the shape map does not need to hold a `unique_ptr` but can hold the
      62             :  * transition function directly as a member.
      63             :  */
      64           1 : class ShapeMapTransitionFunction : public PUP::able {
      65             :  public:
      66           0 :   ShapeMapTransitionFunction() = default;
      67             : 
      68             :   /// @{
      69             :   /*!
      70             :    * Evaluate the transition function at the Cartesian coordinates
      71             :    *`source_coords`.
      72             :    */
      73           1 :   virtual double operator()(
      74             :       const std::array<double, 3>& source_coords) const = 0;
      75           1 :   virtual DataVector operator()(
      76             :       const std::array<DataVector, 3>& source_coords) const = 0;
      77             :   /// @}
      78             : 
      79             :   /*!
      80             :    * Given the mapped coordinates `target_coords` and the corresponding
      81             :    * spherical harmonic expansion \f$\sum_{lm} \lambda_{lm}(t)Y_{lm}\f$,
      82             :    * `distorted_radius`, this method evaluates the original radius from the
      83             :    * mapped radius by inverting the domain::CoordinateMaps::TimeDependent::Shape
      84             :    * map. It also divides by the mapped radius to simplify calculations in the
      85             :    * shape map.
      86             :    */
      87           1 :   virtual std::optional<double> original_radius_over_radius(
      88             :       const std::array<double, 3>& target_coords,
      89             :       double distorted_radius) const = 0;
      90             : 
      91             :   /*!
      92             :    * Evaluate the gradient of the transition function with respect to the
      93             :    * Cartesian coordinates x, y and z at the Cartesian coordinates
      94             :    * `source_coords`.
      95             :    */
      96             :   /// @{
      97           1 :   virtual std::array<double, 3> gradient(
      98             :       const std::array<double, 3>& source_coords) const = 0;
      99           1 :   virtual std::array<DataVector, 3> gradient(
     100             :       const std::array<DataVector, 3>& source_coords) const = 0;
     101             : 
     102           1 :   virtual std::unique_ptr<ShapeMapTransitionFunction> get_clone() const = 0;
     103             : 
     104             :   /// @}
     105           0 :   virtual bool operator==(const ShapeMapTransitionFunction& other) const = 0;
     106           0 :   virtual bool operator!=(const ShapeMapTransitionFunction& other) const = 0;
     107             : 
     108           0 :   WRAPPED_PUPable_abstract(ShapeMapTransitionFunction);
     109           0 :   explicit ShapeMapTransitionFunction(CkMigrateMessage* m) : PUP::able(m) {}
     110             : };
     111             : }  // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions