Binary.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <limits>
8 #include <optional>
9 
10 #include "DataStructures/CachedTempBuffer.hpp"
12 #include "DataStructures/TempBuffer.hpp"
13 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
15 #include "Elliptic/Systems/Xcts/Tags.hpp"
17 #include "Options/Auto.hpp"
18 #include "Options/Options.hpp"
20 #include "Parallel/PupStlCpp17.hpp"
21 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
22 #include "PointwiseFunctions/AnalyticData/Xcts/CommonVariables.hpp"
23 #include "PointwiseFunctions/AnalyticSolutions/Xcts/AnalyticSolution.hpp"
24 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
25 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
26 #include "Utilities/Requires.hpp"
27 #include "Utilities/TMPL.hpp"
28 #include "Utilities/TaggedTuple.hpp"
29 
30 /// \cond
31 namespace PUP {
32 class er;
33 } // namespace PUP
34 /// \endcond
35 
36 namespace Xcts::AnalyticData {
37 
38 namespace detail {
39 
40 template <typename DataType>
41 struct BinaryVariables;
42 
43 template <typename DataType>
44 using BinaryVariablesCache = cached_temp_buffer_from_typelist<
45  BinaryVariables<DataType>,
47  common_tags<DataType>,
48  Tags::ConformalMetric<DataType, 3, Frame::Inertial>,
50  tmpl::size_t<3>, Frame::Inertial>,
53  Tags::ShiftBackground<DataType, 3, Frame::Inertial>,
55  tmpl::size_t<3>, Frame::Inertial>,
56  Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
57  DataType, 3, Frame::Inertial>,
58  Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
59  Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
60  Tags::Conformal<gr::Tags::MomentumDensity<3, Frame::Inertial, DataType>,
61  0>,
62  // For initial guesses
63  Tags::ConformalFactor<DataType>,
64  Tags::LapseTimesConformalFactor<DataType>,
65  Tags::ShiftExcess<DataType, 3, Frame::Inertial>>>;
66 
67 template <typename DataType>
68 struct BinaryVariables
69  : CommonVariables<DataType, BinaryVariablesCache<DataType>> {
70  static constexpr size_t Dim = 3;
71  using Cache = BinaryVariablesCache<DataType>;
72  using CommonVariables<DataType, BinaryVariablesCache<DataType>>::operator();
73 
74  using superposed_tags = tmpl::list<
75  Tags::ConformalMetric<DataType, Dim, Frame::Inertial>,
77  tmpl::size_t<Dim>, Frame::Inertial>,
80  Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
81  Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
82  Tags::Conformal<gr::Tags::MomentumDensity<Dim, Frame::Inertial, DataType>,
83  0>,
84  Tags::ConformalFactor<DataType>,
85  Tags::LapseTimesConformalFactor<DataType>,
86  Tags::ShiftExcess<DataType, Dim, Frame::Inertial>>;
87 
88  const tnsr::I<DataVector, Dim>& x;
89  const double angular_velocity;
90  const std::optional<std::array<double, 2>> falloff_widths;
91  const std::array<tnsr::I<DataVector, Dim>, 2> x_isolated;
92  const std::array<DataVector, 2> windows;
93  const tuples::tagged_tuple_from_typelist<superposed_tags> flat_vars;
95  isolated_vars;
96 
97  template <typename Tag,
99  void operator()(gsl::not_null<typename Tag::type*> superposed_var,
100  gsl::not_null<Cache*> /*cache*/,
101  Tag /*meta*/) const noexcept {
102  for (size_t i = 0; i < superposed_var->size(); ++i) {
103  (*superposed_var)[i] =
104  get<Tag>(flat_vars)[i] +
105  windows[0] *
106  (get<Tag>(isolated_vars[0])[i] - get<Tag>(flat_vars)[i]) +
107  windows[1] * (get<Tag>(isolated_vars[1])[i] - get<Tag>(flat_vars)[i]);
108  }
109  if constexpr (std::is_same_v<
110  Tag, ::Tags::deriv<Tags::ConformalMetric<DataType, 3,
112  tmpl::size_t<3>, Frame::Inertial>>) {
113  add_deriv_of_window_function(superposed_var);
114  }
115  }
116  void operator()(
117  gsl::not_null<tnsr::I<DataType, Dim>*> shift_background,
119  Tags::ShiftBackground<DataType, Dim, Frame::Inertial> /*meta*/)
120  const noexcept;
121  void operator()(
122  gsl::not_null<tnsr::iJ<DataType, Dim>*> deriv_shift_background,
124  ::Tags::deriv<Tags::ShiftBackground<DataType, Dim, Frame::Inertial>,
125  tmpl::size_t<Dim>, Frame::Inertial> /*meta*/)
126  const noexcept;
127  void operator()(gsl::not_null<tnsr::II<DataType, Dim, Frame::Inertial>*>
128  longitudinal_shift_background_minus_dt_conformal_metric,
129  gsl::not_null<Cache*> cache,
130  Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
131  DataType, Dim, Frame::Inertial> /*meta*/) const noexcept;
132 
133  private:
134  void add_deriv_of_window_function(gsl::not_null<tnsr::ijj<DataType, Dim>*>
135  deriv_conformal_metric) const noexcept;
136 };
137 } // namespace detail
138 
139 /// \cond
140 template <typename IsolatedObjectRegistrars, typename Registrars>
141 class Binary;
142 
143 namespace Registrars {
144 template <typename IsolatedObjectRegistrars>
145 struct Binary {
146  template <typename Registrars>
148 };
149 } // namespace Registrars
150 /// \endcond
151 
152 /*!
153  * \brief Binary compact-object data in general relativity, constructed from
154  * superpositions of two isolated objects.
155  *
156  * This class implements background data for the XCTS equations describing two
157  * objects in a quasi-equilibrium orbit, i.e. with \f$\bar{u}=0\f$ and
158  * \f$\partial_t K=0\f$. Both objects can be chosen from the list of
159  * `IsolatedObjectRegistrars`, e.g. they can be black-hole or neutron-star
160  * solutions in different coordinates. Most quantities are constructed by
161  * superposing the two isolated solutions (see e.g. Eq. (8-9) in
162  * \cite Varma2018sqd or Eq. (45-46) in \cite Lovelace2008tw):
163  *
164  * \f{align}
165  * \bar{\gamma}_{ij} &= f_{ij} + \sum_{\alpha=1}^2
166  * e^{-r_\alpha^2 / w_\alpha^2}\left(\gamma^\alpha_{ij} - f_{ij}\right) \\
167  * K &= \sum_{\alpha=1}^2 e^{-r_\alpha^2 / w_\alpha^2}K^\alpha
168  * \f}
169  *
170  * where \f$\gamma^\alpha_{ij}\f$ and \f$K^\alpha\f$ denote the spatial metric
171  * and extrinsic-curvature trace of the two individual solutions, \f$r_\alpha\f$
172  * is the Euclidean coordinate-distance from the center of each object and
173  * \f$w_\alpha\f$ are parameters describing Gaussian falloff-widths. The
174  * Gaussian falloffs facilitate that the influence of either of the two objects
175  * at the position of the other is strongly damped, and they also avoid
176  * logarithmic scaling of the solution at large distances where we would
177  * typically employ an inverse-radial coordinate map and asymptotically-flat
178  * boundary conditions. The falloff-widths are chosen in terms of the Newtonian
179  * Lagrange points of the two objects in \cite Varma2018sqd and
180  * \cite Lovelace2008tw, and they are input parameters in this implementation.
181  * The falloff can be disabled by passing `std::nullopt` to the constructor, or
182  * `None` in the input file.
183  *
184  * The remaining quantities that this class implements relate to the orbital
185  * motion of the two objects. To obtain initial data in "co-rotating"
186  * coordinates where the two objects are initially at rest we prescribe the
187  * background shift
188  *
189  * \f{equation} \beta^i_\mathrm{background} = (-\Omega y, \Omega x, 0) \f}
190  *
191  * where \f$\Omega\f$ is the angular-velocity parameter.
192  */
193 template <typename IsolatedObjectRegistrars,
194  typename Registrars = tmpl::list<
195  Xcts::AnalyticData::Registrars::Binary<IsolatedObjectRegistrars>>>
196 class Binary : public ::AnalyticData<3, Registrars> {
197  private:
199 
200  public:
201  using IsolatedObjectBase =
203 
204  struct XCoords {
205  static constexpr Options::String help =
206  "The coordinates on the x-axis where the two objects are placed";
207  using type = std::array<double, 2>;
208  };
209  struct ObjectA {
210  static constexpr Options::String help =
211  "The object placed on the negative x-axis";
213  };
214  struct ObjectB {
215  static constexpr Options::String help =
216  "The object placed on the positive x-axis";
218  };
220  static constexpr Options::String help = "Orbital angular velocity";
221  using type = double;
222  };
223  struct FalloffWidths {
224  static constexpr Options::String help =
225  "The widths for the window functions around the two objects, or 'None' "
226  "to disable the Gaussian falloff.";
228  };
229  using options =
230  tmpl::list<XCoords, ObjectA, ObjectB, AngularVelocity, FalloffWidths>;
231  static constexpr Options::String help =
232  "Binary compact-object data in general relativity, constructed from "
233  "superpositions of two isolated objects.";
234 
235  Binary() = default;
236  Binary(const Binary&) = delete;
237  Binary& operator=(const Binary&) = delete;
238  Binary(Binary&&) = default;
239  Binary& operator=(Binary&&) = default;
240  ~Binary() = default;
241 
244  std::unique_ptr<IsolatedObjectBase> object_b, double angular_velocity,
245  std::optional<std::array<double, 2>> falloff_widths) noexcept
246  : xcoords_(xcoords),
247  superposed_objects_({std::move(object_a), std::move(object_b)}),
248  angular_velocity_(angular_velocity),
249  falloff_widths_(falloff_widths) {}
250 
251  explicit Binary(CkMigrateMessage* m) noexcept : Base(m) {}
252  using PUP::able::register_constructor;
254 
255  template <typename DataType, typename... RequestedTags>
256  tuples::TaggedTuple<RequestedTags...> variables(
257  const tnsr::I<DataType, 3, Frame::Inertial>& x,
258  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
259  return variables_impl<DataType>(x, std::nullopt, std::nullopt,
260  tmpl::list<RequestedTags...>{});
261  }
262  template <typename... RequestedTags>
263  tuples::TaggedTuple<RequestedTags...> variables(
264  const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
265  const InverseJacobian<DataVector, 3, Frame::Logical, Frame::Inertial>&
266  inv_jacobian,
267  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
268  return variables_impl<DataVector>(x, mesh, inv_jacobian,
269  tmpl::list<RequestedTags...>{});
270  }
271 
272  // NOLINTNEXTLINE
273  void pup(PUP::er& p) noexcept override {
274  Base::pup(p);
275  p | xcoords_;
276  p | superposed_objects_;
277  p | angular_velocity_;
278  p | falloff_widths_;
279  }
280 
281  const std::array<double, 2>& x_coords() const noexcept { return xcoords_; }
282  const std::array<std::unique_ptr<IsolatedObjectBase>, 2>& superposed_objects()
283  const noexcept {
284  return superposed_objects_;
285  }
286  double angular_velocity() const noexcept { return angular_velocity_; }
287  const std::optional<std::array<double, 2>>& falloff_widths() const noexcept {
288  return falloff_widths_;
289  }
290 
291  private:
292  std::array<double, 2> xcoords_{};
293  std::array<std::unique_ptr<IsolatedObjectBase>, 2> superposed_objects_{};
294  Xcts::Solutions::Flatness<> flatness_{};
295  double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
296  std::optional<std::array<double, 2>> falloff_widths_{};
297 
298  template <typename DataType, typename... RequestedTags>
299  tuples::TaggedTuple<RequestedTags...> variables_impl(
300  const tnsr::I<DataType, 3, Frame::Inertial>& x,
303  const InverseJacobian<DataType, 3, Frame::Logical, Frame::Inertial>>>
304  inv_jacobian,
305  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
306  std::array<tnsr::I<DataVector, 3>, 2> x_isolated{{x, x}};
307  std::array<DataVector, 2> euclidean_distance{};
308  std::array<DataVector, 2> windows{};
309  // Possible optimization: Only retrieve those superposed tags from the
310  // isolated solutions that are actually needed. This needs some dependency
311  // logic, because some of the non-superposed tags depend on superposed tags.
312  using VarsComputer = detail::BinaryVariables<DataType>;
313  using requested_superposed_tags = typename VarsComputer::superposed_tags;
315  isolated_vars;
316  for (size_t i = 0; i < 2; ++i) {
317  get<0>(gsl::at(x_isolated, i)) -= gsl::at(xcoords_, i);
318  gsl::at(euclidean_distance, i) = get(magnitude(gsl::at(x_isolated, i)));
319  if (falloff_widths_.has_value()) {
320  gsl::at(windows, i) = exp(-square(gsl::at(euclidean_distance, i)) /
321  square(gsl::at(*falloff_widths_, i)));
322  } else {
323  gsl::at(windows, i) = make_with_value<DataVector>(x, 1.);
324  }
325  gsl::at(isolated_vars, i) =
326  gsl::at(superposed_objects_, i)
327  ->variables(gsl::at(x_isolated, i), requested_superposed_tags{});
328  }
329  auto flat_vars = flatness_.variables(x, requested_superposed_tags{});
330  typename VarsComputer::Cache cache{
331  get_size(*x.begin()),
332  VarsComputer{{std::move(mesh), std::move(inv_jacobian)},
333  x,
334  angular_velocity_,
335  falloff_widths_,
336  std::move(x_isolated),
337  std::move(windows),
338  std::move(flat_vars),
339  std::move(isolated_vars)}};
340  return {cache.get_var(RequestedTags{})...};
341  }
342 };
343 
344 /// \cond
345 template <typename IsolatedObjectRegistrars, typename Registrars>
346 PUP::able::PUP_ID Binary<IsolatedObjectRegistrars, Registrars>::my_PUP_ID =
347  0; // NOLINT
348 /// \endcond
349 
350 } // namespace Xcts::AnalyticData
domain::push_front
CoordinateMap< SourceFrame, TargetFrame, NewMap, Maps... > push_front(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by prepending the new map to the beginning of the old maps.
gr::Tags::TraceExtrinsicCurvature
Definition: Tags.hpp:120
Xcts::AnalyticData::Binary::XCoords
Definition: Binary.hpp:204
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
Xcts::AnalyticData::Binary::AngularVelocity
Definition: Binary.hpp:219
CharmPupable.hpp
Frame::Inertial
Definition: IndexType.hpp:44
Xcts::AnalyticData
Analytic data for the XCTS equations, i.e. field configurations that do not solve the equations but a...
Definition: AnalyticData.hpp:8
Options.hpp
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
Xcts::Solutions::Flatness
Flat spacetime in general relativity. Useful as initial guess.
Definition: Flatness.hpp:44
Xcts::AnalyticData::Binary
Binary compact-object data in general relativity, constructed from superpositions of two isolated obj...
Definition: Binary.hpp:196
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Xcts::AnalyticData::Binary::ObjectA
Definition: Binary.hpp:209
get_size
decltype(auto) get_size(const T &t, SizeFunction size=GetContainerSize{}) noexcept
Retrieve the size of t if t.size() is a valid expression, otherwise if T is fundamental or a std::com...
Definition: ContainerHelpers.hpp:145
Xcts::Tags::ConformalMetric
Conformal< gr::Tags::SpatialMetric< Dim, Frame, DataType >, -4 > ConformalMetric
The conformally scaled spatial metric , where is the Xcts::Tags::ConformalFactor and is the gr::Tag...
Definition: Tags.hpp:44
Options::AutoLabel::None
'None' label
Definition: Auto.hpp:20
std::reference_wrapper
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
array
Xcts::AnalyticData::Binary::ObjectB
Definition: Binary.hpp:214
Xcts::Solutions::AnalyticSolution
Base class for analytic solutions of the XCTS equations.
Definition: AnalyticSolution.hpp:34
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
Options::Auto
A class indicating that a parsed value can be automatically computed instead of specified.
Definition: Auto.hpp:36
AnalyticData
Provides analytic tensor data as a function of the spatial coordinates.
Definition: AnalyticData.hpp:34
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:29
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
limits
std::exp
T exp(T... args)
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
Tensor.hpp
Xcts::AnalyticData::Binary::FalloffWidths
Definition: Binary.hpp:223
Requires.hpp
magnitude
Scalar< DataType > magnitude(const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector) noexcept
Compute the Euclidean magnitude of a rank-1 tensor.
Definition: Magnitude.hpp:27
Tags::deriv
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:52
optional
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
PartialDerivatives.hpp
std::unique_ptr
Prefixes.hpp
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13