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 <limits>
8 : #include <pup.h>
9 :
10 : #include "DataStructures/CachedTempBuffer.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "DataStructures/Tags/TempTensor.hpp"
13 : #include "DataStructures/Tensor/IndexType.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Options/Context.hpp"
16 : #include "Options/String.hpp"
17 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
18 : #include "PointwiseFunctions/AnalyticData/GeneralRelativity/AnalyticData.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : namespace Tags {
24 : template <typename Tag>
25 : struct dt;
26 : } // namespace Tags
27 : namespace gsl {
28 : template <class T>
29 : class not_null;
30 : } // namespace gsl
31 : /// \endcond
32 :
33 : namespace gr::AnalyticData {
34 : /*!
35 : * \brief Brill Lindquist data \cite Brill1963yv corresponding to two black
36 : * holes momentarily at rest
37 : *
38 : * The spatial metric is given by \f$\gamma_{ij} = \psi^4 \delta_{ij}\f$
39 : * where the conformal factor is given by
40 : * \f$\psi = 1 + \frac{m_A}{2 r_A} + \frac{m_B}{2 r_B}\f$ where
41 : * \f$m_{A,B}\f$ are the masses of the black holes and \f$r_{A,B}\f$ are the
42 : * positions of a point relative to the center of each black hole
43 : *
44 : * The data is time symmetric (\f$K_{ij} = 0\f$) and we arbitrarily choose
45 : * unit lapse and zero shift.
46 : */
47 1 : class BrillLindquist : public AnalyticDataBase<3>, public MarkAsAnalyticData {
48 : public:
49 0 : struct MassA {
50 0 : using type = double;
51 0 : static constexpr Options::String help = {"Mass of the black hole A"};
52 0 : static type lower_bound() { return 0.; }
53 : };
54 0 : struct MassB {
55 0 : using type = double;
56 0 : static constexpr Options::String help = {"Mass of the black hole B"};
57 0 : static type lower_bound() { return 0.; }
58 : };
59 0 : struct CenterA {
60 0 : using type = std::array<double, 3>;
61 0 : static constexpr Options::String help = {
62 : "The [x,y,z] center of the black hole A"};
63 : };
64 0 : struct CenterB {
65 0 : using type = std::array<double, 3>;
66 0 : static constexpr Options::String help = {
67 : "The [x,y,z] center of the black hole B"};
68 : };
69 0 : using options = tmpl::list<MassA, MassB, CenterA, CenterB>;
70 0 : static constexpr Options::String help{
71 : "Brill-Lindquist data for two black holes"};
72 :
73 0 : BrillLindquist(double mass_a, double mass_b,
74 : const std::array<double, 3>& center_a,
75 : const std::array<double, 3>& center_b,
76 : const Options::Context& context = {});
77 0 : explicit BrillLindquist(CkMigrateMessage* /*unused*/);
78 :
79 0 : BrillLindquist() = default;
80 0 : BrillLindquist(const BrillLindquist& /*rhs*/) = default;
81 0 : BrillLindquist& operator=(const BrillLindquist& /*rhs*/) = default;
82 0 : BrillLindquist(BrillLindquist&& /*rhs*/) = default;
83 0 : BrillLindquist& operator=(BrillLindquist&& /*rhs*/) = default;
84 0 : ~BrillLindquist() = default;
85 :
86 : /*!
87 : * \brief Computes and returns spacetime quantities for BrillLindquist data at
88 : * a specific Cartesian position
89 : *
90 : * \param x Cartesian coordinates of the position at which to compute
91 : * spacetime quantities
92 : */
93 : template <typename DataType, typename Frame, typename... Tags>
94 1 : tuples::TaggedTuple<Tags...> variables(
95 : const tnsr::I<DataType, volume_dim, Frame>& x,
96 : tmpl::list<Tags...> /*meta*/) const {
97 : static_assert(
98 : tmpl2::flat_all_v<
99 : tmpl::list_contains_v<tags<DataType, Frame>, Tags>...>,
100 : "At least one of the requested tags is not supported. The requested "
101 : "tags are listed as template parameters of the `variables` function.");
102 : IntermediateVars<DataType, Frame> cache(get_size(*x.begin()));
103 : IntermediateComputer<DataType, Frame> computer(*this, x);
104 : return {cache.get_var(computer, Tags{})...};
105 : }
106 :
107 : // NOLINTNEXTLINE(google-runtime-references)
108 0 : void pup(PUP::er& p);
109 :
110 : /*!
111 : * \brief Return the mass of black hole A
112 : */
113 1 : SPECTRE_ALWAYS_INLINE double mass_a() const { return mass_a_; }
114 : /*!
115 : * \brief Return the mass of black hole B
116 : */
117 1 : SPECTRE_ALWAYS_INLINE double mass_b() const { return mass_b_; }
118 : /*!
119 : * \brief Return the center of black hole A
120 : */
121 1 : SPECTRE_ALWAYS_INLINE const std::array<double, 3>& center_a() const {
122 : return center_a_;
123 : }
124 : /*!
125 : * \brief Return the center of black hole B
126 : */
127 1 : SPECTRE_ALWAYS_INLINE const std::array<double, 3>& center_b() const {
128 : return center_b_;
129 : }
130 :
131 : /*!
132 : * \brief Tags defined for intermediates specific to BrillLindquist data
133 : */
134 1 : struct internal_tags {
135 : /*!
136 : * \brief Tag for the position of a point relative to the center of
137 : * black hole A
138 : *
139 : * \details Defined as \f$X_A^i = \left(x^i - C_A^i\right)\f$, where
140 : * \f$C_A^i\f$ is the Cartesian coordinates of the center of black hole A
141 : * and \f$x^i\f$ is the Cartesian coordinates of the point where we're
142 : * wanting to compute spacetime quantities.
143 : */
144 : template <typename DataType, typename Frame>
145 1 : using x_minus_center_a = ::Tags::TempI<0, 3, Frame, DataType>;
146 :
147 : /*!
148 : * \brief Tag for the radius corresponding to the position of a point
149 : * relative to the center of black hole A
150 : *
151 : * \details Defined as \f$r_A = \sqrt{\delta_{ij} X_A^i X_A^j}\f$, where
152 : * \f$X_A^i\f$ is defined by `internal_tags::x_minus_center_a`.
153 : */
154 : template <typename DataType>
155 1 : using r_a = ::Tags::TempScalar<1, DataType>;
156 :
157 : /*!
158 : * \brief Tag for the position of a point relative to the center of
159 : * black hole B
160 : *
161 : * \details Defined as \f$X_B^i = \left(x^i - C_B^i\right)\f$, where
162 : * \f$C_B^i\f$ is the Cartesian coordinates of the center of black hole B
163 : * and \f$x^i\f$ is the Cartesian coordinates of the point where we're
164 : * wanting to compute spacetime quantities.
165 : */
166 : template <typename DataType, typename Frame>
167 1 : using x_minus_center_b = ::Tags::TempI<2, 3, Frame, DataType>;
168 :
169 : /*!
170 : * \brief Tag for the radius corresponding to the position of a point
171 : * relative to the center of black hole B
172 : *
173 : * \details Defined as \f$r_B = \sqrt{\delta_{ij} X_B^i X_B^j}\f$, where
174 : * \f$X_B^i\f$ is defined by `internal_tags::x_minus_center_b`.
175 : */
176 : template <typename DataType>
177 1 : using r_b = ::Tags::TempScalar<3, DataType>;
178 : /*!
179 : * \brief Tag for the conformal factor
180 : *
181 : * \details Defined as \f$\psi = 1 + \frac{m_A}{2 r_A} + \frac{m_B}{2
182 : * r_B}\f$ where \f$m_{A,B}\f$ are the masses of the black holes and
183 : * \f$r_{A,B}\f$ are the positions of a point relative to the center of each
184 : * black hole
185 : */
186 : template <typename DataType>
187 1 : using conformal_factor = ::Tags::TempScalar<4, DataType>;
188 : /*!
189 : * \brief Tag for the deriatives of the conformal factor
190 : *
191 : * \details Defined as \f$d_i\psi = -\frac{m_A X_A^j}{2 r_A^3} \delta_{ij}
192 : * - \frac{m_B X_B^j}{2 r_B^3} \delta_{ij}\f$ where \f$m_{A,B}\f$ are the
193 : * masses of the black holes and \f$r_{A,B}\f$ are the positions of a point
194 : * relative to the center of each black hole. (Note we are free to
195 : * raise/lower coordinate indices with a Eucledian metric)
196 : */
197 : template <typename DataType, typename Frame>
198 1 : using deriv_conformal_factor = ::Tags::Tempi<5, 3, Frame, DataType>;
199 : };
200 :
201 : template <typename DataType, typename Frame>
202 0 : using DerivSpatialMetric =
203 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, volume_dim, Frame>,
204 : tmpl::size_t<volume_dim>, Frame>;
205 :
206 : /*!
207 : * \brief Buffer for caching computed intermediates and quantities that we do
208 : * not want to recompute across the solution's implementation
209 : *
210 : * \details See `internal_tags` documentation for details on what quantities
211 : * the internal tags represent
212 : */
213 : template <typename DataType, typename Frame>
214 1 : using CachedBuffer =
215 : CachedTempBuffer<internal_tags::x_minus_center_a<DataType, Frame>,
216 : internal_tags::r_a<DataType>,
217 : internal_tags::x_minus_center_b<DataType, Frame>,
218 : internal_tags::r_b<DataType>,
219 : internal_tags::conformal_factor<DataType>,
220 : internal_tags::deriv_conformal_factor<DataType, Frame>,
221 : gr::Tags::SpatialMetric<DataType, 3, Frame>,
222 : DerivSpatialMetric<DataType, Frame>>;
223 :
224 : /*!
225 : * \brief Computes the intermediates and quantities that we do not want to
226 : * recompute across the solution's implementation
227 : */
228 : template <typename DataType, typename Frame>
229 1 : class IntermediateComputer {
230 : public:
231 0 : using CachedBuffer = BrillLindquist::CachedBuffer<DataType, Frame>;
232 :
233 : /*!
234 : * \brief Constructs a computer for spacetime quantities of a given
235 : * `gr::AnalyticData::BrillLindquist` solution at at a specific
236 : * Cartesian position
237 : *
238 : * \param analytic_data the given `gr::AnalyticData::BrillLindquist` data
239 : * \param x Cartesian coordinates of the position at which to compute
240 : * spacetime quantities
241 : */
242 1 : IntermediateComputer(const BrillLindquist& analytic_data,
243 : const tnsr::I<DataType, 3, Frame>& x);
244 :
245 : /*!
246 : * \brief Computes the intermediate defined by
247 : * `internal_tags::x_minus_center_a`
248 : */
249 1 : void operator()(
250 : gsl::not_null<tnsr::I<DataType, 3, Frame>*> x_minus_center_a,
251 : gsl::not_null<CachedBuffer*> /*cache*/,
252 : internal_tags::x_minus_center_a<DataType, Frame> /*meta*/) const;
253 :
254 : /*!
255 : * \brief Computes the radius defined by `internal_tags::r_a`
256 : */
257 1 : void operator()(gsl::not_null<Scalar<DataType>*> r_a,
258 : gsl::not_null<CachedBuffer*> cache,
259 : internal_tags::r_a<DataType> /*meta*/) const;
260 :
261 : /*!
262 : * \brief Computes the intermediate defined by
263 : * `internal_tags::x_minus_center_b`
264 : */
265 1 : void operator()(
266 : gsl::not_null<tnsr::I<DataType, 3, Frame>*> x_minus_center_b,
267 : gsl::not_null<CachedBuffer*> /*cache*/,
268 : internal_tags::x_minus_center_b<DataType, Frame> /*meta*/) const;
269 :
270 : /*!
271 : * \brief Computes the radius defined by `internal_tags::r_b`
272 : */
273 1 : void operator()(gsl::not_null<Scalar<DataType>*> r_b,
274 : gsl::not_null<CachedBuffer*> cache,
275 : internal_tags::r_b<DataType> /*meta*/) const;
276 :
277 : /*!
278 : * \brief Computes the intermediate defined by
279 : * `internal_tags::conformal_factor`
280 : */
281 1 : void operator()(gsl::not_null<Scalar<DataType>*> conformal_factor,
282 : gsl::not_null<CachedBuffer*> /*cache*/,
283 : internal_tags::conformal_factor<DataType> /*meta*/) const;
284 :
285 : /*!
286 : * \brief Computes the intermediate defined by
287 : * `internal_tags::deriv_conformal_factor`
288 : */
289 1 : void operator()(
290 : gsl::not_null<tnsr::i<DataType, 3, Frame>*> deriv_conformal_factor,
291 : gsl::not_null<CachedBuffer*> /*cache*/,
292 : internal_tags::deriv_conformal_factor<DataType, Frame> /*meta*/) const;
293 :
294 : /*!
295 : * \brief Computes the spatial metric
296 : *
297 : * \details Computed as \f$\gamma_{ij} = \delta_{ij} \psi^4\f$
298 : * where \f$\psi\f$ is the conformal factor defined by
299 : * `internal_tags::conformal_factor`.
300 : */
301 1 : void operator()(gsl::not_null<tnsr::ii<DataType, 3, Frame>*> spatial_metric,
302 : gsl::not_null<CachedBuffer*> cache,
303 : gr::Tags::SpatialMetric<DataType, 3, Frame> /*meta*/) const;
304 :
305 : /*!
306 : * \brief Computes the spatial derivative of the spatial metric
307 : *
308 : * \details Computed as \f$\partial_k \gamma_{ij} = 4 \psi^3 \partial_k
309 : * \psi \delta_{ij} \f$ where \f$\psi\f$ is the conformal factor defined by
310 : * `internal_tags::conformal_factor`.
311 : */
312 1 : void operator()(
313 : gsl::not_null<tnsr::ijj<DataType, 3, Frame>*> deriv_spatial_metric,
314 : gsl::not_null<CachedBuffer*> cache,
315 : DerivSpatialMetric<DataType, Frame> /*meta*/) const;
316 :
317 : private:
318 : /*!
319 : * \brief The BrillLindquist data
320 : */
321 1 : const BrillLindquist& analytic_data_;
322 : /*!
323 : * \brief Cartesian coordinates of the position at which to compute
324 : * spacetime quantities
325 : */
326 1 : const tnsr::I<DataType, 3, Frame>& x_;
327 : };
328 :
329 : /*!
330 : * \brief Computes and returns spacetime quantities of interest
331 : */
332 : template <typename DataType, typename Frame>
333 1 : class IntermediateVars : public CachedBuffer<DataType, Frame> {
334 : public:
335 0 : using CachedBuffer = BrillLindquist::CachedBuffer<DataType, Frame>;
336 : using CachedBuffer::CachedBuffer;
337 1 : using CachedBuffer::get_var;
338 :
339 : /*!
340 : * \brief Returns the lapse, which is 1
341 : */
342 1 : Scalar<DataType> get_var(
343 : const IntermediateComputer<DataType, Frame>& computer,
344 : gr::Tags::Lapse<DataType> /*meta*/);
345 :
346 : /*!
347 : * \brief Returns the time derivative of the lapse, which is 0
348 : */
349 1 : Scalar<DataType> get_var(
350 : const IntermediateComputer<DataType, Frame>& computer,
351 : ::Tags::dt<gr::Tags::Lapse<DataType>> /*meta*/);
352 :
353 : /*!
354 : * \brief Returns the spatial derivative of the lapse, which is 0
355 : */
356 1 : tnsr::i<DataType, 3, Frame> get_var(
357 : const IntermediateComputer<DataType, Frame>& computer,
358 : DerivLapse<DataType, Frame> /*meta*/);
359 :
360 : /*!
361 : * \brief Returns the shift, which is 0
362 : */
363 1 : tnsr::I<DataType, 3, Frame> get_var(
364 : const IntermediateComputer<DataType, Frame>& computer,
365 : gr::Tags::Shift<DataType, 3, Frame> /*meta*/);
366 :
367 : /*!
368 : * \brief Returns the time derivative of the shift, which is 0
369 : */
370 1 : tnsr::I<DataType, 3, Frame> get_var(
371 : const IntermediateComputer<DataType, Frame>& computer,
372 : ::Tags::dt<gr::Tags::Shift<DataType, 3, Frame>> /*meta*/);
373 :
374 : /*!
375 : * \brief Returns the spatial derivative of the shift, which is 0
376 : */
377 1 : tnsr::iJ<DataType, 3, Frame> get_var(
378 : const IntermediateComputer<DataType, Frame>& computer,
379 : DerivShift<DataType, Frame> /*meta*/);
380 :
381 : /*!
382 : * \brief Returns the time derivative of the spatial metric, which is 0
383 : */
384 1 : tnsr::ii<DataType, 3, Frame> get_var(
385 : const IntermediateComputer<DataType, Frame>& computer,
386 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3, Frame>> /*meta*/);
387 :
388 : /*!
389 : * \brief Computes and returns the square root of the determinant of the
390 : * spatial metric
391 : */
392 1 : Scalar<DataType> get_var(
393 : const IntermediateComputer<DataType, Frame>& computer,
394 : gr::Tags::SqrtDetSpatialMetric<DataType> /*meta*/);
395 :
396 : /*!
397 : * \brief Computes and returns the inverse spatial metric
398 : */
399 1 : tnsr::II<DataType, 3, Frame> get_var(
400 : const IntermediateComputer<DataType, Frame>& computer,
401 : gr::Tags::InverseSpatialMetric<DataType, 3, Frame> /*meta*/);
402 :
403 : /*!
404 : * \brief Computes and returns the extrinsic curvature which is 0
405 : */
406 1 : tnsr::ii<DataType, 3, Frame> get_var(
407 : const IntermediateComputer<DataType, Frame>& computer,
408 : gr::Tags::ExtrinsicCurvature<DataType, 3, Frame> /*meta*/);
409 : };
410 :
411 : private:
412 : /*!
413 : * \brief Mass of black hole A
414 : */
415 1 : double mass_a_{std::numeric_limits<double>::signaling_NaN()};
416 : /*!
417 : * \brief Mass of black hole B
418 : */
419 1 : double mass_b_{std::numeric_limits<double>::signaling_NaN()};
420 : /*!
421 : * \brief Center of black hole A
422 : */
423 1 : std::array<double, 3> center_a_ =
424 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
425 : /*!
426 : * \brief Center of black hole B
427 : */
428 1 : std::array<double, 3> center_b_ =
429 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
430 : };
431 :
432 : /*!
433 : * \brief Return whether two BrillLindquist data are equivalent
434 : */
435 1 : SPECTRE_ALWAYS_INLINE bool operator==(const BrillLindquist& lhs,
436 : const BrillLindquist& rhs);
437 :
438 : /*!
439 : * \brief Return whether two BrillLindquist data are not equivalent
440 : */
441 1 : SPECTRE_ALWAYS_INLINE bool operator!=(const BrillLindquist& lhs,
442 : const BrillLindquist& rhs);
443 : } // namespace gr::AnalyticData
|
