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 <pup.h>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/CachedTempBuffer.hpp"
      13             : #include "DataStructures/DataBox/Prefixes.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "Elliptic/Systems/Punctures/Tags.hpp"
      16             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
      19             : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
      20             : #include "Utilities/Gsl.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : #include "Utilities/TaggedTuple.hpp"
      23             : 
      24           0 : namespace Punctures::AnalyticData {
      25             : 
      26             : /// A puncture representing a black hole
      27           1 : struct Puncture {
      28           0 :   struct Position {
      29           0 :     using type = std::array<double, 3>;
      30           0 :     static constexpr Options::String help{"The position C of the puncture"};
      31             :   };
      32           0 :   struct Mass {
      33           0 :     using type = double;
      34           0 :     static constexpr Options::String help{"The puncture mass (bare mass) M"};
      35           0 :     static double lower_bound() { return 0.; }
      36             :   };
      37           0 :   struct Momentum {
      38           0 :     using type = std::array<double, 3>;
      39           0 :     static constexpr Options::String help{
      40             :         "The dimensionless linear momentum P / M, where M is the bare mass of "
      41             :         "the puncture."};
      42             :   };
      43           0 :   struct Spin {
      44           0 :     using type = std::array<double, 3>;
      45           0 :     static constexpr Options::String help{
      46             :         "The dimensionless angular momentum S / M^2, where M is the bare mass "
      47             :         "of the puncture."};
      48             :   };
      49           0 :   using options = tmpl::list<Position, Mass, Momentum, Spin>;
      50           0 :   static constexpr Options::String help{"A puncture representing a black hole"};
      51             : 
      52           0 :   Puncture();
      53           0 :   Puncture(std::array<double, 3> position_in, double mass_in,
      54             :            std::array<double, 3> dimensionless_momentum_in,
      55             :            std::array<double, 3> dimensionless_spin_in);
      56           0 :   void pup(PUP::er& p);
      57             : 
      58           0 :   std::array<double, 3> position{
      59             :       {std::numeric_limits<double>::signaling_NaN()}};
      60           0 :   double mass = std::numeric_limits<double>::signaling_NaN();
      61           0 :   std::array<double, 3> dimensionless_momentum{
      62             :       {std::numeric_limits<double>::signaling_NaN()}};
      63           0 :   std::array<double, 3> dimensionless_spin{
      64             :       {std::numeric_limits<double>::signaling_NaN()}};
      65             : };
      66             : 
      67           0 : bool operator==(const Puncture& lhs, const Puncture& rhs);
      68             : 
      69           0 : bool operator!=(const Puncture& lhs, const Puncture& rhs);
      70             : 
      71             : namespace detail {
      72             : 
      73             : namespace Tags {
      74             : struct OneOverAlpha : db::SimpleTag {
      75             :   using type = Scalar<DataVector>;
      76             : };
      77             : }  // namespace Tags
      78             : 
      79             : struct MultiplePuncturesVariables {
      80             :   static constexpr size_t Dim = 3;
      81             : 
      82             :   using Cache =
      83             :       CachedTempBuffer<detail::Tags::OneOverAlpha, Punctures::Tags::Alpha,
      84             :                        Punctures::Tags::TracelessConformalExtrinsicCurvature,
      85             :                        Punctures::Tags::Beta,
      86             :                        ::Tags::FixedSource<Punctures::Tags::Field>>;
      87             : 
      88             :   const tnsr::I<DataVector, 3>& x;
      89             :   const std::vector<Puncture>& punctures;
      90             : 
      91             :   void operator()(gsl::not_null<Scalar<DataVector>*> one_over_alpha,
      92             :                   gsl::not_null<Cache*> cache,
      93             :                   detail::Tags::OneOverAlpha /*meta*/) const;
      94             :   void operator()(gsl::not_null<Scalar<DataVector>*> alpha,
      95             :                   gsl::not_null<Cache*> cache,
      96             :                   Punctures::Tags::Alpha /*meta*/) const;
      97             :   void operator()(
      98             :       gsl::not_null<tnsr::II<DataVector, 3>*>
      99             :           traceless_conformal_extrinsic_curvature,
     100             :       gsl::not_null<Cache*> cache,
     101             :       Punctures::Tags::TracelessConformalExtrinsicCurvature /*meta*/) const;
     102             :   void operator()(gsl::not_null<Scalar<DataVector>*> beta,
     103             :                   gsl::not_null<Cache*> cache,
     104             :                   Punctures::Tags::Beta /*meta*/) const;
     105             :   void operator()(gsl::not_null<Scalar<DataVector>*> fixed_source_for_field,
     106             :                   gsl::not_null<Cache*> cache,
     107             :                   ::Tags::FixedSource<Punctures::Tags::Field> /*meta*/) const;
     108             : };
     109             : }  // namespace detail
     110             : 
     111             : /*!
     112             :  * \brief Superposition of multiple punctures
     113             :  *
     114             :  * This class provides the source fields $\alpha$ and $\beta$ for the puncture
     115             :  * equation (see \ref ::Punctures) representing any number of black holes. Each
     116             :  * black hole is characterized by its "puncture mass" (or "bare mass") $M_I$,
     117             :  * position $\mathbf{C}_I$, linear momentum $\mathbf{P}_I$, and angular momentum
     118             :  * $\mathbf{S}_I$. The corresponding Bowen-York solution to the momentum
     119             :  * constraint for the conformal traceless extrinsic curvature is:
     120             :  *
     121             :  * \begin{equation}
     122             :  * \bar{A}^{ij} = \frac{3}{2} \sum_I \frac{1}{r_I^2} \left(
     123             :  * 2 P_I^{(i} n_I^{j)} - (\delta^{ij} - n_I^i n_I^j) P_I^k n_I^k
     124             :  * + \frac{4}{r_I} n_I^{(i} \epsilon^{j)kl} S_I^k n_I^l\right)
     125             :  * \end{equation}
     126             :  *
     127             :  * From it, we compute $\alpha$ and $\beta$ as:
     128             :  *
     129             :  * \begin{align}
     130             :  * \frac{1}{\alpha} &= \sum_I \frac{M_I}{2 r_I} \\
     131             :  * \beta &= \frac{1}{8} \alpha^7 \bar{A}_{ij} \bar{A}^{ij}
     132             :  * \end{align}
     133             :  *
     134             :  * \see ::Punctures
     135             :  */
     136           1 : class MultiplePunctures : public elliptic::analytic_data::Background,
     137             :                           public elliptic::analytic_data::InitialGuess {
     138             :  public:
     139           0 :   struct Punctures {
     140           0 :     static constexpr Options::String help =
     141             :         "Parameters for each puncture, representing black holes";
     142           0 :     using type = std::vector<Puncture>;
     143             :   };
     144           0 :   using options = tmpl::list<Punctures>;
     145           0 :   static constexpr Options::String help = "Any number of black holes";
     146             : 
     147           0 :   MultiplePunctures() = default;
     148           0 :   MultiplePunctures(const MultiplePunctures&) = default;
     149           0 :   MultiplePunctures& operator=(const MultiplePunctures&) = default;
     150           0 :   MultiplePunctures(MultiplePunctures&&) = default;
     151           0 :   MultiplePunctures& operator=(MultiplePunctures&&) = default;
     152           0 :   ~MultiplePunctures() = default;
     153             : 
     154           0 :   MultiplePunctures(std::vector<Puncture> punctures)
     155             :       : punctures_(std::move(punctures)) {}
     156             : 
     157           0 :   explicit MultiplePunctures(CkMigrateMessage* m)
     158             :       : elliptic::analytic_data::Background(m),
     159             :         elliptic::analytic_data::InitialGuess(m) {}
     160             :   using PUP::able::register_constructor;
     161           0 :   WRAPPED_PUPable_decl_template(MultiplePunctures);
     162             : 
     163             :   template <typename... RequestedTags>
     164           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     165             :       const tnsr::I<DataVector, 3, Frame::Inertial>& x,
     166             :       tmpl::list<RequestedTags...> /*meta*/) const {
     167             :     typename detail::MultiplePuncturesVariables::Cache cache{
     168             :         get_size(*x.begin())};
     169             :     const detail::MultiplePuncturesVariables computer{x, punctures_};
     170             :     return {cache.get_var(computer, RequestedTags{})...};
     171             :   }
     172             : 
     173             :   template <typename... RequestedTags>
     174           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     175             :       const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& /*mesh*/,
     176             :       const InverseJacobian<DataVector, 3, Frame::ElementLogical,
     177             :                             Frame::Inertial>& /*inv_jacobian*/,
     178             :       tmpl::list<RequestedTags...> /*meta*/) const {
     179             :     return variables(x, tmpl::list<RequestedTags...>{});
     180             :   }
     181             : 
     182             :   // NOLINTNEXTLINE
     183           0 :   void pup(PUP::er& p) override {
     184             :     elliptic::analytic_data::Background::pup(p);
     185             :     elliptic::analytic_data::InitialGuess::pup(p);
     186             :     p | punctures_;
     187             :   }
     188             : 
     189           0 :   const std::vector<Puncture>& punctures() const { return punctures_; }
     190             : 
     191             :  private:
     192           0 :   std::vector<Puncture> punctures_{};
     193             : };
     194             : 
     195             : }  // namespace Punctures::AnalyticData