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/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Elliptic/Systems/Punctures/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
18 : #include "Options/String.hpp"
19 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
20 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
21 : #include "Utilities/Gsl.hpp"
22 : #include "Utilities/TMPL.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
|
