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